Addressing Neutrino Mixing Models with DUNE and T2HK

We consider schemes of neutrino mixing arising within the discrete symmetry approach to the well-known flavour problem. We concentrate on $3\nu$ mixing schemes in which the cosine of the Dirac CP violation phase $\delta_\mathrm{CP}$ satisfies a sum rule by which it is expressed in terms of three neutrino mixing angles $\theta_{12}$, $\theta_{23}$, and $\theta_{13}$, and a fixed real angle $\theta^\nu_{12}$, whose value depends on the employed discrete symmetry and its breaking. We consider five underlying symmetry forms of the neutrino mixing matrix: bimaximal (BM), tri-bimaximal (TBM), golden ratio A (GRA) and B (GRB), and hexagonal (HG). For each symmetry form, the sum rule yields specific prediction for $\cos\delta_\mathrm{CP}$ for fixed $\theta_{12}$, $\theta_{23}$, and $\theta_{13}$. In the context of the proposed DUNE and T2HK facilities, we study (i) the compatibility of these predictions with present neutrino oscillation data, and (ii) the potential of these experiments to discriminate between various symmetry forms.


Introduction and Motivation
All compelling neutrino oscillation data are compatible with 3ν mixing [1], i.e., with existence of three light neutrino states ν 1,2,3 with definite masses m 1,2,3 , three orthogonal linear combinations of which form the three flavour neutrino states ν e , ν µ and ν τ . The flavour neutrino (flavour antineutrino) oscillation probabilities in the case of 3ν mixing are characterised, as is well known, by six fundamental parameters. In the standard parametrisation of the Pontecorvo, Maki, Nakagawa, Sakata (PMNS) neutrino mixing matrix [1], these are the solar, reactor, and atmospheric mixing angles, θ 12 , θ 13 , and θ 23 , respectively, the Dirac CP violation (CPV) phase, δ CP , and the two independent mass squared differences, ∆m 2 21 ≡ m 2 2 − m 2 1 and ∆m 2 31 ≡ m 2 3 − m 2 1 . The mixing angles θ 12 and θ 13 as well as the solar and the absolute value of the atmospheric neutrino mass squared differences, ∆m 2 21 and |∆m 2 31 |, respectively, have been measured in neutrino oscillation experiments with a relatively high precision [2][3][4]. The precision on θ 23 is somewhat worse, the relative 1σ uncertainty on sin 2 θ 23 being approximately 10%. At the same time, the octant of θ 23 and the value of δ CP remain unknown. The sign of ∆m 2 31 , which is also undetermined, allows, as is well known, for two possible types of the neutrino mass spectrum: (i) with normal ordering (NO) if ∆m 2 31 > 0, and (ii) with inverted ordering (IO) if ∆m 2 31 < 0. In Table 1, we summarise the best fit values and 1σ, 2σ, 3σ ranges of the oscillation parameters obtained in one of the latest global analysis of the neutrino oscillation data [3] 1 . In the case of IO spectrum, instead of ∆m 2 31 , we use the largest positive mass squared difference ∆m 2 23 ≡ m 2 2 − m 2 3 , which corresponds to ∆m 2 31 of the NO spectrum. In particular, we see from Table 1 that for the NO mass spectrum, the values of δ CP ∈ (31 • , 137 • ) are already disfavoured at more than 3σ C.L., while the values of δ CP between 180 • and 342 • are allowed at 2σ.
Understanding the origin of the patterns of neutrino oscillation parameters revealed by the data is one of the most challenging problems in neutrino physics. It is a part of the more general fundamental problem in particle physics of understanding the origins of flavour, i.e., the patterns of quark, charged lepton, and neutrino masses, and quark and lepton mixing. There exists a possibility that the high-precision measurements of the oscillation parameters may shed light on the origin of the observed pattern of neutrino mixing and lepton flavour. This would be the case if the observed form of neutrino (and possibly quark) mixing were determined by an underlying discrete flavour symmetry. One of the most striking features of the discrete symmetry approach to neutrino mixing and lepton flavour (see, e.g., [5][6][7] for reviews), is that it leads to (i) fixed predictions of the values of some of the neutrino mixing angles and the Dirac CPV phase δ CP , and/or (ii) existence of correlations between some of the mixing angles and/or between the mixing angles and δ CP (see, e.g., [8][9][10][11][12][13][14][15][16][17][18][19]). These correlations are often referred to as neutrino mixing sum rules 2 . Most importantly, these sum rules can be tested using oscillation data [8-11, 15, 19, 29-31].
Within the discrete flavour symmetry approach, the PMNS matrix is predicted to have an underlying symmetry form, where θ 12 , θ 23 , and θ 13 Table 1: The best fit values and 1σ, 2σ, 3σ ranges of the neutrino oscillation parameters obtained in the global analysis of the neutrino oscillation data performed in [3]. NO (IO) stands for normal (inverted) ordering of the neutrino mass spectrum.
perturbative corrections from their respective measured values. The approach seems very natural in view of the fact that U PMNS = U † e U ν , where U e and U ν are 3 × 3 unitary matrices which diagonalise the charged lepton and neutrino mass matrices. Typically (but not universally) the matrix U ν has a certain symmetry form, while the matrix U e provides the corrections necessary to bring the symmetry values of the angles in U ν to their experimentally measured values. A sum rule which relates cos δ CP with θ 12 , θ 23 , and θ 13 , arising in this approach, depends on the underlying symmetry form of the PMNS matrix and on the form of the "correcting" matrix U e .
In the case of the BM symmetry form, the obtained best fit value of cos δ CP = −1.26 is unphysical. This reflects the fact that the BM symmetry form does not provide a good description of the present best fit values of the neutrino mixing angles, as discussed in [13]. As can be seen from Table 3, current uncertainties on the mixing angles allow us to accommodate physical values of cos δ CP for the BM symmetry form. For instance, fixing sin 2 θ 13 and sin 2 θ 23 to their best fit values, cos δ CP = −1 requires sin 2 θ 12 = 0.3343, which is the upper bound of the corresponding 2σ allowed range of sin 2 θ 12 (see Table 1).
A rather detailed analysis of the predictions for cos δ CP of the sum rule in eq. (1.1) has been performed in Refs. [9,10]. In particular, likelihood profiles for cos δ CP for each symmetry form have been presented using the current and prospective precision on the neutrino mixing parameters (see Figs. 12 and 13 in [9]). In the present work, using the potential of the future long-baseline (LBL) neutrino oscillation experiments 3 , namely, Deep Underground Neutrino Experiment (DUNE) and Tokai to Hyper-Kamiokande (T2HK), we study in detail (i) to what degree the sum rule predictions for cos δ CP are compatible with the present neutrino oscillation data, and (ii) how well the considered symmetry forms, BM, TBM, GRA, GRB, and HG, can be discriminated from each other.
The layout of the article is as follows. In Section 2, we take a first glance at the sum rule predictions. In Section 3, we give a short description of the planned DUNE and T2HK experiments and provide expected event rates for both the set-ups. Section 4 contains details of the statistical analysis. In Section 5, we present and discuss results of this analysis. More specifically, in subsection 5.1, we test the compatibility between the considered symmetry forms and present oscillation data. Next, in subsection 5.2, we explore the potential of DUNE, T2HK, and their combination to distinguish between the symmetry forms in question under the assumption that one of them is realised in Nature. In subsection 5.3, we consider the BM symmetry form using the values of the mixing angles for which this form is viable, and study at which C.L. it can be distinguished from the other symmetry forms considered. We conclude in Section 6. Appendix A discusses the issue of external priors on sin 2 θ 12 and sin 2 θ 13 . In Appendix B, we show the impact of marginalisation over ∆m 2 31 . Finally, in Appendix C, we study the compatibility of the considered symmetry forms with any potentially true values of sin 2 θ 23 and δ CP in the context of DUNE and T2HK.

