Integral solutions to the one-loop renormalization-group equations for lepton flavor mixing parameters and the Jarlskog invariant

Working in the basis where the charged-lepton Yukawa matrix is diagonal and making the $\tau$-dominance approximations, we analytically derive integral solutions to the one-loop renormalization-group equations (RGEs) for neutrino masses, flavor mixing angles, CP-violating phases and the Jarlskog invariant under the standard parametrization of the PMNS matrix in the standard model or its minimal supersymmetric extension for both Majorana and Dirac neutrinos. With these integral solutions, we carry out numerical calculations to investigate the RGE running of lepton flavor mixing parameters and the Jarlskog invariant, and also compare these integral solutions with the exact results obtained by numerically solving the one-loop RGEs. It is shown that these integral solutions coincide with the exact results and can well describe the evolution of lepton flavor mixing parameters and the Jarlskog invariant in most cases. Some important features of our integral solutions and the evolution behaviors of relevant flavor parameters are also discussed in detail both analytically and numerically.


Introduction
In the last two decades, compelling evidences obtained from a number of successful neutrino oscillation experiments have proved that neutrinos are massive and lepton flavor mixing exists [1], and this demonstrates that the standard model (SM) of particle physics is incomplete. In order to understand the origin of small neutrino masses and the observed lepton flavor mixing pattern, many models with extra flavor symmetries (see reviews [2][3][4][5][6] and references therein) or some new degrees of freedom (e.g., the right-handed neutrino fields in the type-I seesaw mechanism [7][8][9][10][11]) have been put forward at superhigh energy scales, as well as many phenomenological textures of lepton mass matrices (e.g., texture zeros of the neutrino mass matrix [6,[12][13][14]). With the help of the renormalization-group equations (RGEs), one can confront phenomenological consequences of those models or textures with current experimental data at the low energy scale. Thus, it is very important and useful to investigate the evolution of relevant flavor parameters or the stability of some specific textures against the energy scale by means of the RGEs, especially in the cases where nearly degenerate neutrino masses or large tan β in the minimal supersymmetric standard model (MSSM) is taken into account, so as to establish some correlations between physical phenomena at high and low energy scales and reveal some underlying structures of lepton mass matrices or flavor mixing pattern which are instructive for model building.
So far, whether neutrinos are Dirac or Majorana particles (which is referred to as the Dirac case or the Majorana case) is still an open question. If the neutrinoless double-beta (0νββ) decay is observed, we shall conclude that neutrinos are Majorana particles thanks to the Schechter-Valle theorem [15]. However, we can not claim that neutrinos are Dirac particles even though the 0νββ decay is not observed in experiments [6,16,17]. There still exist some rooms and interesting models for Dirac neutrinos, and it is worth considering these two possibilities before the nature of massive neutrinos is convincingly determined by future experiments. In the Majorana case, the small neutrino masses can be generated by the unique dimension-5 Weinberg operator [18] in an effective field theory. This operator can be obtained by integrating out heavy degrees of freedom in some extended models [19], such as the type-I seesaw mechanism. Then the one-loop RGE for the effective coupling matrix κ of Majorana neutrinos is given by [20][21][22][23] where t ≡ ln (µ/Λ EW ) with µ being an arbitrary renormalization scale between the electroweak scale Λ EW ∼ 100 GeV and the cutoff scale Λ, Y l is the charged-lepton Yukawa coupling matrix, and C = −3/2 , in SM 1 , in MSSM , α κ −3g 2 2 + 6y 2 t + λ , in SM −6g 2 1 /5 − 6g 2 2 + 6y 2 t , in MSSM (2) with g 2 , y t and λ being the SU (2) L gauge coupling, the top-quark Yukawa coupling and the Higgs self-coupling constant respectively. After spontaneous gauge symmetry breaking, the effective Majorana neutrino mass matrix is given by M ν = κv 2 (SM) or M ν = κv 2 tan 2 β/(1 + tan 2 β) (MSSM) with v 174 GeV being the vacuum expectation value of the SM Higgs field and tan β denoting the ratio of the vacuum expectation values of two Higgs doublets in the MSSM. In the Dirac case, the one-loop RGE for the Yukawa coupling matrix Y ν of Dirac neutrinos is [24][25][26][27] where C has been given in Eq. (2), and C = 3/2 , in SM 3 , in MSSM , α ν −9g 2 1 /20 − 9g 2 2 /4 + 3y 2 t , in SM −3g 2 1 /5 − 3g 2 2 + 3y 2 t , in MSSM (4) with g 1 being the gauge coupling. Note that in this case, the Yukawa coupling matrix Y ν must be much smaller so as to be accordant with the smallness of neutrino masses. Thus Y ν Y l holds pretty well and it is quite safe to ignore Y ν Y † ν in Eq. (3). The Dirac neutrino mass matrix is given by M ν = Y ν v (SM) or M ν = Y ν v tan β/ 1 + tan 2 β (MSSM). In the basis where Y l is diagonal, namely Y l = Diag{y e , y µ , y τ }, the neutrino mass matrix can be diagonalized by the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) lepton flavor mixing matrix U [28][29][30], i.e., U † M ν U * = D ν ≡ Diag{m 1 , m 2 , m 3 } in the Majorana case or U † H ν U = D 2 ν in the Dirac case with m i (for i = 1, 2, 3) being the neutrino masses and H ν being defined as H ν ≡ M ν M † ν . A popular parametrization of the PMNS matrix U is given by [1] in which c ij ≡ cos θ ij and s ij ≡ sin θ ij (for ij = 12, 13, 23) with θ ij lying in the first quadrant, δ is the so-called Dirac CP phase, P l ≡ Diag{e iφ e , e iφ µ , e iφ τ } with φ e , φ µ and φ τ being the unphysical phases associated with the charged-lepton fields, and P ν ≡ Diag{e iρ , e iσ , 1} with ρ and σ being the Majorana phases which become unphysical in the Dirac case. Based on Eqs. (1)-(4), the one-loop RGE running effects on neutrino masses and flavor mixing parameters have been extensively studied, which can be seen in the review [31] and references therein. With some specific parametrizations of the PMNS matrix U , the individual RGEs for neutrino masses, flavor mixing angles and CP-violating phases have been derived in Refs. [27,[32][33][34][35] and threshold effects in seesaw models have also been discussed, such as those in the type-I seesaw mechanism [36][37][38][39]. Furthermore, the running effects on leptonic CP-violating phases have been studied in detail [40][41][42][43][44], showing that the evolution of three CP-violating phases is entangled in the Majorana case so that the Dirac CP-violating phase can be radiatively generated even if it is initially assumed to be zero (in the Dirac case, the Dirac CP-violating phase keeps vanishing during the RGE evolution if it is initially zero). And in particular, some recent works [45][46][47] find that the RGE effects can play a significant role in establishing a direct link between the low energy CP violation and the CP-violating asymmetries at a superhigh energy scale, and making the leptogenesis mechanism [48] work successfully.
In the present work, our main purpose is to analytically derive integral solutions to the one-loop RGEs for neutrino masses, flavor mixing angles, CP-violating phases and the Jarlskog invariant of U both in the Majorana case and in the Dirac case, with the τ -dominance approximations. Some previous attempts in this connection have been made to some extent [49][50][51][52][53][54][55][56]. But our work differs from them at least in the following aspects: • Ref. [49] mainly focuses on the seesaw threshold effects in the low-scale seesaw model, and has only derived analytical results for neutrino masses and flavor mixing angles with a special parametrization of U by assuming that CP is preserved. While in our work, we consider the most general case with the popular parametrization of U given in Eq. (5) below the cutoff scale, and the integral results for neutrino masses, flavor mixing angles, CP-violating phases and the Jarlskog invariant are exhaustively derived. In particular, the case for Dirac neutrinos is also taken into account in this work.
• Unlike those works done in Refs. [50][51][52][53]55,56] where some special situations are considered, such as the µ-τ reflection symmetry leading to θ 23 = π/4 and δ = ±π/2 [51,57] or the lightest neutrino being massless, our work is essentially independent of the specific textures of lepton Yukawa coupling matrices or flavor mixing patterns, and gives the most general results for relevant flavor parameters without further assumptions. Thus, the corresponding results in Refs. [50][51][52][53]55,56] can be easily reproduced from our results under some further constraints on flavor mixing parameters or neutrino masses.
• The general integral results for the Jarlskog invariant in the Majorana and Dirac cases have been acquired in Ref. [54] without any specific parametrization of U . Instead, our work takes the widely used parametrization of U given in Eq. (5) so that one can clearly see properties of the evolution of the Jarlskog invariant from one scale to another, especially its dependence on flavor mixing angles and CP-violating phases. In the Majorana case, our result clearly shows that there is an additional term mainly dominated by two Majorana CP-violating phases, from which the Jarlskog invariant can be radiatively generated even if it is initially vanishing at a specific energy scale. This observation is very intuitive to understand CP violation in a long-baseline neutrino oscillation experiment, whose strength is uniquely governed by the Jarlskog invariant.
• Moreover, in our work, we discuss the evolution behaviors of lepton flavor parameters and the Jarlskog invariant both analytically and numerically in great detail by using the latest experimental data and global-fit inputs, including the T2K measurement of CP violation [58]. This is also a merit of our work.
Compared with the differential forms of RGEs for lepton flavor parameters, the integral solutions can explicitly reveal their evolution behaviors, since they are only dominated by two quantities I β (for β = κ or ν) and ∆ τ which are integrals of α β (for β = κ or ν) and y 2 τ to the energy scale and almost independent of the initial inputs. Therefore, given the values of these two quantities against the energy scale, we can easily obtain the RGE running of relevant flavor parameters with some inputs. That is why our integral solutions are expected to be very useful to study radiative corrections to some interesting flavor mixing patterns or textures of lepton mass matrices, shed light on some possible underlying flavor symmetries at a superhigh energy scale, and establish a direct link between physical phenomena at the low and high energy scales.
The rest part of this paper is organized as follows. In section 2, the integral solutions to RGEs for lepton flavor mixing parameters and the Jarlskog invariant in the SM or MSSM for both Majorana and Dirac neutrinos are analytically derived. The numerical calculations are carried out to illustrate the evolution behaviors of relevant flavor parameters in section 3. We summarize our main results in section 4.

