Non-standard neutrino interactions and mu-tau reflection symmetry

Non-standard interactions (NSI), possible sub-leading effects originating from new physics beyond the Standard Model, may affect the propagation of neutrinos and eventually contributes to measurements of neutrino oscillations. Beside this, $ \mu-\tau $ reflection symmetry, naturally predicted by non-Abelian discrete flavour symmetries, has been very successful in explaining the observed leptonic mixing patterns. In this work, we study the combined effect of both. We present a $S_4$ flavour model with $\mu-\tau$ reflection symmetry realised in both neutrino masses and NSI. Under this formalism, we perform a detailed study for the upcoming neutrino experiments DUNE and T2HK. Our simulation results shows that under the $\mu-\tau $ reflection symmetry, NSI parameters are further constrained and the mass ordering sensitivity is less affected by the presence of NSI.


INTRODUCTION
The observation of neutrino oscillations by various experiments using solar, atmospheric, reactor, and accelerator neutrinos indubitably established the existence of neutrino mass, which provides a clear evidence of new physics beyond the Standard Model (SM) [1]. The effective Hamiltonian to describe neutrino oscillation is in general written as Here, M ν is the neutrino mass matrix and V represents the effective potential for neutrino propagating in matter. Considering only the standard matter effect, V = diag{A, 0, 0} with A = 2 √ 2G F N e E, N e the electron number density in matter and E the neutrino energy. Neutrino oscillation data have been inspiring theoretical studies of lepton mixing and mechanisms behind. Traditionally, the observation of maximal atmospheric mixing motivates the notion of the µ − τ interchange symmetry, i.e., the permutation between the mu neutrino and the tau neutrino ν µ ↔ ν τ , which predict θ 23 = 45 • but θ 13 = 0 [2]. Later, with more precise oscillation data, a critical non-vanishing reactor angle and large CP violation were measured. The so-called the µ − τ reflection symmetry, originally proposed in Ref. [3] (for recent review see [4] and the references therein), have been caught more attentions. This symmetry, on the phenomenological side, satisfies |U µi | = |U τ i | (for i = 1, 2, 3), which predicts four cases [5]: (a) θ 23 = 45 • , θ 13 = 0; (b) θ 23 = 45 • , θ 12 = 0; (c) θ 23 = 45 • , θ 12 = 90 • ; (d) θ 23 = 45 • , δ = ±90 • . The last case, specifically, θ 23 = 45 • , δ = −90 • , is still in excellent agreement with the latest global analysis of neutrino oscillation data [6].
Various non-Abelian discrete flavour models have been proposed to realise the µ − τ reflection symmetry. Among them, the littlest µ − τ seesaw model is the most predictive one, where all mixing parameters and the ratio ∆m 2 21 /∆m 2 31 are dependent upon a single RG running parameter [7,8]. More models have been worked out in the framework of generalised CP symmetry [9,10]. Combining the µ − τ interchange with CP transformation, we arrive at the µ − τ reflection transformation, Requiring neutrino mass matrix to be invariant under this transformation is equivalent to imposing the following constraint (assuming Majorana neutrinos) where This immediately leads to where a, d are real and b, c are complex. It predicts not just θ 23 = 45 • , δ = ±90 • with non-zero θ 13 , but also the Majorana phases ρ, σ = 0, 90 • [11]. Symmetry arguments on how to realise the mass texture in Eq. (5) have been given in based on generalised CP combined with different non-Abelian discrete symmetries, e.g., A 4 [12], S 4 [10,13], ∆(48) [14], and ∆(96) [15]. Examples of explicit model construction can be found in, e.g., [12,13,16,17]. Besides this, new physics beyond the SM may lead to corrections to the effective neutrino interactions through higher-dimensional operators. These operators are often described in the framework of non-standard interactions (NSI). Here, we focus on NSI that can impact the propagation of neutrinos through matter. For neutrinos with kinetic energy around or below GeV scale, they are described by the dimension-6 four-fermion operators of the form [18], where f C αβ represent NSI parameters with α, β = e, µ, τ , C = L, R, f = e, u, d, and G F is the Fermi constant.
in the presence of NSI, the effective potential V in Eq. (1) is modified, i.e., where the 3 × 3 NSI matrix αβ with α, β = e, µ, τ is a Hermitian matrix, αβ = * βα . Nextgeneration long baseline (LBL) accelerator neutrino oscillation experiments, DUNE [19] and T2HK [20], that aim to do a precision test of standard three neutrino oscillation paradigm will also reach the sensitivity to probe NSI.
In recent times, there have been many interesting studies involving NSI and their impacts on the measurement of standard neutrino oscillation parameters, for detailed reviews see Refs. [21] and the references therein. In [22] authors have addressed various parameter degeneracies between standard and nonstandard interactions for the determination of neutrino mass ordering and the atmospheric mixing angle in case of DUNE and T2HK. In particular, Ref. [23] shows that when | eµ | = tan θ 23 | eτ |, a cancellation between leading order terms in the appearance channel probabilities will strongly affect the sensitivities to NSI parameters at T2HK. Also, generalized mass ordering degeneracy has been studied in great detail in [24]. Besides this, some exhaustive analysis in the presence of complex NSI parameters on the determination of the Dirac CP-violating phase have been performed in [25].
Notice that the above mentioned works are mostly based on model independent studies, whereas some model dependent analysis considering heavy charged singlet and/or doublet scalars have been performed in [26]. In these studies, all NSI parameters are assumed to be free complex parameters. Such a large redundancy leads to a difficulty of making definite predictions.
As pointed out in [27], the NSI effects may provide important information to extend our understanding of discrete flavour symmetries. A combined study of the flavour symmetry and NSI can be pursued in the future neutrino experiments. As a result, once the µ − τ reflection symmetry is a true symmetry in neutrino mass matrix, one may wonder if NSI also satisfy this symmetry. Furthermore, imposing the µ − τ reflection symmetry on the NSI matrix can efficiently reduce the parameter redundancy and satisfy the condition | eµ | = tan θ 23 | eτ | found in [23]. In this paper, we make an attempt to study the importance of NSI for the upcoming LBL neutrino oscillation experiments DUNE and T2HK in the presence of µ − τ reflection symmetry. More specifically, we consider NSI matrix that obeys µ − τ reflection symmetry and within this framework, the NSI effect has been examined in the context of LBL experiments.
We outline the rest of this work as follows. In Sec. 2, for illustration, we show how to realise µ − τ reflection symmetry in both the neutrino mass matrix and NSI matrix in a S 4 flavour model. We follow the technique developed in Ref. [27], which is the first article in the literature which combines NSI with flavour symmetries together. Furthermore, within the S 4 flavour symmetry and the CP symmetry, we show that it is possible to realise NSI matrix that follows µ − τ reflection symmetry. In Sec. 3, we perform a phenomenological study of the model both analytically and numerically. Starting from the analysis of the oscillation probabilities, we study the constraints on αβ , and then measure the impacts of NSI on the mass ordering sensitivity. We summarise our results in Sec. 4.