A First Glance at the Sum Rule Predictions
Let us first take a closer look at the sum rule predictions for δ CP summarised in Tables 2 and  3. Several important points that can be made from these two tables and the current global data (see Table 1) are in order.
• The parameter θ ν 12 has a fixed value for each symmetry form as given in Table 2. For fixed choices of θ 12 , θ 23 , θ 13 , and θ ν 12 , eq. (1.1) predicts a certain value of cos δ CP , which gives rise to two values of δ CP (see fourth column of Table 2).
• The complementarity [52] between the modern reactor (Daya Bay, RENO, and Double Chooz) and accelerator (T2K and NOνA) data has enabled us to probe the parameter space for δ CP . Already, the latest global data have disfavoured values of δ CP ∈ (31 • , 137 • ) at 3σ C.L. for NO (see Table 1). From Table 3, we see that out of the three mixing angles, the 3σ allowed range of θ 12 causes the largest uncertainty in δ CP predicted by eq. (1.1). But, note that, all the symmetry forms except BM predict one of the ranges of δ CP in the interval of 31 • to 137 • , which has been already ruled out at 3σ. Therefore, we will not consider these ranges further in our study, and we will only consider the values of δ CP in the interval of 180 • to 360 • .
• Further, we notice from Tables 1 and 3 that for TBM and GRB, the predicted intervals of δ CP lie within the 1σ experimentally allowed range. For BM, GRA, and HG, the intervals of interest fall within the 2σ range. Now, it would be quite interesting to assess the sensitivity of the future LBL experiments in discriminating various symmetry forms, which is the main thrust of the present work.
• T2HK and DUNE will not be able to constrain θ 12 , which causes the largest uncertainty in predicting the range of δ CP (see Table 3). However, the proposed JUNO experiment will provide a high-precision measurement of sin 2 θ 12 with a relative 1σ uncertainty of 0.7%. Therefore, we impose a prior on sin 2 θ 12 expected from JUNO, which we will discuss later in detail in Section 4 and in Appendix A.
• Table 3 shows that the sum rule predictions depend to some extent on θ 23 . Therefore, we vary this angle both in data and in fit. The LBL experiments themselves are sensitive to θ 23 . Thus, we do not impose any external prior on this angle (see details in Section 4).
• The experiments under discussion are sensitive to θ 13 through the appearance channels ν µ → ν e andν µ →ν e . Therefore, the role of an external prior on sin 2 θ 13 is negligible for the physics case under study, as we will see later in Fig. 4. However, we put a prior on sin 2 θ 13 as expected from Daya Bay to speed up our simulations (see Section 4 for details).

Experimental Features and Event Rates
In this section, we first briefly describe the key experimental features of the proposed highprecision DUNE and T2HK facilities that we use in our numerical simulation.