Integral solutions to one-loop RGEs
Working in the basis where Y l is diagonal, we find that Y l remains diagonal 1 during the RGE running so that we can integrate Eq. (1) or (3) from an arbitrary lower energy scale µ to the superhigh energy scale Λ and obtain or corresponding to Majorana neutrinos or Dirac neutrinos, where T l = Diag{I e , I µ , I τ }, and I β (for β = κ, ν) and I γ (for γ = e, µ, τ ) are defined as Considering y 2 e y 2 µ y 2 τ in the SM or MSSM together with our requirement of tan β 30 and the smallness of the loop factor 1/ (16π 2 ), it is reasonable and safe to make the τ -dominance approximations, i.e., I e I µ 1 and I τ 1 + ∆ τ with being a small quantity. Since the RGE evolution of gauge couplings, Yukawa coupling matrices of charged fermions and the Higgs self-coupling constant are essentially independent of κ in the Majorana case and Y ν in the Dirac case, the values of I β (for β = κ, ν) and ∆ τ against the energy scale µ do not depend on initial inputs from the neutrino sector, as illustrated in Fig. 1. As can be seen from Fig. 1, the magnitude of ∆ τ is of O (10 −5 ) at Λ EW in the SM, and it can be largely 1 Y l keeps diagonal strictly in the Majorana case but approximately in the Dirac case since in the Dirac case, the RGE of Y l contains Y ν Y † ν which can make Y l deviate from the diagonal form during the RGE running. Fortunately, due to Y ν Y l in this case, the off-diagonal elements induced by RGE effects are much smaller than the diagonal ones of Y l , namely, Y l approximately remains diagonal. The complete one-loop RGEs of Yukawa coupling matrices, gauge couplings and the Higgs self-coupling constant in the Majorana and Dirac cases can be found in the latest review [6].  Figure 1: The values of I β (for β = κ or ν) and ∆ τ against the energy scale µ in the Majorana case or the Dirac case with Λ = 10 14 GeV, where SM and MSSM with tan β = 10 or tan β = 30 are considered. "×100" in the right two panels means that the value of ∆ τ in the SM has been enhanced by a factor 100. enhanced in the MSSM with a large value of tan β. It can reach O (10 −3 ) with tan β = 10 or O (10 −2 ) with tan β = 30 at Λ EW . Thus generally, ∆ τ can be regarded as a small quantity in the SM or MSSM with tan β 30. Note that the values of ∆ τ in the Majorana and Dirac cases are practically equal. The reason is that after neglecting Y ν Y † ν in the Dirac case, the one-loop RGEs for charged-fermion Yukawa coupling matrices and gauge couplings are the same in the Majorana and Dirac cases, and with the same initial inputs for these parameters, they evolute equally in these two cases. For the same reason, values of I κ within the MSSM in the Majorana case are the square of the corresponding values of I ν in the Dirac case, as can be observed in Fig. 1.
By the way, at the one-loop level, there are some interesting and exact relations which can be derived from Eqs. (6)- (7). That is, in the Majorana case [34,54], and in the Dirac case [52], where ∆m 2 ij ≡ m 2 i − m 2 j (for i, j = 1, 2, 3) are neutrino mass-squared differences and J is the Jarlskog invariant of CP violation [59], defined as with the Greek and Latin subscripts running over (e, µ, τ ) and (1, 2, 3), respectively. With the parametrization of the PMNS matrix U given in Eq. (5), the Jarlskog invariant is given by J = 1/8 sin 2θ 12 sin 2θ 13 cos θ 13 sin 2θ 23 sin δ. Those relations in Eqs. (10) and (11) are very interesting because they connect neutrino masses, phases in U or the Jarlskog invariant at µ and Λ with each other via some simple ways without any approximation at the one-loop level. The relationships for neutrino masses in the Majorana and Dirac cases are slightly different, but both of them explicitly show that if one of the neutrinos is massless, it will remain massless at the one-loop level. Only after the two-loop effects are taken into account, can nonzero neutrino mass be generated for the initially massless neutrino [56,60]. The relationship for det [U ] in the Majorana case gives us a correlation between phases in U , but such a correlation is dependent on the parametrization of U . Within the parametrization given in Eq. (5), φ e + φ µ + φ τ + ρ + σ µ − φ e + φ µ + φ τ + ρ + σ Λ = 0 or π holds, and a similar relation φ e + φ µ + φ τ + ρ + σ − φ µ − φ e + φ µ + φ τ + ρ + σ − φ Λ = 0 or π holds if the parametrization proposed in Ref. [61] and the phase matrices P l and P ν are used. The differential forms of these relations for phases can be found in Refs. [35,44]. The relationship for the Jarlskog invariant in the Dirac case is instructive. It transparently shows that if J = 0 holds initially, it will keep vanishing during the RGE running, implying that in the Dirac case CP violation at a low energy scale can not be radiatively generated when there is no CP violation at the superhigh energy scale and vice versa, and if CP violation does exist, it will exist at any energy scale below the cutoff scale. Now we are going to perturbatively solve Eqs. (6) and (7) with the parametrization of U given in Eq. (5) 2 . Note that in the Majorana case, three unphysical phases φ α (for α = e, µ, τ ) in P l are all involved in the RGE running and have their own evolutions against the energy scale µ, but in the Dirac case, two phases in P ν and one overall phase in P l can be cancelled in H ν so that only two of the five unphysical phases take part in the RGE running. Thus in the Dirac case, we redefine P l as P ≡ Diag e iφ 1 , e iφ 2 , 1 and take P ν = 1. It is convenient to define with ij = 12, 13, 23 and α = e, µ, τ in the Majorana case, or with ij = 12, 13, 23 and k = 1, 2 in the Dirac case, to describe the evolution of the corresponding parameters. Due to the smallness of above quantities defined in Eq. (13) or (14) which are proportional to ∆ τ , it is reasonable to take them and ∆ τ as small perturbations in the leading order approximation. To make the relevant expressions concise, we can define ∆U ≡ U (µ) − U (Λ) which is a function of the quantities defined in Eq. (13) or (14) and satisfies at the leading order guaranteed by the unitarity of U (µ) and U (Λ). After expanding U (Λ) with respect to quantities defined in Eq. (13) or Eq. (14), we can obtain the explicit expression of ∆U in terms of the low-energy parameters and quantities defined above, as shown in Appendix A.
• All the results for ∆θ ij (for ij = 12, 13, 23), ∆δ, ∆ρ and ∆σ are proportional to ∆ τ so the relevant parameters at µ involved in these results can generally be replaced by their values at Λ at the leading order level. The same observation is also true for terms proportional to ∆ τ in the results for m i (for i = 1, 2, 3) and J .
• Since the signs of ∆ τ are opposite in the SM and MSSM which can be seen from Eq. (9) and Fig. 1, the running directions of flavor mixing angles and CP-violating phases, whose evolutions are determined by Eqs. (20)- (22), are opposite in the SM and MSSM.
Though it is not obvious that Eq. (32) coincides with the relation for the Jarlskog invariant given in Eq. (11), one may use Eq. (29) to check that the relation for J in Eq. (11) actually leads to the results given by Eqs. (32) and (33). All analytical results for flavor mixing angles, the Dirac CP-violating phase and the Jarlskog invariant in the Dirac case are much simpler than those in the Majorana case, owing to two less physical degrees of freedom in the former case. Some comments are as follows: • Similar to that in the Majorana case, ∆θ ij (for ij = 12, 13, 23) and ∆δ are all proportional to ∆ τ , and thus the running directions of these parameters are opposite in the SM and MSSM. Generally, relevant flavor parameters at µ involved in these results can be replaced by their values at Λ, as well as those contained in the terms proportional to ∆ τ in the expressions of m i (for i = 1, 2, 3) and J .
• As shown in Eqs. (31) and (32), ∆δ ∝ sin δ and J (µ) ∝ J (Λ) hold. It means that δ and J can not be radiatively generated if they are initially vanishing. This observation differs from that in the Majorana case.
• Taking all the CP-violating phases to be vanishing and with the replacements ζ −1 ij → ξ ij (i.e., m i → m 2 i for i, j = 1, 2, 3) in the Majorana case, and meanwhile taking the Dirac CPviolating phase to be zero in the Dirac case, one may check that the corresponding results for flavor mixing angles in Eqs. (20) and (30) are exactly the same. It is easy to understand these results, since after unphysical phases being rotated away from M ν and H ν in this case, M ν in the Majorana case and H ν in the Dirac case are both real and can be diagonalized by a real orthogonal matrix, so all things are formally the same in the Majorana and Dirac cases except the eigenvalues of M ν and H ν .