SYMMETRIES
For illustration, we will construct a flavour model in which the µ − τ reflection symmetry is preserved in both the neutrino mass matrix and NSI matrix in the flavour basis. In order to generate the µ − τ reflection symmetry in the neutrino mass terms, we introduce the flavour symmetry S 4 × Z 4 . We introduce necessary flavon fields which gain VEVs. As a consequence, the flavour symmetry is broken and a non-trivial flavour mixing is obtained. In the electroweak space, we imply the idea of scalar doublet-singlet mixing to generate sizeable NSI [26], where an additional global symmetry, Z 2 , is combined to avoid a myriad of experimental constraints from the charged lepton sector.
We give particle contents of the model in Table I. Three lepton gauge doublets L = (L 1 , L 2 , L 3 ) T are arranged as a triplet 3 of S 4 . Three right-handed charged leptons e R , µ R and τ R are singlets of S 4 . Four flavon fields are introduced: ϕ is a triplet 3 of S 4 to generate charged lepton Yukawa coupling; χ, ζ and ξ are triplet 3, doublet 2 and singlet 1, respectively, to generate the neutrino mass texture. The scalar doublet-singlet mixing is realised by introducing new electroweak doublets η and a charged gauge singlet scalar φ + , respectively. Here, we arrange η and φ + as a triplet and a singlet of S 4 , respectively. We discuss how to realise neutrino mass texture and the effective Hamiltonian terms of NSI with the µ − τ reflection symmetry, respectively in the following two subsections.
Fields L e R µ R τ R H η φ + ϕ χ ζ ξ SU (2) L 2 1 1 1 2 2 1 1 1 1 1 We first consider the construction of a simplified neutrino mass model. To be invariant in S 4 , the effective Lagrangian terms to generate charged lepton masses and neutrino masses should take the form (c.f. Eq. (A5) of the Appendix) where L = (L 1 , L 2 , L 3 ),H = iσ 2 H * , Λ is the scale of heavy flavour multiplets decouple and Λ ss is the see-saw scale where heavy mediators (e.g., right-handed neutrinos) decouple. This arrangement can be achieved by arranging additional Z n charges for flavons [13,16].
Combining the CP symmetry with the flavour symmetry together and keeping in mind that particles are flavour multiplets in the flavour space, all coefficients are forced to be real [9,10].
Charged lepton masses are obtained after the flavon ϕ and the Higgs gain VEVs. Given the VEV for the flavon in the charged lepton sector, a diagonal charged lepton mass matrix is obtained with diagonal entries given by where the Higgs VEV GeV has been used. In this case, the three lepton gauge doublets L 1 , L 2 and L 3 are identified with the three flavour states respectively. Neutrinos obtain their masses after the flavons χ, ζ, ξ and the Higgs H gain VEVs. Since we assume that the µ − τ reflection symmetry is a residual symmetry in the neutrino sector, VEVs of ξ, ζ and χ should be invariant under the relevant transformation, namely, where X µτ,2 and X µτ are representation-dependent transformation matrix respect to the µ − τ reflection, This condition leads to real ξ , χ 1 and complex-conjugate pairs ζ 1 We denote them as where v i and u i (for i = 1, 2, 3) are all real. We then arrive at the neutrino mass matrix in Eq. (5) with This is the most general neutrino mass matrix preserving the µ − τ reflection symmetry. It is straightforward to check its invariance under the µ − τ reflection transformation [11] ν e → ν c e , ν µ → ν c τ , ν τ → ν c µ .
It is equivalent to require This simplified model shows a generic way to achieve flavour mixing with the µ − τ reflection symmetry. For an explicit flavour model construction, one should carefully consider how to achieve the required flavon VEVs, and how to avoid the unnecessary higher dimensional operators. In order to solve these problems, supersymmetry and extra fields have to been included, extra Z n symmetries may be imposed, and UV completion of the effective theory may be considered. For relevant discussions, see e.g., [13,16,17]. As a consequence, correlations between oscillation parameters are induced in these models, which can be tested in the future neutrino oscillation experiments. In this paper, we will not consider these parameter correlations.