The Next Generation Experiments: DUNE and T2HK
The planned Deep Underground Neutrino Experiment (DUNE) aims to achieve new milestones in the intensity frontier with a new, high-intensity, on-axis, wide-band neutrino beam from Fermilab directed towards a massive liquid argon time-projection chamber (LArTPC) far detector housed at the Homestake Mine in South Dakota over a baseline of 1300 km [53][54][55][56][57].
In our simulation, we consider a fiducial mass of 35 kt for the far detector and the detector characteristics have been taken from Table 1 of Ref. [58]. As far as beam specifications are concerned, we assume a modest proton beam power of 708 kW in its initial phase with 120 GeV proton energy, which can supply 6 × 10 20 protons on target (p.o.t.) in 188 days per calendar year. In our calculation, we have used the fluxes which were generated assuming a decay pipe length of 200 m and 200 kA horn current [59]. We assume that DUNE will collect data for ten years (5 years in ν mode and 5 years inν mode), which is equivalent to a total exposure of 248 kt · MW · year 4 . In our simulation, we consider the reconstructed ν energy range to be 0.5 GeV to 10 GeV for both appearance and disappearance channels. We take the same energy range for antineutrino as well.
The proposed Hyper-Kamiokande (HK) water Cherenkov detector will serve as the far detector of a long-baseline neutrino experiment using an upgraded neutrino beam from the J-PARC facility, commonly known as "T2HK" (Tokai to Hyper-Kamiokande) experiment [60][61][62]. This set-up is highly sensitive to the Dirac CPV phase δ CP of the PMNS 3ν mixing matrix and holds promise to resolve the mystery of leptonic CP violation in neutrino oscillations at an unprecedented confidence level [61]. We perform the simulation for T2HK according to Refs. [61,62]. To produce an intense ν/ν beam for HK, we consider an integrated proton beam power of 7.5 MW × 10 7 seconds, which can deliver in total 15.6 × 10 21 p.o.t. with a 30 GeV proton beam. We assume that these total p.o.t. will be shared among ν andν modes with a run-time ratio of 1:3 to have almost equal statistics in both the modes. The huge 560 kt (fiducial) HK detector will be placed in the Tochibora mine, at a distance of 295 km from J-PARC at an off-axis angle of ∼ 2.5 • , which will produce a narrow band beam with a sharp peak around the first oscillation maximum of 0.6 GeV. The total exposure that we consider for T2HK is 4200 kt · MW · year. In our simulation, we take the reconstructed ν e andν e energy range of 0.1 GeV to 1.25 GeV for the appearance channel. As far as the disappearance channel is concerned, the assumed energy range is 0.1 GeV to 7 GeV for both ν µ andν µ candidate events.
Recently, the baseline design for T2HK has been revised [63]. According to this latest publication [63], the total beam exposure is 27 × 10 21 p.o.t. and the HK design proposes the construction of two identical water Cherenkov detectors in stage with fiducial mass of 187 kt per detector. The possibility of placing the first detector near the Super-Kamiokande site, 295 km away and 2.5 • off-axis from the J-PARC neutrino beam and the second detector in Korea having a baseline of 1100 km from J-PARC at an off-axis angle of ∼ 2.5 • has also been explored in Ref. [63], and this set-up has been referred as "T2HKK". We follow the details as given in Ref. [63] to simulate the T2HKK set-up.

Event Rates
In this subsection, we present the expected total event rates for both the set-ups under consideration. We compute the number of expected electron events 5 in the i-th energy bin in the detector using the following well-known expression: , T is the total running time (year), n n is the number of target nucleons in the detector, is the detector efficiency, and R(E, E A ) is the Gaußian energy resolution function (GeV −1 ) of the detector. σ νe is the neutrino interaction cross section (m 2 ) which has been taken from Refs. [64,65], where the authors calculated the cross section for water and isoscalar targets. In order to have LAr cross sections for DUNE, we have scaled the inclusive charged current (CC) cross sections of water by a factor of 1.06 (0.94) for neutrino (antineutrino) [66,67]. We denote the true and reconstructed neutrino energies by the quantities E and E A , respectively. Table 4 shows a comparison between the expected total signal and background event rates 6 Mode (Channel)  Table 1). Here "Int" means intrinsic beam contamination, "Mis-id" represents mis-identified muon events, and "NC" stands for neutral current. See text for other details.
in the appearance and disappearance modes for DUNE and T2HK set-ups. We compute the same for both neutrino and antineutrino runs assuming a total exposure of 248 kt · MW · year for DUNE and 4200 kt · MW · year for T2HK. We consider a run-time ratio of 1:1 among neutrino and antineutrino modes in DUNE and the corresponding ratio is 1:3 in T2HK. The total event rates are calculated assuming NO, δ CP = 248 • , and sin 2 θ 23 = 0.425. For all other oscillation parameters, we consider the best fit values which are applicable for NO (see Table 1). To compute the full three-flavour neutrino oscillation probabilities in matter, we take the line-averaged constant Earth matter density of 2.80 (2.87) g/cm 3 for the T2HK (DUNE) baseline following the Preliminary Reference Earth Model (PREM) [68].
The main sources of backgrounds while selecting the ν e andν e candidate events are the intrinsic ν e /ν e component in the beam, the muon events which will be mis-identified as electron events, and the neutral current (NC) events. Table 4 clearly depicts that in case of appearance searches, the dominant background component is the intrinsic ν e /ν e in the beam. For the ν µ /ν µ candidate events, the main backgrounds are the NC events. Though we present the total event rates in Table 4, but, in our simulation, we have performed a full spectral analysis using the binned events spectra for both the DUNE and T2HK set-ups.