Numerical analysis and discussion
In the numerical analysis, we use the best-fit values of neutrino parameters obtained from the latest global analysis of currently available neutrino oscillation data [62,63] including the T2K measurements of the Dirac CP-violating phase [58], namely, and where both the normal neutrino mass ordering (upper values) and the inverted one (lower values) have been taken into account, and the two neutrino mass-squared differences are defined as δm 2 ≡ m 2 2 − m 2 1 and ∆m 2 ≡ m 2 3 − (m 2 1 + m 2 2 ) /2. These experimental data are all given at the eletroweak scale Λ EW ∼ 100 GeV. We consider the following two neutrino mass spectra: • The normal mass ordering (NMO) with m 1 (Λ EW ) = 0.001 eV. In the Dirac case, ξ ij (µ) ξ ij (Λ EW ) 1 (for ij = 21, 31, 32) hold pretty well, and thus the results for flavor mixing angles and the Dirac CP-violating phase can be largely simplified and their evolution behaviors are more transparent. But similar approximations for ζ ij (µ) (for ij = 21, 31, 32) are not good in the Majorana case.
• The inverted mass ordering (IMO) with m 3 (Λ EW ) = 0.001 eV. In the Majorana case, Here, we do not consider the case of nearly degenerate neutrino masses in which ζ −1 21 (µ) or ξ 21 (µ) is so strongly enhanced that ∆θ 12 , ∆δ, ∆ρ and ∆σ containing ζ −1 21 (µ) or ξ 21 (µ) are significantly enlarged and the approximations that they are small quantities become bad especially in the MSSM with a large tan β. Hence the analytical results for ∆θ 12 , ∆δ, ∆ρ and ∆σ can remarkably deviate from the corresponding exact results in this case. Actually, we also confront this situation in the MSSM with a sizeable tan β for the inverted neutrino mass spectrum but it is not severe and thus acceptable when tan β 30. It is worth mentioning that the normal neutrino mass ordering is currently favored over the inverted one at the level of around 3σ indicated by a globle analysis of current neutrino oscillation data and the total neutrino mass is constrained to be m ν < 0.12 eV by some cosmology observations [62][63][64][65]. The latter infers that nearly degenerate neutrino masses are disfavored by the CMB anisotropies at 2.4σ level or at 5.9σ level after the BAO data are added [66]. In the numerical analysis, we only exhibit the numerical results in the MSSM with tan β = 10 or 30 and neglect those in the SM because in the SM, the RGE effects are extremely small and hence all analytical results coincide with the exact ones very well even in the case where neutrino masses are nearly degenerate. But this does not mean that these small RGE-induced effects in the SM are inessential, on the contrary, sometimes they can play a greatly important role, such as establishing a direct connection between the CP-violating asymmetries at a superhigh energy scale and CP violation at the electroweak scale via the seesaw bridge [45,47].
To compute the evolution of flavor mixing parameters and the Jarlskog invariant from Λ down to µ, we choose the initial values of relevant parameters at Λ in such a way that the best-fit values of θ 12 , θ 13 , θ 23 , δm 2 and ∆m 2 at Λ EW shown in Eqs. (34) and (35) can be satisfied, and some given values of δ (together with ρ and σ in the Majorana case) and the lightest neutrino mass at Λ or Λ EW are initially input or can be achieved. Therefore the initial inputs at Λ may not be the same in different cases. With these initial inputs, we calculate both the exact results by numerically solving the RGEs and the approximate ones by means of the analytical results we have obtained above.