The mu-tau reflection symmetry realised in nonstandard interactions
We further discuss how to achieve NSI satisfying the µ−τ reflection symmetry. We couple flavons to the NSI mediator charged scalars η and φ + . Since VEVs of flavons contributing to neutrino masses satisfies the µ − τ reflection symmetry, any interactions between these flavons and charged scalars satisfy the µ − τ reflection symmetry.
In order to see the behaviour of the µ − τ reflection symmetry in NSI, we first consider, the following benchmark Lagrangian terms for illustration, where all coefficients are real by imposing CP symmetry. The first row gives rise to charged scalar masses after the SM Higgs obtain the VEV. In the basis S + = (η + 1 , η + 2 , η + 3 , φ + ) T , the charged scalar mass matrix is given by, Here, the Z 2 symmetry is not broken after flavour symmetry breaking and electroweak symmetry breaking. The second row of Eq. (17) gives rise to the interaction of charged scalars, neutrinos and the electron. It is explicitly written out as where the coefficient matrix P λ is given by The charged scalars have masses O(100) GeV. For GeV-scale neutrino beams propagating in the Earth, these scalars can be integrated and dimension-six operators are left, With the help of Fierz identity, we express it in the form where , and the Higgs VEV again has been replaced with (2 √ 2G F ) −1/2 . In this benchmark, the µ − τ reflection symmetry in NSI originates from the contribution of χ to the off-diagonal mass term between η and φ + . A correlation between NSI parameters, 2Arg( eµ ) + Arg( µτ ) = 0, since only one single source of the µ−τ reflection symmetry, i.e. the single coupling with χ, is included. This correlation is easily to be removed by including additional interactions between flavons (ξ, ζ and χ) and charged scalars. Assuming no interactions between the flavon ϕ and electric-charged scalars, which can be achieved by assuming Z n charges in the supersymmetry framework, all interactions related to NSI does not induce any source to break the µ − τ reflection symmetry, and eventually the µ − τ reflection symmetry is always preserved.
This model implies the advantage of using the scalar mixing between electroweak doublets and charged singlets to realise relatively sizeable NSI, αβ ∼ O(0.1). As originally proposed in [26], sizeable NSI require a price of fine-tuning between a large mixing and hierarchical masses for scalars, i.e., the mixing angle around 0.3, a low scalar mass around 100 GeV, and the heavier scalar mass around 10 TeV. Compared with the original model in [26], we extend the gauge doublet η from a flavour singlet to a flavour triplet of S 4 . With the price of fine tuning, NSI αβ ∼ O(0.1) can still be achieved.