Details of Statistical Analysis
This section is devoted to describe the strategy that we adopt for the statistical treatment to quantify the sensitivities of DUNE and T2HK in testing various lepton mixing schemes. To produce our results, we take the help of the widely used GLoBES software [69,70] which calculates the median sensitivity of the experiment without considering the statistical fluctuations. Unless stated otherwise, we generate our simulated data considering the best fit values of the oscillation parameters obtained in the global analysis assuming NO for the neutrino mass spectrum (see second column of Table 1). We also keep the choice of the neutrino mass ordering to be fixed to NO in the fit 7 . The solar and atmospheric mass squared differences are already very well measured [2][3][4] and moreover, they do not appear in the sum rule (see eq. (1.1)) that relates cos δ CP with the mixing angles. Therefore, we also keep them fixed in the fit at their best fit values while showing our results. Only in Appendix B, we give a plot, where we marginalise over test ∆m 2 31 in the fit in its present 3σ allowed range of (2.45 − 2.69) × 10 −3 eV 2 . The mixing angles play an important role in the sum rule and we treat them very carefully in our analysis. In the fit, we marginalise over test sin 2 θ 23 in its present 3σ allowed range of 0.381 to 0.615. We show few results where we also vary the true value of sin 2 θ 23 or marginalise over it in the same 3σ range. We do not impose any external prior on sin 2 θ 23 as it will be directly measured by the DUNE and T2HK experiments. The sum rule as given in eq. (1.1) is very sensitive to the value of sin 2 θ 12 . We vary this parameter in the fit in its present 3σ allowed range of 0.25 to 0.354. Since both DUNE and T2HK cannot constrain the solar mixing angle (see the probability expressions in [74]), we impose an external Gaußian prior of 0.7% (at 1σ) on this parameter as the proposed medium-baseline reactor oscillation experiment JUNO will be able to measure sin 2 θ 12 with this precision [75]. We also marginalise over test sin 2 θ 13 in its present 3σ allowed range of 0.019 to 0.024. While doing so, we apply an external Gaußian prior of 3% (at 1σ) on this parameter expecting that the Daya Bay experiment would be able to achieve this precision by the end of its run [76]. Both DUNE and T2HK can measure θ 13 with high precision using ν µ → ν e andν µ →ν e oscillation channels and therefore, the prior on θ 13 is not very crucial in our study (see Fig. 4 and related discussion in Appendix A). Still, we use this prior in our analysis to speed up the marginalisation procedure. For a given test choice of θ 12 , θ 13 , and θ 23 , the test value of δ CP is calculated using the sum rule (see eq. (1.1)) for a particular choice of the lepton mixing scheme which is characterised by a certain value of θ ν 12 . Since the best fit value of δ CP may change in the future, we also show some results varying the true choice of δ CP in the range 180 • to 360 • .
To perform our statistical analysis, we follow the procedure outlined in Refs. [77,78], and use the following Poissonian χ 2 function: where n is the total number of reconstructed energy bins and Above, N th i ({ω}) denotes the predicted number of CC signal events (estimated using eq. (3.1)) in the i-th energy bin for a set of oscillation parameters ω. N b i ({ω}) stands for the total number of background events 8 . The quantities π s and π b in eq. (4.2) represent the systematic errors on the signal and background events, respectively. For DUNE (T2HK), we consider π s = 5% (5%) and π b = 5% (10%) in the form of normalisation errors for both the appearance and disappearance channels. We take the same uncorrelated systematic uncertainties for both the neutrino and antineutrino modes. The quantities ξ s and ξ b denote the "pulls" due atmospheric neutrinos at more than 3σ C.L. for both NO and IO provided sin 2 θ23 > 0.45 [73]. Combining beam and atmospheric neutrinos in HK, the mass ordering can be determined at more than 3σ (5σ) with five (ten) years of data [73]. 8 Note that we consider both CC and NC background events in our analysis and the NC backgrounds are independent of oscillation parameters.
to the systematic error on the signal and background, respectively. The data in eq. (4.1) are included through the variable x i = N ex i + N b i , where N ex i is the number of observed CC signal events and N b i is the background as discussed earlier. To obtain the total χ 2 , we add the χ 2 contributions coming from all the relevant oscillation channels for both neutrino and antineutrino modes in a given experiment in the following fashion: In the above expression, we assume that all the oscillation channels for both neutrino and antineutrino modes are completely uncorrelated, all the energy bins in a given channel are fully correlated, and the systematic errors on signal and background are fully uncorrelated. The fact that the flux normalisation errors in ν µ → ν e and ν µ → ν µ oscillation channels are same (i.e., they are correlated) is taken into account in the error budget for each of the two channels. However, there are other uncertainties which contribute to the total normalisation error for each of the two channels, like the uncertainties in cross sections, detector efficiencies, etc., which are uncorrelated. For this reason, we simply assume that the total normalisation errors in these two channels are uncorrelated. The same is true forν µ →ν e andν µ →ν µ oscillation channels. In our opinion, with the current understanding of the two detectors in DUNE and T2HK experiments, it is premature to perform a very detailed analysis taking into account such fine effects as, e.g., the correlation between the flux normalisation errors in the appearance and disappearance channels. In eq. (4.3), the last term appears due to the external Gaußian priors that we impose on sin 2 θ 12 and sin 2 θ 13 in the following way: where we take σ(sin 2 θ 12 ) = 0.007 × sin 2 θ true 12 and σ(sin 2 θ 13 ) = 0.03 × sin 2 θ true 13 as mentioned earlier. In our analysis, we assume that the true values of sin 2 θ 12 and sin 2 θ 13 will remain unchanged in the future. While implementing the minimisation procedure, χ 2 total is first minimised with respect to the "pull" variables ξ s , ξ b , and then marginalised over the various oscillation parameters in their allowed ranges in the fit as discussed above to obtain ∆χ 2 min . In Fig. 6 in Appendix C, we quote the statistical significance of our results for 1 d.o.f. in terms of nσ, where n = ∆χ 2 min , which is valid in the frequentist method of hypothesis testing [79]. Table 1 shows the current best fit values of the neutrino mixing angles and the CPV phase δ CP . Equation (1.1) relates δ CP with the parameter θ ν 12 which characterises the symmetry forms under consideration. The first question we want to answer is how much compatibility there is between the mixing symmetry forms and the present best fit values of the oscillation parameters. To this aim, we assume the current best fit values of θ 12 , θ 13 , θ 23 , and δ CP to be their true values and generate prospective data using the DUNE, T2HK, and T2HKK experimental set-ups according to the procedure explained in Section 4. Further, in the test, we fix the mixing angles to their best fit values and let δ CP to vary between 180 • and 360 • . For each value of δ CP , first we calculate sin 2 θ ν 12 according to eq. (1.1), and then, we estimate the corresponding ∆χ 2 using the prospective data from combined DUNE + T2HK and DUNE + T2HKK set-ups. This is ∆χ 2 between a given value of δ CP and its best fit value, i.e., δ CP = 248 • . In Fig. 1, we plot ∆χ 2 as a function of sin 2 θ ν 12 and show on the resulting curves:

Compatibility between Various Symmetry Forms and Present Neutrino Oscillation Data
• the black dot corresponding to the current best fit value of δ CP = 248 • which translates to sin 2 θ ν 12 = 0.364 (∆χ 2 = 0); • the coloured dots corresponding to the values of sin 2 θ ν 12 which characterise the GRB (violet), TBM (red), GRA (blue) and HG (green) symmetry forms.
From this figure we see that, if the present best fit values of θ 12 , θ 13 , θ 23 , and δ CP were the true values of these parameters, the GRB (TBM) symmetry form would be compatible with them at slightly less (more) than 1σ C.L., while the GRA and HG schemes would be disfavoured at more than 2.7σ and 3.7σ, respectively, for both the combined set-ups.
However, at present the CPV phase δ CP is not severely constrained, and as can be seen from Table 1, any value between 180 • and 342 • is allowed at 2σ C.L., and any value except for the ones between 31 • and 137 • is allowed at 3σ. Fixing the three mixing angles to their best fit values, we find from eq. (1.1) that the full range of cos δ CP ∈ [−1, 1] (allowed at present Thus, in principle, any value from this range may turn out to be favoured in the future. For instance, imagine that in the future the best fit value of δ CP will shift from 248 • to 290 • , while the best fit values of the mixing angles will remain the same. Then, the value of sin 2 θ ν 12 = 0.250, and thus the HG symmetry form, will be favoured. With this said, one should keep in mind that the position of the black dot in Fig. 1 is likely to change in the future, but having more precise measurements of δ CP and the mixing angles at our disposal, we will be able to repeat this analysis favouring some symmetry forms and disfavouring the others. Having obtained an idea of how much the mixing symmetry forms in question are compatible with the present best fit values of the oscillation parameters, we go next to a more involved analysis which will allow us to see the compatibility of the studied symmetry forms with any value of δ CP between 180 • and 360 • , should it turn out to be the true value. To this aim, we fix the true value of the CPV phase, δ true CP , to be between 180 • and 360 • , the true value of the atmospheric mixing angle, θ true 23 , to a value from its 3σ range, and the true values of the solar and reactor mixing angles, θ true 12 and θ true 13 , to their corresponding best fit values. Then, we generate data with this input using the DUNE and T2HK set-ups. In the test, we assume a given symmetry form to hold and fix the three test values θ test 12 , θ test 13 , and θ test 23 to values from the corresponding 3σ ranges. Using these test values and known for a given symmetry form θ ν 12 , we calculate δ test CP from eq. (1.1). Each couple of the true and test oscillation vectors, y true = (θ true 12 , θ true 13 , θ true 23 , δ true CP ) and y test = (θ test 12 , θ test 13 , θ test 23 , δ test CP ), provides a certain value of ∆χ 2 . Note that in calculating this value, we use external priors on sin 2 θ 12 and sin 2 θ 13 from JUNO and Daya Bay, as explained in Section 4. A detailed discussion on the impact of these two priors is presented in Appendix A. Further, for each y true we marginalise over θ test ij (over y test ). Finally, for each δ true CP we marginalise also over θ true 23 . We repeat this procedure for each of the four symmetry forms in study. The obtained results are shown in Fig. 2. Two comments on Fig. 2 are in order.
• For each symmetry form a significant part of the parameter space gets disfavoured at more than 3σ. Should the true value of δ CP lie in this part of the parameter space, the corresponding symmetry form will be disfavoured at 3σ confidence level.
• Now we can see at which C.L. any given symmetry form is compatible with any potentially true value of δ CP . We just need to draw a vertical line at δ true CP of interest. The points where it crosses the ∆χ 2 curves will provide the confidence levels at which the TBM, GRA, GRB, and HG forms are compatible with this δ true CP . In particular, for δ true CP = 248 • , we find numbers which correspond to those extracted from Fig. 1

How well can DUNE and T2HK Separate between Various Symmetry
Forms?
In this subsection, we will answer the question of how well DUNE and T2HK can distinguish the discussed symmetry forms under the assumption that one of them is realised in Nature. Given the fact that the BM form is not compatible with the current best fit values of the neutrino mixing angles, which we are going to use first in our analysis, we end up with four best fit values of interest. Namely, from Table 2, we read δ CP = 256 • , 261 • , 282 • , and 293 • for the GRB, TBM, GRA, and HG symmetry forms, respectively. Assuming one of them to be the true value of δ CP , we will test the remaining three values against the assumed true value using DUNE, T2HK, and their combination. Overall, we have 12 pairs of the values we want to compare. We start with DUNE. After performing a statistical analysis of simulated data, as described in Section 4, we obtain that for all the 12 cases ∆χ 2 does not exceed approximately 3.5. This value of ∆χ 2 is found when the value of δ CP predicted in the HG (GRB) case is tested against the value of δ CP for the GRB (HG) form, which is assumed to be the true one. Therefore, the sensitivity of DUNE alone is not enough to make a 3σ claim on discriminating between the symmetry mixing forms under investigation, and we will test next all the cases using simulated data from the T2HK experiment, whose overall sensitivity to CPV is better than that of DUNE.
Performing a statistical analysis for T2HK, we find that it can discriminate the GRB case from the HG case at approximately 2.5σ confidence level. More specifically, if δ CP = 256 • (293 • ) turned out to be the true value of the CPV phase, then T2HK could disfavour the value of δ CP = 293 • (256 • ) with ∆χ 2 ≈ 7.5. We also find that the TBM and HG symmetry forms, in turn, occur to be resolvable at slightly less C.L. with ∆χ 2 being around 5.5. Thus, the sensitivity of T2HK is not sufficient as well to discriminate between the cases of interest at 3σ C.L. For that reason, we will test them further using the potential of combining DUNE and T2HK.  Table 5: Confidence levels (in number of σ) at which the symmetry forms under consideration can be distinguished from each other assuming that one of them is realised in Nature. The result is obtained using the combination DUNE + T2HK. All the mixing angles have been fixed to their NO best fit values both in data and in test.
The combination of DUNE and T2HK provides better sensitivity to the CPV phase δ CP than either of these two experiments in isolation (see, e.g., [72]). A combined analysis performed by us leads to the results described below. Firstly, the GRB and HG mixing forms can be now distinguished at more than 3σ confidence level. If δ CP = 256 • (293 • ) is the true value, then δ CP = 293 • (256 • ) will be disfavoured with ∆χ 2 ≈ 11 (10.5). Secondly, the TBM and HG cases can be resolved at slightly less than 3σ, the corresponding values of ∆χ 2 being around 8. Thirdly, discriminating between the GRA and GRB forms can be claimed with ∆χ 2 ≈ 5.5. Finally, the sensitivity of the combination of these two experiments is not enough to discern TBM from GRA, GRA from HG, and TBM from GRB at even 2σ. For these three pairs, we find ∆χ 2 ≈ 3.5, 1.2, and 0.2, respectively, when the corresponding predictions for δ CP are compared between themselves.
We have checked that adding NOνA and T2K data to sets of simulated data obtained using the DUNE and T2HK set-ups leads to increase in the values of ∆χ 2 of only several tenths. Thus, inclusion of these data does not help to improve differentiating between the considered mixing schemes. We summarise the obtained results in Table 5, in which we present confidence levels (in number of σ) at which the symmetry forms under consideration can be distinguished from each other assuming that one of them is realised in Nature and using the potential of combining DUNE with T2HK.
Further, performing the more involved analysis described in Appendix A, we obtain the results summarised in Fig. 3. This figure allows us to see immediately at which C.L. a given pair of the symmetry forms can be distinguished, under the assumption that one form in the pair is realised in Nature. In particular, the numbers presented in Table 5 get clear graphic representation. Indeed, we see that using the combination DUNE + T2HK, GRB and HG can be resolved at more than 3σ C.L., while TBM and HG can be distinguished at almost 3σ.
As we see from Appendix A, the external prior on sin 2 θ 12 from JUNO is very important for the analyses performed in the present study. Usually, the present precision on sin 2 θ 12 is sufficient for the LBL experiments to achieve their goals on determination of δ CP , neutrino mass ordering, and the octant of θ 23 . However, in our case, the role of θ 12 is very important, since, as we have mentioned earlier, eq. (1.1), and thus predictions for δ CP provided by different symmetry forms, are very sensitive to the value of the solar angle. Thereby, there is a nice synergy between JUNO on the one hand and the LBL experiments on the other: DUNE and T2HK will be much more sensitive in addressing the questions posed in the present study, if they are provided with a precise measurement of θ 12 performed by JUNO. Finally, we would like to notice that the ∆χ 2 values obtained in the case of DUNE + T2HK in Fig. 3 can also be inferred from Fig. 2. Namely, drawing a vertical line at the minimum of ∆χ 2 curve for a given symmetry form in Fig. 2, we can assess how much the other forms are disfavoured with respect to the chosen form. For example, let us assume that the HG form is realised in Nature. Then, we have δ true CP = 293 • (see Table 2). Drawing a vertical line at this value of δ true CP , we read from the intersections with the GRA, TBM, and GRB curves:  Table 6: The values of cos δ CP and δ CP for different symmetry forms obtained from the sum rule in eq. (1.1) fixing sin 2 θ 12 = 0.3343 (its upper 2σ bound) and two other mixing angles to their NO best fit values. ∆χ 2 ≈ 1, 7, and 10, respectively. These are to be compared with the bottom right panel of Fig. 3.

The BM Symmetry Form
As we have mentioned in the Introduction, even though the BM symmetry form is not compatible with the current best fit values of the neutrino mixing angles, it turns out to be viable, if the current 2σ ranges of the mixing angles are taken into account. For example, if we keep sin 2 θ 13 and sin 2 θ 23 fixed to their best fit values for NO, we find that the value of sin 2 θ 12 = 0.3343, which is the upper bound of the corresponding 2σ range (see Table 1), is required to obtain cos δ CP = −1.00 and thus, recover viability of the BM mixing form. For this choice of values of the mixing angles, the values of cos δ CP (and δ CP ), predicted in the TBM, GRA, GRB, and HG cases, change. We summarise them in Table 6.
We perform the analysis in this case and find that the BM form can be distinguished from all the other forms at more than 5σ by DUNE alone. The corresponding ∆χ 2 are between 25 and 31, and they translate to the numbers of σ presented in Table 7. T2HK provides even better results, which we also show in Table 7.

Summary and Conclusions
In the present study, we have explored in detail the sensitivity of the future LBL experiments DUNE and T2HK to test various lepton mixing schemes predicted by flavour models with non-Abelian discrete symmetries. These models provide a natural explanation of the observed pattern of neutrino mixing. We have concentrated on a particular sum rule for cos δ CP given in eq. (1.1), which holds for a rather broad class of discrete flavour symmetry models. We have considered five different underlying symmetry forms of the neutrino mixing matrix, namely, bimaximal (BM), tri-bimaximal (TBM), golden ratio type A (GRA), golden ratio type B (GRB), and hexagonal (HG). Each of these mixing schemes is characterised by a specific value of the angle θ ν 12 entering into the sum rule in eq. (1.1). The values of θ ν 12 for the BM, TBM, GRA, GRB, and HG forms are 45 • , 35 • , 32 • , 36 • , and 30 • , respectively. The BM symmetry form is disfavoured by the present best fit values of the mixing angles. Table 2 summarises the predictions for δ CP for the other symmetry forms assuming the current best fit values of θ 12 , θ 23 , and θ 13 . In our analysis, we have considered only the predicted values of δ CP lying around 270 • , since they are preferred by the present oscillation data (see Table 1).  Table 7: Confidence levels (in number of σ) at which the symmetry forms under consideration can be distinguished from each other by different experiments in the case of possibility to have viable BM mixing in the neutrino sector. "D" and "T" stand for DUNE and T2HK, respectively. When not explicitly specified, the results are for DUNE + T2HK. Both in data and in test, sin 2 θ 12 has been set to 0.3343, while sin 2 θ 23 and sin 2 θ 13 have been fixed to their NO best fit values.

True
Based on the prospective DUNE + T2HK data, the GRB and TBM symmetry forms are compatible with the current best fit values of the mixing parameters at around 1σ confidence level. Under the same condition, the GRA and HG forms are disfavoured at around 3σ and 4σ, respectively (see Fig. 1). Next, in Fig. 2, we show up to what extent any given symmetry form is compatible with any true value of δ CP lying in the range 180 • to 360 • . In our analysis, we impose an external Gaußian prior of 0.7% (at 1σ) on sin 2 θ 12 as expected from the upcoming JUNO experiment, which improves our results significantly, as shown in Fig. 4 in Appendix A. This demonstrates a very important synergy between JUNO and LBL experiments like DUNE and T2HK, while testing various lepton mixing schemes in light of oscillation data.
The combined data from DUNE and T2HK can discriminate among GRB and HG at more than 3σ, if one of them is realised in Nature and the other form is tested against it (see Table 5). The same is true for TBM and HG at almost 3σ. Note, in these two cases, the differences between the predicted best fit values of δ CP are 37 • and 32 • , respectively (see Table 2). Similarly, the GRA symmetry form can be distinguished from GRB and TBM at around 2σ. The corresponding differences in these cases are 26 • and 21 • , respectively. At the same time, there is a difference of 11 • for GRA and HG, which can be discriminated only at 1σ. For TBM and GRB, the difference is only 5 • and therefore, the significance of separation is very marginal (around 0.5σ).
In conclusion, the detailed analyses performed in the present work can be applied to any flavour model leading to a sum rule which predicts δ CP . In this regard, our article can serve as a useful guidebook for further studies. In this appendix, we discuss the role of external priors on various mixing angles that we considered in our analysis. First of all, we do not consider any prior on sin 2 θ 23 since both DUNE and T2HK will be able to measure this parameter with sufficient precision. However, since these experiments are not sensitive to θ 12 (see probability expressions in [74]), we consider an external Gaußian prior of 0.7% (at 1σ) on sin 2 θ 12 as expected from the proposed JUNO experiment [75]. Even though both DUNE and T2HK can provide high precision measurement of θ 13 using their appearance channels, but to speed up our computation, we also apply a Gaußian prior of 3% (at 1σ) on sin 2 θ 13 as expected by the end of Daya Bay's run [76]. Fig. 4 shows the impact of these priors in our analysis. We obtain this figure in the following way. First, we assume one of the symmetry forms, e.g., TBM, to be realised in Nature. For this symmetry form, we estimate the true value of δ CP from eq. (1.1) using the current best fit values of the mixing angles for NO (see Table 1) as their true values. We generate data with this input. Then, in the test, we vary sin 2 θ ij and δ CP in their corresponding 3σ allowed ranges. For each set of test values of sin 2 θ ij and δ CP , we estimate ∆χ 2 and also calculate the corresponding value of sin 2 θ ν 12 using eq. (1.1). For the same sin 2 θ ν 12 , there can be several values of ∆χ 2 . From there, for each sin 2 θ ν 12 , we choose the minimal value of ∆χ 2 . Finally, we plot this minimum ∆χ 2 as a function of sin 2 θ ν 12 in Fig. 4. We present the results for T2HK considering the following cases: (i) without assuming any priors on sin 2 θ 12 and sin 2 θ 13 , (ii) assuming priors on both the parameters, (iii) only the prior on sin 2 θ 12 , and (iv) only the prior on sin 2 θ 13 . From this figure, we can make the following few important observations.
• First, we see that the curves corresponding to cases (i) and (iv) almost overlap with each other. It suggests that for the physics case under study, T2HK does not need an external prior on sin 2 θ 13 since it can provide a necessary precision on this parameter.
• Secondly, we observe that the curves corresponding to cases (ii) and (iii) also overlap with each other. It indicates that for our purpose, T2HK needs an external prior on sin 2 θ 12 from JUNO since it has a very mild sensitivity on this mixing parameter.
We have checked that the above observations are also valid for GRA, GRB, and HG symmetry forms and for DUNE as well. We have also seen that cases (ii) and (iii) are almost equivalent to the fixed parameter scenario, where we keep all the mixing angles to be fixed  : Impact of external Gaußian priors on sin 2 θ 12 and sin 2 θ 13 on the resulting ∆χ 2 function in the case of T2HK. We use a prior on sin 2 θ 12 of 0.7% (at 1σ) from JUNO and a prior on sin 2 θ 13 of 3% (at 1σ) from Daya Bay. While doing so, we assume the current best fit values of these two parameters to be their true values. The black dot corresponds to sin 2 θ ν 12 = 1/3, which characterises the TBM symmetry form.
at their best fit values in the fit. Unless mentioned otherwise, we always impose both these priors in our statistical analysis, as described in Section 4. B Impact of Marginalisation over ∆m 2 31 In this appendix, we give Fig. 5 to show the impact of the present 3σ uncertainty on ∆m 2 31 while testing the compatibility between the considered symmetry forms and present oscillation data. In Fig. 5, we show the potential of the combined DUNE + T2HK set-up for the two different cases: (i) fixed parameter scenario where we keep all the oscillation parameters fixed at their benchmark values in the fit, and (ii) we only marginalise over ∆m 2 31 in the fit in its present 3σ allowed range. The fixed parameter curve is exactly similar to what we have already presented in Fig. 1 for the combined DUNE + T2HK set-up. For the GRB and TBM schemes, we do not see much difference between the fixed parameter case and the case where we marginalise over ∆m 2 31 in the fit. For the GRA (HG) symmetry form, ∆χ 2 gets reduced by ∼ 11% (13%), when we marginalise over ∆m 2 31 instead of keeping it fixed in the fit. for the TBM, GRA, GRB, and HG symmetry forms. We show the results for the two different cases: (i) fixed parameter scenario where all the oscillation parameters are kept fixed in the fit, and (ii) we only marginalise over ∆m 2 31 in the fit in its present 3σ allowed range.
C Agreement between Various Mixing Schemes and Oscillation Data in (sin 2 θ true 23 , δ true CP ) Plane Here we will see which regions of the parameter space in the plane of true values of sin 2 θ 23 and δ CP will be compatible at less than 3σ C.L. for each symmetry form of interest, if that form is realised in Nature. To this aim, for each symmetry form (fixed θ ν 12 ), we calculate δ CP using eq. (1.1) with the test values of the mixing angles θ test ij . Then, we marginalise ∆χ 2 over θ test ij for given true values θ true 23 and δ true CP . The coloured bands in Fig. 6 represent potentially true values of δ CP as well as sin 2 θ 23 with which the form under consideration is compatible at 1σ, 2σ, 3σ confidence levels in the context of DUNE, T2HK and their combination. If the true value of δ CP turned out to lie outside these bands, this would imply that the given symmetry form is disfavoured at more than 3σ C.L. For all the symmetry forms, a significant part of the parameter space gets disfavoured at more than 3σ.
For each symmetry form, the black dashed line has been obtained using eq. (1.1) with θ 12 and θ 13 fixed to their NO best fit values. Note that a given symmetry form is well compatible with any point close to this line. The star denotes the present best fit values of sin 2 θ 23 and Figure 6: Compatibility of various symmetry forms with any potentially true values of sin 2 θ 23 and δ CP in the context of DUNE (left panels), T2HK (middle panels), and DUNE + T2HK (right panels). For a given symmetry form (fixed θ ν 12 ), the black dashed line has been obtained using eq. (1.1) and fixing θ 12 and θ 13 to their NO best fit values. The star indicates the present best fit values for NO as given in Table 1. For all the symmetry forms, a significant part of the parameter space gets disfavoured at more than 3σ C.L. for DUNE + T2HK. δ CP for the NO spectrum as given in Table 1, namely, (sin 2 θ 23 , δ CP ) = (0.425, 248 • ). As we can see from the right panels of Fig. 6, i.e., in the context of DUNE + T2HK, the HG (GRA) form is disfavoured at more than (precisely) 3σ C.L., while the GRB and TBM forms are compatible with the star at 1σ and 2σ, respectively. If the star moves in the future to a different point, we will immediately conclude which symmetry forms are (dis)favoured. Let us assume, e.g., that the future best fit values are (sin 2 θ 23 , δ CP ) = (0.58, 300 • ). Then, in the context of DUNE + T2HK, the GRB and TBM forms would be disfavoured at more than 3σ, while the HG (GRA) symmetry form would be compatible with this hypothetical position of the star at 1σ (2σ) confidence level.
Finally, we would like to notice the compatibility between this figure and the numbers in Table 5. Let us consider an example, in which GRB is the true form and HG is tested against it. In this case, from Table 5, we read the C.L. at which these two symmetry forms can be distinguished by DUNE + T2HK, namely, 3.3σ. We recall that the results in this table have been obtained assuming the current best fit values of the mixing angles to be their true values. Thus, for GRB we have (sin 2 θ 23 , δ CP ) = (0.425, 256 • ), which are the true values of these parameters in the case under consideration. Now, we put this point on the right bottom panel of Fig. 6, corresponding to the HG symmetry form for DUNE + T2HK, and find that this point falls just outside the dark green band representing δ CP values compatible with the HG form at 3σ C.L. It means that if δ CP predicted by GRB together with the present best fit values of θ 23 , θ 12 , and θ 13 are realised in Nature, then the HG symmetry form will be disfavoured by DUNE + T2HK at > 3σ. The same message is conveyed from Table 5 as well.