Neutrino masses and flavor mixing angles
First, we compute the evolution of neutrino masses and flavor mixing angles. We require that m 1 (Λ EW ) = 0.001 eV in the NMO case or m 3 (Λ EW ) = 0.001 eV in the IMO case be satisfied, δ (Λ EW ) take its best-fit value shown in Eq. (34), and ρ (Λ) = σ (Λ) = 0 be initially input in the Majorana case. The results for neutrino masses and flavor mixing angles are illustrated in Figs. 2 and 3, respectively. In particular, the values of neutrino masses and flavor mixing angles at Λ EW obtained with the help of the analytical expressions derived in section 2 are explicitly listed in Table 1, where the corresponding numbers shown in the parentheses are the relative errors compared to the exact results acquired by numerically solving the one-loop RGEs.
From Fig. 2 and Table 1, one can see that the approximate results for neutrino masses are consistent with the exact ones very well in all the cases we have considered. The evolution of neutrino masses is dominated by m i (µ) I β m i (Λ) (for i = 1, 2, 3 and β = κ or ν). Considering the values of I κ and I ν against µ shown in Fig. 1, it is easy to understand that the evolution of neutrino masses in the Majorana case (or that with tan β = 30) is slightly severer than that in the Dirac case (or that with tan β = 10). And the running behaviors of neutrino masses in the SM are more violent than those in the MSSM, especially with Majorana neutrinos, indicated by I β (for  β = κ or ν) shown in Fig. 1.
As for the evolution of θ ij (for i = 12, 13, 23) shown in Fig. 3 and Table 1, the relative errors of approximate results for them are around 1% or smaller at Λ EW , except those for ∆θ 12 with tan β = 30 in the IMO case and ∆θ 13 in the Majorana case with tan β = 30 and the normal neutrino mass ordering. Some discussions and comments on the evolution behaviors are as follows: • Comparing the corresponding results with tan β = 10 and tan β = 30, such as those for Dirac neutrinos with the normal mass ordering in the MSSM where tan β = 10 and tan β = 30 are respectively considered, the results for ∆θ ij (for ij = 12, 13, 23) with tan β = 30 are about ten times larger than those with tan β = 10 during the RGE running, mainly due to ∆θ ij ∝ ∆ τ (for ij = 12, 13, 23) which are approximately proportional to tan 2 β.
• As can be seen from the left upper and lower panels of Fig. 3, the running of ∆θ 12 can be largely enhanced in the IMO case, but the running direction keeps unchanged compared to that in the NMO case. The reason is that the expression for ∆θ 12 given in Eq. (20) or Eq. (30) contains ζ −1 21 or ξ 21 whose value can be strongly enhanced in the IMO case. To make it more distinct, we notice that ∆θ 12 roughly approximates to ∆θ 12 −1/2∆ τ ζ −1 21 (or ξ 21 ) sin 2θ 12 sin 2 θ 23 in the Majorana (or Dirac) case, where terms contain sin θ 13 have been neglected owing to the smallness of θ 13 . It is clear that ∆θ 12 is always positive and enhanced by the largeness of ζ −1 21 or ξ 21 in the IMO case. Especially, when tan β = 30 is taken in the IMO case, ∆θ 12 is extremely enlarged and no longer a small quantity, thus the approximate results for ∆θ 12 in those cases severely deviate from the exact ones as shown in the left lower panel of Fig. 3 and Table 1.