TESTING MU-TAU REFLECTION SYMMETRY AT DUNE AND T2HK
Here we study the phenomenological consequences of the µ − τ reflection symmetry preserved in both the neutrino mass matrix and the NSI matrix. We start our phenomenological study from the analysis on oscillation probabilities with NSI under the µ − τ reflection symmetry. In the following, we briefly introduce the experimental setup for experiments, and the analysis details, e.g. the true values and priors for oscillation parameters and αβ . Before viewing the impact of this symmetry on the mass ordering sensitivity, we investigate the restriction of the µ − τ reflection symmetry on the allowed region of αβ .

Neutrino oscillation probabilities
In the framework of the µ − τ reflection symmetry, M ν takes Eq. (5) and αβ satisfy the correlations * Matrices M ν M † ν and V , appearing in the Hamiltonian, are expressed in forms of with a, c are real and b, d are complex. We have re-defined V with˜ ee ≡ ee − τ τ = ee − µµ to absorb an overall diagonal phase. Note that in this symmetry based formalism V contains five real free NSI parameters, one real (˜ ee ) and two complex ( eµ , µτ ). The appearance probability for long baseline neutrino oscillation experiments in the presence of NSI is complicated in general [23]. However, under the µ − τ reflection symmetry, we can write the appearance probability for the normal ordering (NO) in a simple form, i.e., where φ αβ ≡ Arg( αβ ) and following Ref. [28], x ≡ √ 2c 13 s 13 , y ≡ √ 2rs 12 c 12 , r = |δm 2 21 /δm 2 31 | , The antineutrino probability P µe , which is given by Eq. (27) with δ → −δ, φ αβ → −φ αβ , andÂ → −Â (and hence f →f ). For the inverted ordering (IO), ∆ → −∆, y → −y, A → −Â (and hence f ↔ −f , and g → −g). For small , we can neglect the y 2 , 2 and higher order terms in Eq. (27). For δ = ±90 • and the NO, we have (P µe ) NO x 2 f 2 ∓ 2xf g(y sin ∆ + 2Â| eµ | sin φ eµ cos ∆) , For δ = ±90 • and the IO, we have From Eqs. (29) and (30), we see that if ee = 0 and φ eµ = 0 or 180 • , the differences between the NSI and SM appearance probabilities disappear in both the neutrino and antineutrino modes, which is consistent with the degeneracy between | eµ | and | eτ | found in [23]. Note that for large | eµ |, higher order terms cannot be neglected, and the degeneracy is broken. From Eq. (27), we see that µτ does not appear in the appearance probability up to second order in . Hence, they mainly affect the disappearance channel. Taking ee , eµ and eτ equal to zero, the disappearance probability can be written as The above equation shows that if φ µτ = 90 or 270 • , the disappearance probability only depends on the second order in .

