$ \mu-\tau $ Reflection Symmetry Embedded in Minimal Seesaw

We embed $\mu-\tau$ reflection symmetry into the minimal seesaw formalism, where two right-handed neutrinos are added to the Standard Model of particle physics. Assuming that both the left- and right-handed neutrino fields transform under $\mu-\tau$ reflection symmetry, we obtain the required forms of the neutrino Dirac mass matrix and the Majorana mass matrix for the right-handed neutrinos. To investigate the neutrino phenomenology at low energies, we first consider the breaking of $\mu-\tau$ reflection symmetry due to the renormalization group running, and then systematically study various breaking schemes by introducing explicit breaking terms at high energies.


I. INTRODUCTION
Over past two decades, phenomenal neutrino oscillation experiments have established the formalism of three flavor neutrino oscillations and determined two mass squared differences and three mixing angles. At present, unknowns in neutrino oscillation physics are: the neutrino mass hierarchy, i.e., whether neutrinos obey normal hierarchy (NH, m 1 < m 2 < m 3 with m i 's being neutrino masses) or inverted hierarchy (IH, m 3 < m 1 ∼ m 2 ); the octant of atmospheric mixing angle θ 23 , and the determination of Dirac CP-violating phase δ. 1 Were neutrinos the Majorana particles, there would then exist two additional Majorana phases, which do not affect neutrino oscillation probabilities but can be probed by the neutrinoless double beta-decay (0νββ) experiments [4].
In the Standard Model (SM) of particle physics neutrinos are massless. One economical way to incorporate non-zero neutrino masses is to add two right-handed neutrinos to the SM and allow lepton number violation, 2 resulting in the so-called minimal seesaw scenario (see Ref. [5] for a review) within the context of the Type-I seesaw mechanism [6][7][8][9][10]. Integrating out the heavy right-handed neutrino fields results in the light neutrino mass matrix M ν as In this minimal seesaw set-up, M D is the (3 × 2) neutrino Dirac mass matrix whereas M R is the (2 × 2) Majorana mass matrix for the right-handed neutrinos.
In the basis where the charged lepton Yukawa matrix Y l is diagonal, diagonalizing the light neutrino mass matrix M ν then leads to the lepton mixing matrix, which is found to be sharply different from the quark mixing matrix. Namely, the former is highly non-diagonal while the latter almost diagonal. To explain the peculiar patterns in the lepton mixing matrix, various flavor symmetry models have been considered, e.g., in Refs. [11][12][13][14][15].
Other than assigning the µ − τ reflection symmetry only to the left-handed neutrino fields, in this work we apply the same symmetry to the right-handed neutrino fields as well.
Consequently, both the neutrino Dirac mass matrix M D and the Majorana mass matrix M R need to satisfy certain relations among their entries. While the resultant light neutrino mass matrix M ν still obeys the relations given in Eq. (2), the µ − τ reflection symmetry is now embedded in the minimal seesaw formalism, and both the left-and right-handed neutrinos are treated on the same footing. Similar ideas have been studied for the µ − τ permutation symmetry [43,44], while for the µ − τ reflection symmetry a detailed study on the scenario as ours is still missing. 3 Some recent studies on the minimal seesaw model can be found in Refs. [48][49][50][51][52][53][54][55][56][57][58][59][60][61][62].
In this paper we investigate the implications of the above embedding on neutrino phenomenology at low energies, while the possible discussions on the ultraviolet (UV) aspects, such as explaining the baryon asymmetry via leptogenesis [63] deferred to the future work.
For the low energy neutrino phenomenology, since the resultant light neutrino mass matrix M ν in Eq. (10) still preserves the usual µ − τ reflection symmetry, one may conclude that there is no new prediction in this seesaw embedded setup. However, we want to point out here at least two deserving issues which need careful scrutiny. The first one is to study the breaking of µ − τ reflection symmetry due to the renormalization group (RG) running. As now we impose the µ − τ reflection symmetry above the seesaw mass thresholds, an investigation on the RG-running effects is then inevitable in order to confront with the current global-fit of neutrino oscillation data [1][2][3] at low energies. Secondly, although the latest experimental data of T2K [64] and NOνA [65] are in good agreement with the predictions of the exact µ − τ reflection symmetry, 4 there still exist large uncertainties in the measurement of θ 23 . For instance, the best-fit values of θ 23 for the lower and higher octants are 43.6 • and 48.3 • in Ref. [65], respectively. Thus, it is tenacious to believe the exactness of µ−τ reflection symmetry, especially that there may exist large discrepancies when upcoming experimental data will be included. Therefore, it is also worthwhile to study how we can perturb such an exact µ − τ reflection symmetry, so that remarkable deviations can be observed in the lepton mixing parameters.
We organize our paper as follows. In Section II, we introduce the µ−τ reflection symmetry transformations to the left-and right-handed neutrino fields, and discuss the required forms of M D and M R . In Section III, we proceed to discuss the breaking of µ−τ reflection symmetry due to the RG running, followed by the systematic investigation on all the possible explicit breaking patterns of M D and M R in Section IV. Finally, we summarize our findings in Section V. Details of derivations, explanations and numerical results are relegated to the Appendices.

II. µ − τ REFLECTION SYMMETRY EMBEDDED IN MINIMAL SEESAW
In the minimal seesaw formalism, two right-handed neutrinos, collectively denoted as N R = (N µR , N τ R ) T in the flavor basis, is added to the SM. The relevant Lagrangian containing the neutrino Yukawa matrix and the Majorana mass term for the right-handed neutrinos are given by where L is the lepton doublet in the SM, Y ν stands for the neutrino Yukawa matrix, and H = iσ 2 H * with H denoting the SM Higgs field. After the Higgs field acquiring its vacuum expectation value, i.e., v = H ≈ 174 GeV, we obtain the neutrino Dirac mass term as ν L M D N R + h.c., where M D = vY ν is the neutrino Dirac mass matrix, and ν L = (ν eL , ν µL , ν τ L ) T stands for the left-handed neutrino fields in the flavor basis.
To embed the µ − τ reflection symmetry into the minimal seesaw formalism, we first 4 We note that the rejection of maximal mixing in θ 23 has moved from 2.6σ [66] to 0.8σ [65].
propose the following transformations for the left-and right-handed neutrino fields, where ν c L = Cν L T and N c R = CN R T are the charge-conjugated fields of ν L and N R , respectively, and the transformation matrices S and S are given by Applying the above transformations to the mass terms of neutrinos yields Then, if the neutrino mass terms M D and M R obey the following relations, we state that M D and M R are invariant under the transformations given in Eq. (4), or, µ − τ reflection symmetry is embedded in both M D and M R .
Without loss of generality, the µ − τ reflection symmetric limit of M D and M R can be parameterized as, where the phases and m 23 are all real. According to the seesaw mass formula, we then obtain the mass matrix M ν for the light neutrinos as, Here m 23 = m 23 /(|m 22 | 2 −m 2 23 ) and m 22 = |m 22 |/(|m 22 | 2 −m 2 23 ). We note that the parameters A and D in M ν are real, and M ν preserves the usual µ − τ reflection symmetry in Eq. (2).
The light neutrino mass matrix M ν can be diagonalized as M ν = V m d ν V T , where m d ν = diag{m 1 , m 2 , m 3 } is the diagonalized neutrino mass matrix. In the standard PDG [67] parameterization, the unitary matrix V can be decomposed as where c ij (s ij ) (for j = 12, 23, 13) stands for cos θ ij (sin θ ij ), P l = diag{e iφe , e iφµ , e iφτ } contains three unphysical phases which can be absorbed by the rephasing of charged lepton fields, and finally P ν = diag{e iρ , e iσ , 1} is the Majorana phase matrix.
A detailed derivation of these predictions is given in Appendix A. Note that in the minimal seesaw framework the lightest neutrino is always massless, and consequently one can always remove one of the Majorana phases. Here we take ρ to be absent.
We notice that in both NH and IH there exists a relation among neutrino masses, i.e., Having introduced the µ − τ reflection symmetry in the minimal seesaw formalism, in the subsequent sections we study the breaking of such a symmetry and its impact on neutrino oscillation parameters at low energies.

III. BREAKING DUE TO RENORMALIZATION GROUP RUNNING
We start with the investigation on the breaking of µ − τ reflection symmetry due to the The RG running towards low energies can then be divided into three stages. The first stage of running starts from the GUT scale and ends at the mass threshold of the heavier right-handed neutrino N 2 , schematically, Λ GUT → M 2 . The one-loop RG equations of relevant Yukawa and mass matrices are given by [68][69][70] dY where t = ln(µ/Λ GUT )/(16π 2 ) with µ being the renormalization scale, and the flavorindependent parameters α l and α ν are defined as In the above Y u , Y d and Y l are the up-quark, down-quark and charged-lepton Yukawa matrices, respectively. We note that unlike the RG running below the seesaw threshold, the charged-lepton Yukawa matrix Y l now could be non-diagonal due to the term of Y ν Y † ν , even if it were diagonal at the high energy boundary. As a result, when extracting the lepton mixing parameters above the seesaw threshold, one also needs to take into account the corrections from Y l . The light neutrino mass matrix at this stage of running is given by

and the RG running of M
(2) ν is found to be [68][69][70] dM (2) Similar to Y l , the evolution of M (2) ν now also involves the contribution from the neutrino Yukawa matrix Y ν .
When the renormalization scale is below the mass threshold of N 2 , the second stage of running starts and it ends at the mass threshold of N 1 , i.e., M 2 → M 1 . At the matching scale µ = M 2 , the light neutrino mass matrix is given by where Y ν and M R are the Y ν and M R with the entries corresponding to N 2 removed, while Y ν is the column corresponding to N 2 in Y ν . Specifically, Y ν (µ = M 2 ) and Y ν (µ = M 2 ) are the first and second columns of Y ν (µ = M 2 ), respectively, and M R (µ = M 2 ) is the (11) entry of M R (µ = M 2 ). Below µ = M 2 , the one-loop running of Y l , Y ν , M R and M ν is still formally dictated by Eqs. (18,19,20,23), except that we replace Y ν and M R by Y ν and M R .
The final stage of RG running starts from the mass threshold of N 1 and stops at a chosen low energy scale. Here we take the low energy scale to be the electroweak scale Λ EW .
This stage of RG running is below the seesaw threshold, and its impact on the lepton mixing parameters has been extensively discussed in the literature, e.g., Refs. [68][69][70][71]. In particular, the breaking of µ − τ reflection symmetry due to this stage of RG running is investigated in Refs. [28,72]. To save space, we then would not elaborate more on this stage of running.
With the above RG equations, in principle one can investigate the breaking of µ − τ reflection symmetry due to the RG running analytically, as was done in Ref. [73] for the littlest seesaw scenario. However, in the current setup all entries of Y ν are non-zero, so that an analytical study turns out to be formidable. We then choose to study this issue numerically.
In the numerical study, we set tan β = 30, and the high and low energy boundary scales are taken to be Λ GUT = 2 × 10 16 GeV and Λ EW = 1 TeV, respectively. We also study the cases considering tan β < 30, and as the modifications on mixing parameters are quite small we do not include those results here. The gauge couplings and various Yukawa couplings at Λ GUT are taken to be the default values in the numerical RG running package REAP [69], although M D (or Y ν ) and M R are set to be the forms given in Eq. (8) and Eq. (9), respectively. Our numerical strategy is to scan all the parameters in Y ν and M R at high energies, and seek the allowed parameter space that yields lepton mixing parameters compatible with current experimental data at low energies. The varied ranges of parameters at the high energy boundary are as follows, To guide the parameter scan, we employ the nested sampling package Multinest [74][75][76], with a χ 2 function built based on the latest global fit results [3]. Details about the χ 2 function can be found in Appendix B.
In Fig. 1 we show the numerical result for the NH case. The black scatter points have χ 2 < 30, and the best-fit (BF) point that has the minimal value of χ 2 , denoted as χ 2 min , is shown in red. In this NH case, we obtain χ 2 min = 19.34 for the BF point. In the M 2 vs. M 1 plot of Fig. 1, we display the spread of two mass thresholds M 1,2 for the right-handed neutrinos. It can be seen that M 1 and M 2 are quite to close each other and both of the order of 10 14 GeV. One possible explanation for such closeness of M 1 and M 2 is that the entries in two columns of Y ν are related by the µ − τ reflection symmetry, particularly due to   the symmetry transformation on N R . Therefore, no large hierarchy exists between the two columns of Y ν , and then in order to yield mild hierarchy in the light neutrino mass matrix, the entries in M R also tend to be close to each other, resulting in similar values of M 1 and M 2 . Because of the closeness of M 1 and M 2 , the second stage of RG running between two mass thresholds turns out to be insignificant, and thus we focus on the first and third stages of running in the following.
For the convenience of quantifying the RG running effects, we introduce the quantities of ∆x LH , ∆x LT and ∆x TH in the rest plots of Fig. 1. Here x stands for the lepton mixing angles and phases, and we define ∆x LH as the difference of the mixing parameter x at the low and high energy scales, i.e., ∆x To compare the RG running effects between the first and third stages of running, in the y-axis of these plots we show the absolute values of the ratios of ∆x LT /∆x TH . By inspecting Fig. 1, we then observe: • In this NH case all ∆x LH 's (shown as the x-axis) are rather small, indicating the mixing angles and phases receive small deviations from the RG running. The small deviations at the third stage of RG running are expected, as it is known that in NH the RG running of mixing angles and phases are insignificant below the seesaw threshold [68][69][70]. For the first stage of running, the corrections to M (2) ν are at the order of ln(Λ GUT /M 2 )/(16π 2 ) ∼ 0.03, assuming Y ν to be of O(1). Therefore, in the NH case the contributions from the first stage of running are also small.
• Regarding the relative contributions between the first and third stages of running, we notice that for the Dirac phase δ, the third stage of running tends to yield larger deviations than the first stage, as |∆x LT /∆x TH | > 1 for most of scatter points. However, for the three mixing angles and the Majorana phase σ, the first stage of running turns out to be more important than, or as important as, the first stage. It then demonstrates that in this seesaw embedded setup, the RG running effects above the seesaw threshold can be comparable with that below the seesaw threshold.
• As for the breaking of µ − τ reflection symmetry, θ 23 at low energies tends to be always larger than 45 • , while δ and σ can receive either positive or negative deviations from RG running. The correlation between the positive deviation of θ 23 and NH is in agreement with the previous RG running studies below the seesaw threshold [77,78].
To have better feeling of the RG running in this NH case, in the left three plots of Fig. 2 we show the detailed RG running of the mixing angles, phases and neutrino masses for the BF point. One can easily see that the two mass thresholds, indicated by the dashed vertical lines, are quite close to each other, and the running of mixing angle and phases are indeed not appreciable. However, significant running is observed for the neutrino masses, and because there exist contributions from Y ν in α ν during the first stage of running, RG running at the first stage is more dramatic than the third stage.
We next turn to the IH case, and the corresponding numerical results are shown in Fig. 3 point that has the least value of χ 2 (χ 2 min = 51.3) in that branch of scatter points. One then identifies a sharp decline in θ 12 during the first stage of running. Such a decline can be traced to the potential crossing of neutrino masses when µ ∼ 5 × 10 15 GeV (red arrow).
As interchanging the order of eigenvalues would lead to a 90 • rotation in the mixing angle, it also explains why the sum of the values of θ 12 before and after the decline is around 90 • .
From the running plot of δ, such interchange of eigenvalues seems to induce a 180 • change in δ as well.
As for the Majorana phase σ, from Fig. 3 we notice that the running of σ is also quite mild in this IH case. Moreover, we find that the obtained values of σ's for the scatter points shown  here are all around 90 • . Although σ = 0 is also predicted by the µ − τ reflection symmetry, having σ = 0 at the high energy boundary would lead to dramatic running in θ 12 . This can be seen from the following RG equation for θ 12 below the seesaw threshold [68][69][70][71]79], If σ ∼ 0, there then exists an enhancement factor of (m 1 + m 2 )/∆m 2 21 . Such dramatic running of θ 12 would hinder the sampling program to map out the allowed parameter space corresponding to σ ∼ 0 at high energies. In contrast, if σ ∼ 90 • and because m 1 ∼ m 2 in IH, the combination of |m 1 + m 2 e 2iσ | becomes vanishingly small. This also explains why even in IH the running of θ 12 is still insignificant.
Lastly, we point out that in this IH case θ 23 tends to be smaller than 45 • at low energy.
As in the NH case, this observed correlation between the negative deviation of θ 23 and IH is in agreement with the previous RG running studies below the seesaw threshold [77,78].
In Appendix D we present the numerical values of the neutrino Yukawa matrices and the Majorana mass matrices for the right-handed neutrinos at Λ GUT for the three scenarios shown in Fig. 2. The lepton mixing parameters at various energy scales are also shown.
From the previous RG running study we notice that in both NH and IH the breaking effects due to the RG running are quite mild. For instance, the deviations in θ 23 are only around one degree. Although such small deviations are in compatible with current experimental data, it may become necessary to consider large deviations when more accurate data will be included. In this section, we set out to discuss the breaking of µ − τ reflection symmetry in the low energy neutrino mass matrix by introducing explicit breaking terms in the neutrino Dirac mass matrix M D and the Majorana mass matrix M R for the right-handed neutrinos. As the RG running effects are found to be mild, for simplicity we choose to ignore them in the following discussion.
A. Breaking µ − τ reflection symmetry in M D We start with assigning an explicit breaking term in the (12) position of M D , so that the neutrino Dirac mass matrix M D after breaking is given by where is a small breaking parameter, taken to be real for simplicity. We name this breaking scenario as S1. The above M D leads to a new mass matrix M ν for the light neutrinos, and the difference between M ν and M ν is given by where B 12 = b * /det(M R ), and A i 's are defined as To evaluate the impact of the above breaking on the neutrino masses and lepton mixing angles, we diagonalize M ν with the mixing matrix V , which coincides with the mixing matrix V when = 0. For simplicity, we consider the scenario in which all entries of ∆M ν are real, and expand the deviations of neutrino masses and lepton mixing angles in terms of small parameters , θ 13 and ζ = m 2 /m 3 (ξ = ∆m 2 21 /m 2 2 ) for NH (IH). In the top block of Table I we show the leading order results for the deviations of neutrino masses ∆m i = m i − m i (for i = 1, 2, 3) and the deviations of lepton mixing angles ∆θ ij = θ ij − θ ij (for ij = 12, 13, 23), where m i 's and θ ij 's are the neutrino masses and mixing angles after breaking. It can be seen that in both NH and IH cases, because of the factor θ 13 in ∆θ 23 and the factor 1/ζ ∼ 5 or 1/ξ ∼ 30 in ∆θ 12 , we have |∆θ 23 | < |∆θ 13 | < |∆θ 12 | in general, barring the cases that accidental cancellations exist among A i 's. The suppression of ∆θ 23 also indicates that even with the breaking term in the (12) Table I we show the results for the other two breaking patterns in M D , namely, assigning breaking terms in the (22) and (32) positions of M D and resulting in the breaking scenarios of S2 and S3, respectively. We notice that the deviations of the neutrino mass matrix ∆M ν can also be expressed in terms of the parameters A i 's, except that the overall breaking parameters are modified to be B 22 = b/det(M R ) and B 32 = c * /det(M R ) for S2 and S3, respectively. We also observe that the analytic expressions for ∆θ ij 's and m i 's in S2 and S3 are quite similar, and in both scenarios there is no suppression factor of θ 13 in ∆θ 23 . As a result, one may expect larger deviation in θ 23 in S2 and S3 than that in S1. Thus, the last two breaking patterns may be distinguishable from the first one via the future precision measurement of θ 23 .
Having discussed some analytical results for the three breaking patterns in M D , we next turn to the detailed numerical analysis. On the one hand, the numerical analysis would extend the analysis to the scenarios where the entries in ∆M ν are not all real. On the other hand, we can also obtain the deviations on the Dirac CP-violating phase δ and the Majorana phases, which are not easy to obtain analytically. In the numerical analysis, for each breaking pattern we treat all the parameters in M D and M R as free parameters, and vary them within the same ranges as in the previous RG running study. For the breaking parameter , we vary it as ∈ [−1, 1]. The package Multinest is again employed to guide the parameter scan, and the same χ 2 function as before is utilized.
In Fig. 4 we show the numerical results in the case of NH for the three breaking patterns discussed above, i.e., S1 (left), S2 (middle) and S3 (right). Black/gray points have χ 2 < 30, and the red point in each case still denotes the best-fit scenario. For S1, S2 and S3, we obtain χ 2 min = 9.75, 0.03 and 0.58, respectively. Moreover, we only show the results that have δ < 0 in S2 and S3, as the results for δ > 0 are quite similar except for a sign change in δ . Lastly, we notice that for S2 and S3 there exist two branches of predictions, which are distinguished by black and gray points. From Fig. 4 we then observe: • According to the top three plots, we find that ∆θ 23 in S1 is less than one degree, much smaller than that in the other two breaking patterns. This numerical finding agrees with the analytical results in Table I, i.e., ∆θ 23 is suppressed by a factor of θ 13 in S1.
In the case with positive correlation δ and θ 23 deviate by almost the same amount, while |∆δ| is about three times larger than |∆θ 23 | for the case with negative correlation.
• From the plots in the second row of Fig. 4 we find that |∆θ 12 | is indeed larger than |∆θ 13 | and |∆θ 23 |, and it can reach around 15 • for all three breaking patterns. This is also in agreement with the analytical results given in Table I The numerical results for S1, S2 and S3 in the case of IH are shown in Fig. 5. Because the deviation in θ 12 now has an enhancement factor of 1/ξ ∼ 30, the numerical program is very sensitive to the initial value of θ 12 before breaking, so that locating the favored parameter space becomes challenging. The lowest values of χ 2 min that we are able to obtain are 28.84, 25.75 and 26.19 for S1, S2 and S3, respectively, and the corresponding scatter points are shown in red in Fig. 5. Moreover, in order to show larger region of parameter space, in this case of IH we require all scatter points (black) to satisfy χ 2 < 50. By inspecting the patterns in Fig. 5, we then observe: • In S1, θ 23 tends to be very close to 45 • , although the BF point has a large deviation in θ 23 , which may originate from some special combination in the input parameters.
The obtained δ , however, has a large spread in [−180 • , 180 • ). Regarding θ 12 and θ 13 , the deviation in θ 13 is quite small in general, while for θ 12 large deviations of O(10 • ) can be easily achieved.
entries in ∆M ν to be non-zero. For example, for S4 the obtained ∆M ν is given by where , and all A i 's are given in Eq. (29). With all non-zero entries in ∆M ν it is difficult to investigate the breaking effects on neutrino masses and lepton mixing angles analytically. Thus, we employ a similar numerical analysis as the previous breaking scenarios, and both the varied ranges of input parameters and the defined χ 2 function are kept to be the same.
In Fig. 6 we show the numerical results for S4 and S5 in NH. The obtained lowest values of χ 2 min are 9.12 and 1.75 for S4 and S5, respectively. Again, in both scenarios only the results that have δ < 0 are shown, as the other case of δ > 0 is quite similar except for a sign change in δ . From Fig. 6, we first notice that the patterns of the favored parameter space in S4 and S5 are quite similar, although in the latter case more extended parameter space is observed. Between θ 23 and δ positive correlations are identified when δ ∼ −90 • , and |∆θ 23 | is about four times smaller than |∆δ|. Such a positive correlation of ∆δ ∼ 4∆θ 23 is different from that in S2 and S3, where the positive correlation gives ∆δ ∼ ∆θ 23 , so that we may distinguish the breaking scenarios S4/S5 from S2/S3 by precisely measuring the correlation between θ 23 and δ in the upcoming experiments. Regarding θ 12 and θ 13 , it seems that ∆θ 13 is negatively correlated to ∆θ 12 , and |∆θ 13 | is about five times larger than |∆θ 12 |.
Lastly, as expected, the favored m ν is around 0.06 eV, and because the Majorana phase σ after breaking is still close to 90 • , we also have m ee ∼ 4 meV as in S2 and S3 under NH.
We now turn to the numerical results for S4 and S5 in IH, see Fig. 7 Above we have discussed various breaking scenarios of exact µ − τ reflection symmetry,     the left-and right-handed neutrinos, resulting in some particular forms of neutrino Dirac mass matrix M D and the Majorana mass matrix M R for the right-handed neutrinos. The obtained light neutrino mass matrix M ν is found to still possess the usual µ − τ reflection symmetry, which predicts maximal atmospheric mixing angle (θ 23 = 45 • ) and Dirac CP phase (δ = ±90 • ) along with the trivial Majorana phases. We later extend our study by incorporating the breaking of such symmetry, keeping in mind that theoretical as well as experimental results may favor non-maximal θ 23 .
The first possible breaking of the symmetry is due to the renormalization group running.
Here we choose the minimal supersymmetric standard model as our UV framework, and assume the symmetry to be exact at the GUT scale. When running towards low energies, we encounter three stages of running: above the two seesaw thresholds, between the thresholds, and lastly below the thresholds. Some noteworthy outcomes of our numerical RG analysis are summarized as follows: • The RG running between the thresholds is insignificant, as the two seesaw mass thresholds are found to be quite close. Such closeness of two thresholds is due to the fact that the two columns of the neutrino Yukawa matrix are related by the µ − τ reflection symmetry, particularly the symmetry on the right-handed neutrinos as proposed here.
• For both NH and IH scenarios, we find that the RG running effects above the seesaw thresholds are comparable to those below the thresholds. This would raise the necessity of considering RG running above the seesaw thresholds, if some flavor symmetry were imposed on the right-handed neutrino fields.
• For the three mixing angles, the deviations due to the RG running are all rather small, e.g., ∆θ 23 1 • , except that for θ 12 in IH can there exist a large deviation. The latter exception arises from the fact that the two light neutrino masses may cross each other, leading to an interchange of the order of two neutrino masses.
• The RG running effects of the Dirac and Majorana phases are also quite mild in NH, while large deviations of O(10 • ) can be observed in the case of IH.
• Lastly, we note that the known correlation between the positive/negative deviation of θ 23 and the neutrino mass hierarchies are again observed in this extended RG running above the seesaw thresholds.
Having shown that the RG running effects are quite mild, we then proceed to introduce explicit breaking terms in M D and M R , aiming at obtaining large deviations from the predictions of the exact µ − τ reflection symmetry. In total, we systematically investigate five possible breaking patterns, namely, assigning breaking terms in the (12), (22) and (32) positions of M D (denoted as S1, S2 and S3 breaking scenarios) and the (12) and (22) positions of M R (denoted as S4 and S5 breaking scenarios). Both analytical and numerical studies are pursued for S1, S2 and S3, while for S4 and S5 only the numerical results are attainable.
The main results of these breaking scenarios are listed as follows: • In NH we find that ∆θ 23 0.5 • in S1 while ∆θ 23 of a few degrees can be easily observed for the other breaking patterns. On the other hand, in IH all breaking patterns tend to have ∆θ 23 0.5 • , especially ∆θ 23 = 0 seems to hold exactly for S4 and S5.
• For the deviations in θ 13 , we obtain ∆θ 13 1 • for S2, S3 and S4 in NH, while for S1 and S5 deviations of a few degrees are possible. However, in the case of IH all breaking patterns tend to have small deviations in θ 13 , and again ∆θ 13 = 0 seems to hold exactly in S4 and S5 as well.
• The deviations in θ 12 are found to be around O(10 • ) in general, except that for S4 and S5 we observe ∆θ 12  Finally, we conclude the paper with the excitement and caution about the µ − τ reflection symmetry. Given the current status of neutrino oscillation data, the µ − τ reflection symmetry seems to stand out as a compelling reason for the bewildering flavor puzzles in the lepton sector. However, past two decades also witnessed the shift of the prevailing symmetry pattern in the lepton mixing matrix when more accurate neutrino oscillation data stepped in. Thus, along with continuing the pursuit of the implications behind the µ − τ reflection symmetry theoretically, one should also pay close attention to the experimental results in the upcoming years, especially from those measuring the value of θ 23 .
It is apparent that with → 0, the above M ν possesses the exact µ − τ reflection symmetry.
Next, we first perform a (23) rotation R 23 of π/4 on M ν , namely, with R 23 given by The phase in R T 23 M ν R 23 can be further removed by a diagonal phase matrix P φ = diag{i, 1, 1}, resulting in a pure real mass matrix as follows, .
Surprisingly, we then notice that the above matrix can be brought in a block diagonal form with a (13) rotation R 13 , whose mixing angle θ 13 coincides with the case without breaking! To be explicit, R 13 is given by It is then apparent that the above scenario corresponds to the IH case, as m 3 = 0, and the above matrix can be finally diagonalized by a (12) rotation. Since in the above diagonalization procedure we follow a (23)-(13)- (12) sequence that is same as the PDG convention, the mixing angles of θ 23 and θ 13 stay the same as the case without breaking. Note that if = 0 the above matrix is already diagonalized, and thus we have θ 12 = 0 without breaking; with breaking we instead require θ 12 = 0, so that θ 12 is not immune to the breaking in M R . On the other hand, we arrive at the NH case when 2b 2 − (c − d) 2 < 0, namely, where sgn(x) stands for the sign of x. Now because we need another (23) rotation, the final (23) rotation angle would deviate from π/4 when adopting the PDG convention, and in the meantime θ 13 would also get modified. As a result, in this NH case we expect that all the three mixing angles can be affected by the breaking terms in M R .
Appendix D: Details of best-fit scenarios in the RG running study In Table IV we provide the neutrino Yukawa matrix and the Majorana neutrino mass matrix at the high energy boundary for both the hierarchies. We also give the detailed numerical values of all the neutrino oscillation parameters at the various energy scales.