Majorana
• To understand the evolution of ∆θ 13 shown in the middle panels of Fig. 3, one may consider the neutrino mass spectra and the smallness of θ 13 to simplify the analytical expression for θ 13 as ∆θ 13 1/2∆ τ sin 2θ 13 cos 2 θ 23 in the IMO case for both Majorana and Dirac neutrinos, and ∆θ 13 −1/2∆ τ sin 2θ 13 cos 2 θ 23 in the NMO case for Dirac neutrinos, where the sign difference is induced by different signs of ζ 3i or ξ 3i (for i = 1, 2) in the normal and inverted neutrino mass ordering cases. It is apparent that the evolution of ∆θ 13 is suppressed by the Table 1: The values of m i (Λ EW ) /m i (Λ) (for i = 1, 2, 3) and ∆θ ij (Λ EW ) (for ij = 12, 13, 23) obtained by means of the analytical expressions derived in section 2 in the MSSM with tan β = 10 or 30. The initial inputs at Λ = 10 14 GeV are chosen to guarantee that m 1 = 0.001 eV or m 3 = 0.001 eV and the best-fit values given in Eqs. (34) and (35) can be achieved at Λ EW 100 GeV, and in addition, ρ = σ = 0 are initially input at Λ in the Majorana case. The corresponding numbers given in the parentheses are the relative errors compared to the exact results obtained by numerically solving the one-loop RGEs. smallness of θ 13 , and the running directions are opposite in the normal and inverted neutrino mass ordering cases but the absolute values of theirs are nearly equal for Dirac neutrinos, as shown in Fig. 3 and Table 1. Additionally, indicated by the approximate analytical results, the corresponding values for ∆θ 13 in the IMO case for Majorana and Dirac neutrinos with the same tan β are negative and roughly equal, and this can also be transparently seen in the middle-lower panel of Fig. 3 and Table 1. The evolution of ∆θ 13 in the NMO case for Majorana neutrinos is relatively exotic since its value seems to be largely reduced compared to other cases, as can be seen in the middle-upper panel of Fig. 3. In fact, carefully checking the analytical expression of ∆θ 13 given in Eq. (20), one may discover that there is a large cancellation in this case, and this is also the reason why the relative error of the approximate result for ∆θ 13 is large and reaches −4% with tan β = 10 or −8.6% with tan β = 30 at Λ EW .
• The analytical results for ∆θ 23 are consistent with the exact ones very well in all the cases under consideration and it is quite easy to understand its evolution behaviors in different cases, as shown in the right panels of Fig. 3. After the specific neutrino mass spectra are taken into account, the analytical expressions for ∆θ 23 given in Eqs. (20) and (30) can be further simplified to ∆θ 23 −1/2∆ τ sin 2θ 23 (or ∆θ 23 1/2∆ τ sin 2θ 23 ) in the normal (or inverted) neutrino mass ordering case. Therefore, ∆θ 23 is positive in the NMO case and negative in the IMO case, but its absolute values for these two neutrino mass spectra with the same tan β are roughly equal not only in the Majorana case but also in the Dirac case, as shown in the right upper and lower panels of Fig. 3 and Table 1. And with the same neutrino mass spectrum and tan β, the values of ∆θ 23 in the Majorana and Dirac cases are also nearly equal. As indicated by the right upper panel of Fig. 3 and Table 1, the results for ∆θ 23 in the NMO case for Majorana and Dirac neutrinos slightly depart from each other to some extent. It is mainly because in the NMO case for Majorana neutrinos, ζ −1 32 1 is not a good approximation.
From the above results and discussions, it is interesting to see that the evolution behaviors of neutrino masses and flavor mixing angles in the Majorana case even with initially vanishing Majorana CP-violating phases (i.e., ρ (Λ) = σ (Λ) = 0) can be distinguished from those in the Dirac case by their RGE running strengths [41]. One may also consider the initially nonvanishing Majorana CP-violating phases in the Majorana case, and find that the cancellation in ∆θ 13 is weakened or disappears in the NMO case.