Experimental and simulation details
In this study, we consider two proposed next-generation superbeam experiments DUNE [19,29] and T2HK [20]. To perform the numerical simulation of both the experiments, we use the GLoBES packages [30,31] along with the required auxiliary files presented in Ref. [29] for DUNE and the detector simulation files in Ref. [23] for T2HK. DUNE will utilise existing Neutrinos at the Main Injector (NuMI) beamline design at Fermilab as a neutrino source. The far detector of DUNE will be placed at Sanford Underground Research Facility (SURF)  in Lead, South Dakota, at a distance of 1300 km (800 mile) and about 1.5 km under the surface from the neutrino source. DUNE collaboration has planned to use four 10-kton liquid argon time-projection chamber (LArTPC) detectors. According to their report, the expected design flux corresponds to 1.07 MW beam power which gives 1.47 × 10 21 protons on target (POT) per year for an 80 GeV proton beam energy. In our simulation, we simply consider 40 kton of fiducial mass for the far detector. In addition, the signal and background normalisation uncertainties for appearance and disappearance channels have been adopted from DUNE CDR [29]. We perform our numerical study considering 3.5 years running time in both neutrino and antineutrino modes.
The T2HK experiment [20] will utilise an upgraded J-PARC beam with a power of 1.3 MW, which gives 2.7 × 10 21 POT per year for an 30 GeV proton beam energy. They are planning to build a water Cherenkov (WC) detector, which is located 295 km away from the source. Making the detector 2.5 • off-axis allows it to measure a narrow band beam with a peak energy around 0.6 GeV. Here we consider the 2TankHK-staged design proposed in the HK design report [20]. We assume one tank will take data alone for the first six years, whereas a second tank will be added for another four years. Each tank will have 0.19 Mton fiducial mass with 40% photocoverage. While performing the numerical analysis, we take 2.5 years running time in neutrino modes and 7.5 years running time in anti-neutrino modes, i.e., 1 : 3 ratio out of the total runtime.
To incorporate NSI, we use the GLoBES extension file snu.c as has been presented in Refs. [32,33]. To implement the µ − τ reflection symmetry, we directly impose the form Eq. (26) on the matter potential in snu.c. We set the true parameter values of NSI as zero, whereas the test parameter ranges of NSI have been considered from [21].
Except for the value of δ and θ 23 , we use NuFit4.0 best-fit results for the true values [6], as shown in Tab. II. We fix the value θ 13 at the true value, as we have checked it does not make a great impact on the measurement on NSI parameters. We also fix θ 12 and ∆m 2 21 at the true values as these two parameters do not enter to the neutrino oscillations of DUNE and T2HK. We vary ∆m 2 31 within a Gaussian prior with the true value and the width which are the best fit and the 1σ uncertainty of NuFit4.0 results. In addition, we consider two orderings: for normal (inverse) mass ordering, the central value and the 1σ relative width of the ∆m 2 31 prior are 2.525 eV 2 (−2.586 eV 2 ) and 1.3% (1.3%), respectively. We do not include any other priors. Furthermore, we vary the phases φ eµ and φ µτ in the range [0, 2π). Note that the µ − τ reflection symmetry supports the phases φ eµ , φ µτ = 0, π. Finally, we do not consider the LMA-dark solution that is allowed by the generalized mass ordering degeneracy [24], since it is difficult to achieve αβ ∼ O(1) in the flavor symmetrical models.

Constraints on NSI parameters
We first study the impact of imposing the µ − τ reflection symmetry on the constraints on NSI parameters at next-generation superbeam experiments. In Fig. 1, we present the allowed parameter region on the | eµ | − | eτ | plane without restricting the µ − τ reflection symmetry on αβ for DUNE (solid) and the synergy of DUNE and T2HK (dashed). Gray, orange and black curves represent 1σ, 2σ, and 3σ contours, respectively. The tendency of contours leaves away from | eµ | = | eτ |, because of the contribution from the higher-order term in P (ν µ → ν e ) and P (ν µ →ν e ). Along with | eµ | = | eτ |, the bounds at 1σ, 2σ and 3σ are ∼ 0.125, ∼ 0.16, and ∼ 0.19, respectively. Obviously, this shows that when we impose the restriction of the µ − τ reflection symmetry on NSI, the uncertainty of | eµ | and | eτ | is much reduced. We see a deficit for all contours. For example, for the 3σ contour for DUNE, the discontinuity appears around | eµ | ∼ 0.25 and | eτ | ∼ 0.5. This is because the flipping of the sign of ∆m 2 31 . In Fig. 2, we discuss on the scenario in the presence of the µ − τ reflection symmetry on NSI. We show 1σ (gray), 2σ (red), and 3σ (black) contours on˜ ee − | eµ | (left), | µτ | − | eµ | (right) planes for DUNE (solid) and the combination of DUNE and T2HK (dashed). We see a correlation between | eµ | and | µτ |. This is because of the flipping of sign(∆m 2 31 ). In the other panel, we see that there is a dip on top of solid curves. This feature can be removed by We vary phases φ eµ and φ µτ in (0, 2π). ∆m 2 31 is allowed to be varied by within a Gaussian prior with a 1σ relative width of ∼ 1.25%. Two mass orderings are considered. No other priors are used.  The comparsion of the constraints on NSI parameters between cases with (w/) and without (w/o) the restriction of the µ − τ reflection symmetry for DUNE.
including T2HK data. Taking about the size of constraints, though T2HK is not sensitive for measuring NSI parameters, including its data significantly improves the | µτ | measurement, for which T2HK data can reduce the size of the uncertainty by ∼ 20%. The improvement on the other parameters is weaker. We note that the size of uncertainties is reduced once the µ − τ reflection symmetry is preserved in NSI. We compare the 1σ constraints between scenarios with and without the µ − τ reflection symmetry for DUNE in Tab. III. We see the improvement of the precision of˜ ee and | µτ | measurements, which is not only because of the fact that the number of degrees of freedom is reduced by imposing the µ − τ reflection symmetry in the neutrino-mass matrix and NSI, but also the size of the allowed region for | eτ | is largely shrunk down.

Impact on neutrino mass ordering sensitivity
The main goals of next-generation neutrino oscillation experiments are to measure the Dirac CP phase δ and to determine the neutrino mass ordering and the octant of θ 23 [19]. Since the µ − τ reflection symmetry is preserved in both the neutrino mass matrix and NSI matrix in our formalism, θ 23 = 45 • and δ = ±90 • will be measured by future neutrino oscillation experiments even in the presence of NSI. We find that the wrong sign of δ can be excluded at more than 5σ CL at DUNE and T2HK. Hence, here we fix δ = −90 • which is favored by the latest global analysis of neutrino oscillation data [6], and only focus on the study of the impact of NSI on the determination of mass ordering sensitivity at DUNE and T2HK within this formalism.
We first examine the potential of both experiments to measure the mass ordering in the presence of standard matter interactions as shown in Fig. 3. The combined potential of both experiments are also investigated. The solid green curve shows the ordering χ 2 for DUNE, the dashed-dotted brown curve represents the same for T2HK. The combined potential of DUNE and T2HK are shown by the dotted green curve. We also mark the 3σ, 5σ significance levels by the solid and dotted black horizontal lines, respectively. Our result is consistent with that in Ref. [34]. We notice that DUNE and T2HK can achieve a maximum ordering sensitivity of around χ 2 ∼ 20 and 5 significance levels, respectively, around the true value of the Dirac CP-violating phase δ = −90 • , respectively. Though the mass ordering sensitivity for T2HK is not good as DUNE, we see an improvement by combining two sets of data. Thus, in what follows we discuss about the DUNE and the combination of DUNE+T2HK.
To investigate the impact of NSI on the determination of neutrino mass ordering under the adopted formalism, we present our results in the left panel of Fig. 4. For simplicity, we first discuss our results for non-zero eµ , eτ , whereas the remaining NSI (i.e.,˜ ee , µτ ) have been fixed at zero. We show our results for the DUNE (DUNE+T2HK) using solid curves (dotted curves). The orange curves show our results in the presence of the symmetry, whereas the blue curves represent the same but in absence of the symmetry. Note that for eµ = * eτ curves, the parameter eτ has been marginalised in our analysis. In the case of DUNE, we notice that the mass ordering sensitivity does not get much affected for eµ 0.15 when the NSI parameters respect the symmetry over the whole parameter space (see solid orange curve). Including T2HK data slightly improves the sensitivity. On the other hand, noticeably sizeable impact has been observed when the NSI parameters fail to respect the symmetry (i.e. for eµ = * eτ ), which we show by the solid blue curve for the DUNE (see figure levels for details). The sensitivity curve gets severely affected for eµ 0.15 as can be seen from the solid blue curve. The impact of NSI on the mass ordering sensitivity is not reduced by including T2HK data. In the right panel of Fig. 4, we also demonstrate the  mass ordering χ 2 similar to the left, but for the most general picture, i.e., in the presence of all the NSI parameters (˜ ee , eµ , µτ ) within the symmetry framework. We see the same behaviour of the mass ordering sensitivity in the right panel of Fig. 4, i.e., by comparing the sensitivity curves in the left of with the right panel. However, we notice the stronger impact of varying˜ ee and µτ on the mass ordering sensitivity with the µ − τ reflection symmetry.