CP-violating phases and the Jarlskog invariant
The evolution behaviors of CP-violating phases and the Jarlskog invariant in the Majorana case are much more complicated than those in the Dirac case as shown in Eqs. (20) and (30), since in the Majorana case there are two additional Majorana CP-violating phases, whose evolution behaviors are entangled with that of the Dirac CP-violating phase. Thus as discussed in section 2, the Dirac CP-violating phase δ and the Jarlskog invariant J in the Majorana case can be radiatively generated via the one-loop RGE running unless both ρ and σ are initially equal to 0 or π/2, while those in the Dirac case can not be radiatively generated by means of the one-loop RGE running. This is important and intuitive for us to distinguish between the evolution of relevant flavor parameters in the Majorana and Dirac cases, and understand CP violation at the eletroweak scale. In this subsection, we are going to compare the evolution of δ and J in the Majorana and Dirac cases, and discuss the entanglements of three CP-violating phases in the Majorana case. For our purposes, we choose some sets of initial inputs in the Majorana case, to guarantee that at the electroweak scale Λ EW , m 1 = 0.001 eV or m 3 = 0.001 eV, the best-fit values of two neutrino mass-squared differences and three flavor mixing angles can be achieved, and they satisfy one of the following requirements: • δ (Λ EW ) takes its best-fit value given in Eq. (34), and ρ (Λ) = σ (Λ) = 0 are initially input; • σ (Λ EW ) = π/4 can be achieved, and δ (Λ) = ρ (Λ) = 0 are initially input; • ρ (Λ EW ) = π/4 can be achieved, and δ (Λ) = σ (Λ) = 0 are initially input.
In the Dirac case, we only consider the initial inputs, from which m 1 (Λ EW ) = 0.001 eV or m 3 (Λ EW ) = 0.001 eV can be satisfied, and two neutrino mass-squared differences, three flavor mixing angles and the Dirac CP-violating phase can take their best-fit values given in Eqs. (34) and (35) at Λ EW . With these initial inputs, the results for CP-violating phases and the Jarlskog invariant are plotted in Fig. 4 for Majorana neutrinos and in Fig. 5 for Dirac neutrinos. For illustration, the values of ∆δ, ∆ρ, ∆σ and J at Λ EW in the Majorana case obtained with the help of Eqs. (21)- (25) are listed in Table 2, and those of ∆δ and J at Λ EW in the Dirac case obtained by the aid of Eqs. (31)-(33) are listed in Table 3. Again, the corresponding numbers given in the parentheses are the relative errors compared to the exact results obtained by numerically solving the one-loop RGEs. At the first sight of Figs. 4 and 5 together with Tables 2 and 3, one can find that results in the NMO case coincide with the exact results pretty well, no matter which type of neutrinos (Majorana or Dirac) is and how large the value of tan β (10 or 30) is, and the relative errors are of O (1%) or smaller at Λ EW . Since generally both the CP-violating phases and the Jarlskog invariant contain ζ −1 21 or ξ 21 which can be largely enlarged in the IMO case, the results in the IMO case are worse, especially those with tan β = 30 in the Majorana case, but the relative errors of results with tan β = 10 in the Majorana case or those in the Dirac case are mostly O (10%) or smaller at Λ EW , which are essentially acceptable. To understand the evolution behaviors of these CP-violating phases and the Jarlskog invariant, we first make some further approximations for those analytical formulas by taking into account the neutrino mass spectra and initial inputs under consideration, as well as the smallness of θ 13 . Though these simplified versions may not be consistent with the original ones very well, they can illustrate the salient properties of the evolution. We begin with those in the Majorana case.
(3) δ (Λ) = σ (Λ) = 0 Table 3: The values of ∆δ (Λ EW ) and J (Λ EW ) obtained by means of the analytical expressions derived in section 2 for Dirac neutrinos in the MSSM with tan β = 10 or 30. The initial inputs at Λ = 10 14 GeV are chosen to guarantee that m 1 = 0.001 eV or m 3 = 0.001 eV and the best-fit values given in Eqs. (34) and (35) can be achieved at Λ EW 100 GeV. The corresponding numbers given in the parentheses are the relative errors compared to the exact results obtained by numerically solving the one-loop RGEs.
Normal neutrino mass ordering (NMO) Inverted neutrino mass ordering (IMO) tan β = 10 tan β = 30 tan β = 10 tan β = 30 With ρ (Λ EW ) = π/4, Eqs. (21)- (22) can be simplified to in the normal neutrino mass ordering case; and in the inverted neutrino mass ordering case. In these two cases, the Jarlskog invariant is also determined by Eq. (42) in the normal neutrino mass ordering case; or ∆δ −∆ τ ξ 21 sin θ 13 sin 2θ 23 sin 2θ 12 sin δ , in the inverted neutrino mass ordering case, where the Jarlskog invariant is governed by Then based on Eqs. (36)- (46), some discussions and comments on the evolution behaviors of these CP-violating phases and the Jarlskog invariant shown in Figs. 4 and 5 or Tables 2 and 3 are in order.
• Since the RGE running effects are dominated by ∆ τ , which is approximately proportional to the value of tan 2 β, in general the results with tan β = 30 are nearly ten times larger than the corresponding results with tan β = 10 as can be seen from Figs. 4 and 5 or Tables 2  and 3. But this does not seem to be true for the results with ρ (Λ) = σ (Λ) = 0 in the IMO case. The formulas in Eq. (38), suppressed by sin θ 13 , can well describe the evolution of ∆δ, ∆ρ and ∆σ with tan β = 10 in this case, but when tan β = 30, an additional term proportional to ∆ τ ζ −1 21 sin 2 (ρ − σ) gradually dominates the evolution of ∆δ, ∆ρ and ∆σ during the running. Since in the IMO case, ζ −1 21 is extremely large and sin 2 (ρ − σ) at low energy scales with tan β = 30 is not very small, this additional term becomes dominant over the evolution. In addition, it is the largeness of ζ −1 21 for Majorana neutrinos or ξ 21 for Dirac neutrinos in the IMO case, together with the huge value of tan β dominating the strength of RGE running, that makes our approximations worse. Thus the results obtained from the analytical expressions with tan β = 30 in the IMO case remarkably deviate from the corresponding exact results, as obviously shown in Figs. 4 and 5 or Tables 2 and 3.
• Comparing the results with the inputs ρ (Λ) = ρ (Λ) = 0 in the Majorana case shown in the first row of Fig. 4 to those in the Dirac case illustrated in Fig. 5, one can find that both the running direction and strength of ∆δ are different. These differences between Majorana and Dirac neutrinos in the IMO case are obviously shown by Eqs. (38) and (46). However, the differences in the NMO case are not obviously indicated by Eqs. (36) and (45), and each term in Eqs. (36) and (45) needs to be carefully checked and compared. Nevertheless, one can conclude that the evolution behaviour of the Dirac CP-violating phase δ with the initially vanishing Majorana CP-violating phases in the Majorana case can be distinguished from that in the Dirac case not only by its RGE running strength but also by its running direction. Additionally, in the Majorana case, the Majorana CP-violating phases can be radiatively generated even they are initially vanishing, namely ρ (Λ) = σ (Λ) = 0. The evolution behaviors of the Jarlskog invariant J in the Majorana and Dirac cases are similar but the running effect on J with the inverted neutrino mass ordering in the Majorana case is stronger than that in the Dirac case due to ζ −1 21 > ξ 21 , as indicated by Eqs. (39) and (46) and clearly shown in Figs. 4 and 5.
• In Fig. 4 and Table 2, we have shown the results with different inputs for CP-violating phases in the Majorana case. It is evident that these three CP-violating phases are entangled with one another during the RGE running, implying that once there is a nonvanishing phase initially, the other two phases may be generated radiatively via the one-loop RGE running. If one wants to gain strong running effects to generate a sizeable phase, usually a large tan β and the nearly degenerate or inverted neutrino mass hierarchy should be taken into consideration [40,44], but under this circumstance, our analytical results are poor to describe the evolution of CP-violating phases. With the help of Eqs. (40)- (44), it is easy to understand the evolution behaviors of three CP-violating phases and the Jarlskog invariant with the Dirac CP-violating phase vanishing initially (i.e., δ (Λ) = 0), shown in the last two rows of Fig. 4. One may check that the running directions of three CP-violating phases in the case of δ (Λ) = ρ (Λ) = 0 are all opposite to those in the case of δ (Λ) = σ (Λ) = 0. And in the former case, ∆δ < 0, ∆ρ > 0 and ∆σ > 0 hold for both the normal and inverted neutrino mass hierarchies. As for the running strengths of these phases, Eqs. (40)- (44) indicate that in the IMO case, the absolute values of ∆δ, ∆ρ and ∆σ with δ (Λ) = ρ (Λ) = 0 are nearly equal to those with δ (Λ) = σ (Λ) = 0; but in the NMO case, the former ones are slightly larger than the latter ones since there is an additional term proportional to ζ 32 − ζ −1 32 in the formulas given by Eq. (40) compared with those in Eq. (43), which can enlarge the absolute values of ∆δ, ∆ρ and ∆σ. These properties are explicitly illustrated by Fig. 4 and Table 2. Actually, the sizes of ∆δ, ∆ρ and ∆σ in the same case can also be understood by taking advantage of the analytical results given in Eqs. (40)- (44). For example, in the NMO case with δ (Λ) = σ (Λ) = 0, Eq. (43) tells us that the absolute value of ∆δ is the largest among those of ∆δ, ∆ρ and ∆σ, which are different from one another only by some overall factors (namely 1 for ∆δ, cos 2 θ 12 for ∆ρ and sin 2 θ 12 for ∆σ). Finally, the evolution of J in the case of δ (Λ) = 0 is controlled by that of ∆δ, namely J (µ) ∝ ∆δ, as shown in Eq. (42), indicating that the evolution of J is similar to ∆δ in this case, and both δ and J can be radiatively generated by the one-loop RGEs even if they are initially vanishing. Before ending this section, it is worth remarking that in order to investigate general evolution behaviors of lepton flavor mixing parameters and the Jarlskog invariant by making full use of the analytical results obtained in section 2, we have only considered the initial inputs which are essentially independent of some specific textures of lepton Yukawa coupling matrices or flavor mixing patterns. These analytical results can certainly be applied to some fantastic models or flavor mixing patterns, such as the well-known tri-bimaximal mixing pattern [67][68][69], democratic mixing pattern [70][71][72] and the µ-τ reflection symmetry [51,57], which have been previously studied either with the differential RGEs [73][74][75][76] or with the special forms of their integral solutions [51][52][53]55].