SUMMARY
The µ − τ reflection symmetry which predicts the maximal value of the atmospheric mixing angle as well as the Dirac CP-violating phase are in well agreement with the latest neutrino oscillation data. In this work, we explore the possibility of embedding the µ − τ reflection symmetry in the NSI matrix and study its importance for the determination of neutrino mass ordering in case of DUNE and T2HK.
We present an S 4 flavour model, with the µ − τ reflection symmetry realised in both neutrino mass matrix and NSI matrix. In this model, The sizeable NSI effects originate from the mediator of scalar doublet-singlet mixing. The essential to achieve the µ − τ reflection symmetry is the introduction of flavon fields and assuming the vacuum satisfies the symmetry. By coupling the flavons to neutrinos, the µ − τ reflection symmetry in neutrino mass matrix is realised. By coupling the flavons to the scalar mediators, that in the NSI matrix is achieved. The model provides a very generic way to construct NSI with the µ − τ reflection symmetry.
We find that imposing the µ − τ reflection symmetry improves the constraints on NSI parameters. In particular, the constraint to˜ ee is improved most, its allowed regime is improved from [−2.5, 1.2] (without the µ − τ reflection symmetry) to [−0.13, − 0.13] (within the µ − τ reflection symmetry). It is not just because of the reduced numbers of NSI parameters, but also due to the fact that the allowed region for eτ is reduced. This exactly points out that NSI under this scenario are much more testable with the upcoming neutrino oscillation experiments.
We also study the impact of NSI on the mass ordering sensitivity. Assuming˜ ee = µτ = 0, we find that in the presence of the µ − τ reflection symmetry, DUNE alone can reach almost 14σ mass ordering sensitivity. However, the mass ordering sensitivity gets compromised in absence of the concerned symmetry as presented in Fig. 4. We find slightly better sensitivity for DUNE+T2HK. Considering the most general scenario under the µ − τ reflection symmetry, we find that DUNE can achieve more than 10σ sensitivity, whereas the combined analysis of DUNE+T2HK data can give relatively better sensitivity than DUNE.
In summary, the existence of the µ − τ reflection symmetry in the neutrino mass matrix and matter effect NSI is attractive in both theoretical and phenomenological points of view. It is naturally realised in the framework of flavour symmetries, and is consistent with the current data. Furthermore, once the µ − τ reflection symmetry is the true symmetry behind, the testability of NSI is enhanced, and the impact of NSI on the mass ordering sensitivity is reduced.
Lagrangian terms contributing to lepton masses and invariant under S 4 × Z 4 are given by Applying the CG coefficients to the above formula, we arrive at Eq. (A5) with y e = −y e1 + 2y e2 . With the help of the CG coefficient, Lagrangian terms contribute to NSI in Eq. (17) can be written in a more compact form, whereη = (η † 1 , η † 3 , η † 2 ) is understood. They invariant under the symmetry S 4 × Z 4 × Z 2 .