Summary
The RGEs as a powerful tool to establish a link between physical phenomena at the high and low energy scales have been extensively studied. In this work, working in the basis where the charged-lepton Yukawa matrix Y l is diagonal and with the standard parametrization of the PMNS Table 4: The values of m i (Λ EW ) /m i (Λ) (for i = 1, 2, 3) and ∆θ ij (Λ EW ) (for ij = 12, 13, 23) obtained by means of the analytical expressions derived in section 2 in the MSSM with tan β = 10 or 30. The initial inputs at Λ = 10 14 GeV are chosen to guarantee that m 1 = 0.03 eV or m 3 = 0.015 eV and the best-fit values given in Eqs. (34) and (35) can be achieved at Λ EW 100 GeV, and in addition, ρ = σ = 0 are initially input at Λ in the Majorana case. The corresponding numbers given in the parentheses are the relative errors compared to the exact results obtained by numerically solving the one-loop RGEs.
Normal neutrino mass ordering (NMO) Inverted neutrino mass ordering (IMO) tan β = 10 tan β = 30 tan β = 10 tan β = 30 Majorana neutrinos matrix U given in Eq. (5), we have analytically derived integral solutions to the one-loop RGEs for neutrino masses, flavor mixing angles, CP-violating phases and the Jarlskog invariant in the SM or MSSM for both Majorana and Dirac neutrinos. In addition, we have also gained some interesting and exact relations between neutrino masses or phases at the superhigh energy scale Λ and those at a lower energy scale µ for Majorana neutrinos at the one-loop level, as well as the reltions between neutrino masses or the Jarlskog invariants at Λ and µ for Dirac neutrinos.  ), ∆ρ (Λ EW ), ∆σ (Λ EW ) and J (Λ EW ) obtained by means of the analytical expressions derived in section 2 for Majorana neutrinos in the MSSM with tan β = 10 or 30. The initial inputs at Λ = 10 14 GeV are chosen to guarantee that m 1 = 0.03 eV or m 3 = 0.015 eV and the best-fit values of two neutrino mass-squared differences and three flavor mixing angles given in Eqs. (34) and (35) can be achieved at Λ EW 100 GeV, and in addition, three CP-violating phases are required to be (1) δ (Λ EW ) = 1.28π or 1.52π in the normal or inverted neutrino mass ordering case and ρ (Λ) = σ (Λ) = 0; (2) σ (Λ EW ) = π/4 and δ (Λ) = ρ (Λ) = 0; (3) ρ (Λ EW ) = π/4 and δ (Λ) = σ (Λ) = 0. The corresponding numbers given in the parentheses are the relative errors compared to the exact results obtained by numerically solving the one-loop RGEs.
tan β = 10 tan β = 30  Table 6: The values of ∆δ (Λ EW ) and J (Λ EW ) obtained by means of the analytical expressions derived in section 2 for Dirac neutrinos in the MSSM with tan β = 10 or 30. The initial inputs at Λ = 10 14 GeV are chosen to guarantee that m 1 = 0.03 eV or m 3 = 0.015 eV and the best-fit values given in Eqs. (34) and (35) can be achieved at Λ EW 100 GeV. The corresponding numbers given in the parentheses are the relative errors compared to the exact results obtained by numerically solving the one-loop RGEs. Different from the differential form of RGEs for lepton flavor mixing parameters, our analytical results are of the integral form, which consist of two energy-dependent quantities, namely I β (for β = κ or ν) and ∆ τ , and the initial or final values of neutrino masses, flavor mixing angles and CP-violating phases at Λ or µ. Thus given the values of I β (for β = κ or ν) and ∆ τ against the energy scale and the initial or final values of relevant flavor parameters, one can easily calculate the RGE effects on these parameters. Moreover, compared with some previous works [49][50][51][52][53][54][55][56], we have acquired the most general and complete results for all the lepton flavor mixing parameters including the Jarlskog invariant in the standard parametrization of U for both Majorana and Dirac neutrinos. Therefore, most results of the previous works can be easily reproduced by use of ours derived in this work with some specific assumptions.
We have also carried out the numerical analysis in the MSSM for Majorana and Dirac neutrinos. Both the approximate results calculated by using the analytical formulas and the exact results obtained by numerically solving the one-loop RGEs have been shown. One can see that our approximate results can well describe the evolution behaviors of neutrino masses, flavor mixing angles, CP-violating phases and the Jarlskog invariant in most cases, and thus make them easy to be analytically understood. But in the case of Majorana neutrinos with the inverted neutrino mass ordering and a sizeable tan β, our approximate results may deviate from the exact ones to some extent. In particular, with our analytical expressions, the differences between evolution behaviors of lepton flavor mixing parameters for Dirac neutrinos and those for Majorana neutrinos with initially vanishing Majorana CP-violating phases are shown transparently, besides the entanglements among three CP-violating phases for Majorana neutrinos during the RGE running.
It is worth remarking that the integral solutions to the one-loop RGEs not only have great advantages of transparently describing the evolution behaviors of lepton flavor mixing parameters and the Jarlskog invariant, but also can be used to describe explicit radiative corrections to some fantastic flavor mixing patterns or mass textures with flavor symmetries. Of course, it may also be possible to establish an explicit connection between the phenomena of CP violation at low and high energy scales by means of such integral solutions. Thus our results are expected to be useful for model building at a superhigh energy scale so as to understand the true origin of neutrino masses and CP violation at the electroweak scale.