# Structure and function in flow networks

Nicolás Rubido  Institute for Complex Systems and Mathematical Biology,
University of Aberdeen, King’s College, AB24 3UE Aberdeen, UK  Instituto de
Física, Facultad de Ciencias, Universidad de la República, Iguá 4225,
Montevideo, 11200, Uruguay  Celso Grebogi  Institute for Complex Systems and
Mathematical Biology, University of Aberdeen, King’s College, AB24 3UE
Aberdeen, UK  Murilo S. Baptista  Institute for Complex Systems and
Mathematical Biology, University of Aberdeen, King’s College, AB24 3UE
Aberdeen, UK

###### Abstract

This Letter presents a unified approach for the fundamental relationship
between structure and function in flow networks by solving analytically the
voltages in a resistor network, transforming the network structure to an
effective all-to-all topology, and then measuring the resultant flows.
Moreover, it defines a way to study the structural resilience of the graph and
to detect possible communities.

Resistor Networks, Maximum Flows, Generalized Inverse Laplacian Matrix,
Communities.

###### pacs:

41.20-q, 95.75.Pq, 89.40.-a, 89.75.Hc, 89.75.Fb

The relation between structure and function is one of the most studied topics
in the theory of complex networks Girvan ; Newman1 ; Newman2 ; Motter ; Lai1 ;
Lai2 ; Havlin . The structure is a topological representation of the
interacting elements forming a network and it is rigorously described by the
theory of graphs. Behaviour is a functional observable of the network and can
be measured by a variety of different approaches and methodologies, e.g., it
can be quantified as a statistical description of the weights of connections
of the structure of the network Girvan ; Newman1 ; Newman2 ; Motter  or depend
on the dynamics of interactions among the elements.

This work deals with flow networks that satisfy conservation laws (Kirchhoff’s
laws). The model for the flow network is stated in terms of resistor networks
(weighted symmetric graphs), which have a source node $s$ and a sink node $t$
feeding the system, and a linear relationship between voltages and currents
VF_problem . The problem models a wide variety of physical phenomena, such as
river flow channels, transport networks, ion channels, and electrical
circuits. Its solutions establish a clear relationship between the topological
structure of the networks (namely, adjacency matrix and edge weights, assumed
known) and the functional flows passing through nodes and edges (that are a
consequence of solving the flow model). The foundation of network flow theory
roots back to Kirchhoff Kirchhoff , answering what are the current flows in an
electrical circuit when a set of voltages is applied. The solution is then
achieved by solving Kirchhoff’s equations Kirchhoff ; Ahuja ; Bollobas .

In this Letter, we present a unified approach for the fundamental relationship
between structure and function in flow networks by calculating voltages
(loads) when input currents (flows) are given. The main results of this
approach are the following.

First, a novel formula for the voltages in the flow network model, i.e., the
_Voltage-Flow (VF) problem_ [Eq. (2)], is analytically derived. The solution
[Eq. (6)] is expressed in terms of the known graph’s weighted Laplacian matrix
eigenspace base. In particular, we show that the behaviour of the flows is not
only dependent on the matrix spectra, but intrinsically connected to the
eigenvectors of this matrix.

Second, we show that the _equivalent resistance_ (also known as two-point
distance resistance) formula Wu-Arpita  is retrieved [Eq. (7)] without further
assumptions from the analytical expression for the voltages induced by the
single $s$-$t$ pair. The equivalent resistance between any two nodes is a
topological measure that gives further information about the structure of any
graph. It is shown here how it constitutes a highly useful tool to reduce the
complexity of any network structure: an arbitrary topology is transformed to
an effective all-to-all complete graph.

Third, using both the voltage differences and the equivalent resistance, an
_effective functional flow adjacency_ (EFFA) between any two points in the
circuit is constructed [Eq. (8)]. These effective flows are the real
measurable observables that are obtained when a probe is placed between any
two nodes in the network. The EFFA matrix expresses the level of connectivity
in terms of flow transmission; the functional relationship among the network
components in the effective all-to-all topology. In particular, if two nodes
are not directly connected, an effective current between them may still exist.

Fourth, as an application of our approach, we show how to identify the nodes
that are mainly responsible for the transmission of flow, how to detect
qualitatively communities of the network, and how to infer maximal node
outflow for different graph topologies and resistor distributions. This
reveals topological issues that are related to node maximum capacities in a
natural manner.

The starting point for our unified approach is the calculation of voltages in
a network by solving the VF problem. The VF problem is set by interpreting the
graph structure as a linear resistor network circuit, which then receives an
input flow $I$ at node $s$ and leaves at node $t$ (a current $I$ enters the
circuit at a source and leaves at a sink). The network structure is given by a
connected strict graph $\mathcal{G}=\\{\mathcal{V},\,\mathcal{E}\\}$ (where
$\mathcal{V}$ and $\mathcal{E}$ are the sets of nodes and edges, respectively)
with symmetric edge weights $W_{ij}\equiv A_{ij}/R_{ij}$ for
$i,\,j=1,\ldots,\,N$ ($N$ being the number of nodes in $\mathcal{V}$, $A_{ij}$
the $ij$-th element of the adjacency matrix, and $R_{ij}$ the edge’s
resistance).

Imposing conservation of charge in every node of the network, i.e., a graph
where at each node the current arriving at it equals the amount leaving it
(first Kirchhoff’s law), the flow vector is
$(\vec{F}^{(st)})_{i}=I\,\left(\delta_{is}-\delta_{it}\right)$ for
$i=1,\,\ldots,\,N$, with $\delta_{ij}$ being the Kronecker delta. Thus, using
the laws of Kirchhoff and Ohm, the net current in any node is given by

$F_{i}^{(st)}=I\,\left(\delta_{is}-\delta_{it}\right)=\sum_{j=1}^{N}\frac{A_{ij}}{R_{ij}}\left(V_{i}^{(st)}-V_{j}^{(st)}\right)\,,$
(1)

where $V_{i}^{(st)}$ is the voltage potential at node $i$ given the particular
$s$-$t$ pair. Equation (1) is rearranged such that the VF problem is expressed
in matrix form

$\mathbf{G}\,\vec{V}^{(st)}=\vec{F}^{(st)}\,,$ (2)

where the upper-indexes indicate that the problem depends on the location of
the $s$-$t$ nodes on $\mathcal{V}$ and $\mathbf{G}$ is the weighted Laplacian
matrix. The entries of $\mathbf{G}$ are

$G_{ij}=\left\\{\begin{array}[]{lcl}\sum_{k=1}^{N}W_{ik}&\text{if}&i=j\,,\\\
-W_{ij}&\text{if}&i\neq j\,.\end{array}\right.$ (3)

Then, the VF solution is achieved once the voltages in each node are found
from inverting $\mathbf{G}$. However, the rank of $\mathbf{G}$ is $N-1$
($\det\left(\mathbf{G}\right)=\prod_{n=0}^{N-1}\lambda_{n}=0$) and direct
inversion is not possible Bollobas ; FanChung .

The Laplacian inversion is addressed here by a different interpretation of the
derivation of the generalized inverse matrix, also known as Moore-Penrose
matrix Bollobas ; Wu-Arpita .

The eigenvalue problem for the Laplacian is given by
$\mathbf{G}\,\mathbf{P}=\mathbf{P}\,\mathbf{\Lambda}$, where
$\Lambda_{ij}=\delta_{ij}\,\lambda_{n}$ ($n=i-1=0,\ldots,N-1$) is a diagonal
matrix containing all eigenvalues and
$\mathbf{P}=\\{\vec{v}_{0},\,\vec{v}_{1},\,\ldots,\,\vec{v}_{N-1}\\}$ is the
matrix whose columns are the eigenvectors of $\mathbf{G}$. We prove here that
to express any element of $\mathbf{G}$, the only elements needed from the
spectral decomposition are the non-null eigenvalues and corresponding
eigenvectors. This results in the following decomposition of $\mathbf{G}$

$\mathbf{G}=\mathbf{P}_{r}\,\mathbf{\Lambda}_{r}\,\mathbf{P}_{r}^{T}\,,$ (4)

where $T$ indicates transpose, $\mathbf{P}_{r}\in\mathbb{R}^{N\times(N-r)}$,
$\mathbf{\Lambda}_{r}\in\mathbb{R}^{(N-r)\times(N-r)}$, and $r$ being the
number of connected components of $\mathcal{G}$, which in this Letter is fixed
at $r=1$.

Equation (4) is derived by analyzing each Laplacian matrix entry using the
full spectral decomposition. In particular,
$G_{ij}=\sum_{k=1}^{N}P_{ik}\lambda_{k-1}P^{-1}_{kj}$, where for $k=1$,
$\lambda_{0}=0$, so no contribution is added to the sum. Thus, observing that
the inverse of an orthogonal matrix is its transpose
($P^{-1}_{kj}=P_{jk}=\left(\vec{v}_{k-1}\right)_{j}$), the $ij$’s entry of
$\mathbf{G}$ is given by:
$G_{ij}=\sum_{k=2}^{N}\left(\vec{v}_{k-1}\right)_{i}\,\lambda_{k-1}\,\left(\vec{v}_{k-1}\right)_{j}$,
which is easily extendible to various connected components making the lower
bound for the summation index $k=r+1$.

However, removing the null eigenspace of the spanning set of eigenvectors, the
rank of the new $\mathbf{G}$ in Eq. (4) is $N$ and inversion can be directly
performed. Consequently, each element of the generalized inverse Laplacian
matrix $Z_{ij}\equiv(\mathbf{G}^{-1})_{ij}$ is found from

$Z_{ij}=\sum_{n=1}^{N-1}\left(\vec{v}_{n}\right)_{i}\frac{1}{\lambda_{n}}\left(\vec{v}_{n}\right)_{j}\,,$
(5)

with the new index $n=k-1$ and $i,\,j=1,\ldots,N$. It is easy to show that
$\sum_{k=1}^{N}G_{ik}\,Z_{kj}=\delta_{ij}-1/N$, where the extra constant
factor comes from the incompleteness of the restricted eigenspace and it
cancels out when calculating voltage differences.

Returning to Eq. (2), the _solution for the VF problem_ at a given node
$i\in\mathcal{V}$ is
$\left(\vec{V}^{(st)}\right)_{i}=\sum_{j=1}^{N}Z_{ij}\,\left(\vec{F}^{(st)}\right)_{j}=I\,\left(Z_{is}-Z_{it}\right)$,

$\left(\vec{V}^{(st)}\right)_{i}=I\,\sum_{n=1}^{N-1}\frac{\left(\vec{v}_{n}\right)_{i}}{\lambda_{n}}\left[\left(\vec{v}_{n}\right)_{s}-\left(\vec{v}_{n}\right)_{t}\right]\,.$
(6)

This formula is the backbone of our unified approach. It is applicable to any
connected weighted graph and is extendible to many sources and sinks with
different inputs and outputs as long as flow conservation is fulfilled
($\sum_{i}F_{i}^{(st)}=0$ and $F_{i}=0\;\forall\,i\neq s$ or $t$).
Furthermore, it allows to derive the equation for the equivalent resistance
between any two nodes, to construct the effective flows, and later apply them
to, e.g., community detection.

The derivation of the _equivalent resistance_ $\rho_{ij}$ between any two
nodes Wu-Arpita  is done as follows. Set the source at a starting node ($i=s$)
and the sink at an ending point ($j=t$), then Eq. (6) provides an exact closed
formula for $\rho_{ij}$ in terms of eigenvalues and eigenvectors of the
weighted Laplacian matrix. The reason for doing this is that the voltage
difference between these two nodes gives the incoming current times the
resistance between them, hence,
$(\vec{V}^{(st)})_{s}-(\vec{V}^{(st)})_{t}=I\,\rho_{st}$. Therefore, using the
same procedure for every pair of nodes in the graph all $\rho_{ij}$’s are
determined,

$\rho_{ij}=\sum_{n=1}^{N-1}\frac{1}{\lambda_{n}}\left[\left(\vec{v}_{n}\right)_{i}-\left(\vec{v}_{n}\right)_{j}\right]^{2}\,,$
(7)

providing a solution that is independent of where the $s$-$t$ pair is placed.

The matrix element $\rho_{ij}$ is a topological measure that weighs all paths
between any pair of nodes $i$ and $j$. All paths between these nodes that are
allowed by the adjacency structure of the graph are collapsed to one
equivalent link with weight $\rho_{ij}$. The result is that a single effective
link is created between $i$ and $j$. The important contribution is that _this
quantity presents the possibility to diminish the complexity of any graph to a
simpler all-to-all equivalent topology_ with a betweenness score given by
$\rho_{ij}$.

As a working example, results of applying the solution given by Eqs. (6) and
(7) to a ring graph, namely $C_{N}$, of $N=16$ nodes with resistors equal to
one ($R_{ij}=1\;\forall\,i,j$) and input current $I=1$ are shown in Fig. 1.
For this network, the set of unordered eigenvalues is given by
$\lambda_{n}=2-2\cos(2\pi\,n/N)$, with $n=0,\ldots,N-1$ FanChung .

The top left panel in Fig. 1 shows that for the $N\,(N-1)=240$ possibilities
of $s$-$t$ pairs in $C_{N}$ the maximum voltages achieved in the ring do not
surpass $2$ (arbitrary units). This is the maximum difference that the
eigenvector coordinates can achieve in the non-normalized frame (top right
panel in Fig. 1), hence giving a maximum voltage difference of $4$. This
result shows that not only the eigenvalues of the Laplacian are important to
infer the network’s functional behaviour, but eigenvectors are also needed.
Also, the maximum (minimum) is periodic, in the sense that the largest
difference happens every time the source and sink are separated by $8$ edges.

Figure 1: The top left panel has the maximum (black) and minimum (red)
voltages that are achieved in an unweighed ($R_{ij}=1\;\forall\,i,j$) ring
graph ($C_{N}$) with $N=16$ nodes for every $s$-$t$ pair in units of $I=1$ for
the input current at node $s$. The top right panel shows the eigenvectors
($\sqrt{N}\vec{v}_{n}$) coordinate values (color scale) of the corresponding
weighed Laplacian matrix. The bottom left panel exhibits in color coded, for
all these possible input-output configurations, the voltages $V_{i}^{(st)}$ of
every node. The bottom right panel shows the equivalent resistor matrix
($\rho_{ij}$) for every pair of nodes in $C_{N}$.

The bottom panels in Fig. 1 show the voltage solutions for all the $s$-$t$
pairs (left) and all equivalent resistances (right). The resulting values show
the functional and structural symmetry that this simple network has. As all
edge weights are unity ($R_{ij}=1$), the equivalent resistance $\rho_{ij}$
maximum can only be $4$. This is easily seen when calculating it directly by
means of series and parallel resistor calculations. For instance, for nodes
$i=1$ and $j=9$ there are two parallel paths with $8$ edges in series, hence:
$\rho_{1,9}^{-1}=\frac{1}{8}+\frac{1}{8}\;\Rightarrow\;\rho_{1,9}=4$, which is
the same value as element $\rho_{1,9}$ in the bottom right panel of Fig. 1
(found from Eq. (7)).

Having the voltage difference and equivalent resistance given by Eqs. (6) and
(7), the third result of our unified approach is addressed: the construction
of the _effective edge flows_ , namely, $\phi_{ij}^{(st)}$. We introduce these
flows by finding the voltage differences among any two nodes in the effective
topology given by $\rho_{ij}$. Thus, the EFFA matrix elements for each $s$-$t$
configuration are given by

$\phi_{ij}^{(st)}=\frac{1}{\rho_{ij}}\left[\left(\vec{V}^{(st)}\right)_{i}-\left(\vec{V}^{(st)}\right)_{j}\right]\,,$
(8)

where $(\vec{V}^{(st)})_{i}>(\vec{V}^{(st)})_{j}$ such that an effective
current is flowing from node $i$ to $j$, otherwise $\phi_{ij}^{(st)}=0$.
Hence, the $s$-$t$ dependent EFFA matrix defines an effective directed
network, which breaks the symmetry of the effective all-to-all topology
structure given by $\rho_{ij}$. In particular, the physical current on each
edge $ij\in\mathcal{E}$ of the graph is given by
$I_{ij}^{(st)}=\frac{A_{ij}}{R_{ij}}\,[(\vec{V}^{(st)})_{i}-(\vec{V}^{(st)})_{j}]$.

The voltage differences define a functional relationship among nodes and the
equivalent resistance an effective topology structure. The introduction of
_the EFFA relates both structure and functional characteristics of the graph_
, thus, it enables to find functional hubs that might not be there in the
topological structure, and detect highly connected areas (communities) in a
functional way. In this sense, the $s$-$t$ analysis acts similarly as how a
dye does when introduced to a circulatory system of a person: it highlights a
certain part of the network structure to detect what is being sought, e.g.,
clog blocking, arteries stiffness, etc.

Figure 2: The left panel shows the effective flow functional adjacency (EFFA)
matrix for the unweighed $C_{N}$ graph of Fig. 1, where all source-sink pairs
have been averaged, i.e., $\phi_{ij}$. The right panel exhibits the same EFFA
matrix for a particular source-sink pair, namely, $s=2$ and $t=10$
($\phi_{ij}^{(2,10)}$). Each matrix element’s magnitude is associated with a
color scale.

The mean EFFA matrix is defined to detect statistically relevant functional
observables. It is calculated by all the effective currents (Eq. (8)) among
the different pairs averaged over all possible realizations of source-sink
pairs. Its matrix elements are

$\phi_{ij}\equiv\left\langle\phi_{ij}^{(st)}\right\rangle_{(st)}=\sum_{s=1}^{N}\sum_{t=1\,(t\neq
s)}^{N}\frac{\phi_{ij}^{(st)}}{N\,(N-1)}\,,$ (9)

which define a matrix that constitutes an effective structural flow quantity
and is related to the equivalent resistance matrix. The EFFA matrix and its
averaged counterpart for the example of the ring graph are shown in Fig. 2,
where the left panel shows the mean EFFA matrix and the right panel exhibits
the symmetry breaking for a particular $s$-$t$ pair.

With the three quantities analytically expressed (voltages in Eq. (6),
equivalent resistances in Eq. (7), and the effective currents in Eq. (8)), we
investigate graphs  random graphs (RN) Erdos , small-world (SW) networks
Strogatz , scale-free (SF) topologies Barabasi , and numerically generated
clustered networks (CN) of $N=2^{9}$ nodes. The following conclusions
constitute the application of our unified approach to these concrete
topologies. In particular, the differences between the respective adjacency
and effective topologies (equivalent resistance) are shown in the left and
centre columns of Figs. 3 and 4.

To study other features of the EFFA matrix that are not represented in its
averaged version and to answer how the flow is effectively distributed in any
given topology, the probability density function (PDF) of the effective edge
flows $\phi_{ij}^{(st)}$ for every particular $s$-$t$ is calculated. This does
not only tells how diffusive or localized the transport of currents is as a
function of the effective topology, but also gives the maximum flow values a
particular graph can have, i.e., the maximum edge capacity.

Results show that the EFFA PDFs for RN, SW, and SF networks are roughly
invariant under changes of the $s$-$t$ pair. However, SW and SF networks
produce more high flow magnitudes ($\phi_{ij}^{(st)}\simeq 1$) than RN for a
broad range of parameters graphs , hence, RNs distribute flows in a more
uniform way than SF or SW networks. These observations are seen in the example
shown on the right column of Fig. 3.

Figure 3: The figure shows the following unweighted graphs: a random (top
row), a scale-free (centre row), and a small-world (bottom row) network of
$N=2^{9}$ nodes and connecting probability $p=0.05$ graphs . The adjacency
matrix (left column), equivalent resistance matrix (centre column in color
coded), and effective functional flow adjacency PDFs (right column) for all
$s$-$t$ configurations (gray scale) are shown for each graph.

The detection of communities is done by analysis of the mean EFFA matrix and
the equivalent resistance. In general, _the application of the EFFA analysis
on numerically generated CNs qualitatively exhibits the communities that the
networks have_. This is seen in the centre and right panels of the example
exhibited in Fig. 4, where the communities are spotted as dark areas, both in
$\rho_{ij}$ and $\phi_{ij}$, though the contrast in the mean EFFA matrix is
higher than in the equivalent resistance, implying that communities are better
identified by $\phi_{ij}$. _The $\phi_{ij}$ matrix_ not only detects the
communities but also _shows that the effective flows between communities are
larger than the ones within each community_. Hence, the intra-edges connecting
communities must have a high flow capacity as they represent the vulnerable
links in the network. Removal of any of these edges generate drastic changes
in the flows, leading to the isolation of communities.

Figure 4: The figure shows the adjacency matrix (left column), equivalent
resistance (middle column), and the average effective functional adjacency
matrix (right column) for an unweighed clustered network formed by $4$ random
graphs of $N_{c}=2^{7}$ graphs . The network has a total of $N=2^{9}$ nodes
and each matrix element magnitude is associated with a color scale.

The observations on flow transmission and community detection have been
previously shown by doing structural analysis of networks, e.g., in Refs.
Girvan ; Newman1 ; Newman2 ; Motter ; Lai1 ; Lai2 ; Havlin . In particular,
the possibility of using voltages as a betweenness measure score is pointed
out in Refs. Girvan ; Newman1 ; Newman2  and is related to shortest-paths and
random walks. With our unified approach, and the four results that are
addressed in this Letter, we show the practicality of it. Moreover, the linear
relationship between loads and flows, and conservation of charge, are
restrictions needed to define the flow network on top of the structure but not
for the structure itself. Thus, the flow analysis is applicable to any network
structure outside the domain of flow networks as long as the graph to be
analysed is connected.

An estimate of the computational time that our approach requires is as
follows. To calculate the voltages one needs the eigenvectors and eigenvalues
of the Laplacian matrix. This is calculated in time $\mathcal{O}(N^{\alpha})$,
with $\alpha$ being an exponent that nowadays rounds $2$ and $N$ the number of
nodes in the network. Then, if all possible source-sink pairs are analysed,
the flows require $\mathcal{O}(N^{3})$ to be found. This means that the
algorithm is not as effective as other methods Andrea  but its mathematical
formulas still provide great information for obtaining upper and lower bounds
for currents without the need for simulations. Furthermore, community
detection can be quantitatively described if other techniques are used, such
as PDF analysis of the EFFA matrix elements.

Authors acknowledge the Scottish University Physics Alliance (SUPA).

## References

  * (1) M. Girvan and M. E. J. Newman, PNAS 99(12), 7821-7826 (2002).
  * (2) M. E. J. Newman and M. Girvan, Phys. Rev. E 69, 026113 (2004).
  * (3) M. E. J. Newman, Eur. Phys. J. B 38, 321-330 (2004).
  * (4) A. E. Motter, Phys. Rev. Lett. 93(9), 098701 (2004).
  * (5) Rui Yang, Wen-Xu Wang, Ying-Cheng Lai, and Guanrong Chen, Phys. Rev. E 79, 026112 (2009).
  * (6) Wen-Xu Wang and Ying-Cheng Lai, Phys. Rev. E 80, 036109 (2009).
  * (7) Roni Parshani, Sergey V. Buldyrev, and Shlomo Havlin, Phys. Rev. Lett. 105, 048701 (2010).
  * (8) In particular, the approach that we take is derived from the maximum flow problem, which corresponds to finding a feasible flow through a single-source $s$ and a single-sink $t$ that is maximum when costs and capacities at the edges are given. The solution for this problem is equal to the minimum capacity of an $s$-$t$ cut in the network as stated in the max-flow min-cut theorem Ahuja ; Bollobas . On the other hand, if the $s$-$t$ flow is fixed, then the problem turns into finding how is the transmission given and what is its implications to the structure and function of the network.
  * (9) G. Kirchhoff, Vorlesungen über Mechanik (Wien, Wilhelm, Ed. 1864-1928).
  * (10) Ravindra K. Ahuja, Thomas L. Magnanti, and James B. Orlin, Network Flows: Theory, Algorithms, and Applications (Prentice-Hall, 1993).
  * (11) B. Bollobás, Modern Graph Theory (Springer-Verlag, New York, 1998, chap. 9).
  * (12) Fan R. K. Chung, Spectral Graph Theory (Am. Math. Soc. and CBMS 92, 1997).
  * (13) F. Y. Wu, J. Phys. A: Math. Gen. 37, 6653-6673 (2004); Arpita Ghosh, Stephen Boyd, and Amin Saberi, SIAM Rev. 50(1), 37-66 (2008); and references therein.
  * (14) Starting with a ring graph of $N=2^{9}$ nodes, a tuning parameter $p$ controls the addition of links. For the RN, SW, and SF graphs the process is straightforward (taken from Refs. Erdos ; Strogatz ; Barabasi ) and $p$ is related to the probability of adding an edge. The RN process follows the strategy of Erdos . The SF graph is constructed as in the rewiring procedure of Ref. Barabasi . The SW graph is constructed as in Ref. Strogatz . For the CN, four RN of $N_{c}=2^{7}$ are created with equal statistics (probability $p=0.01$), then approximately $10$ links are added between them randomly. In addition, weights are included from the uniform random distribution such that $R_{ij}\in(0,1]\;\forall\,i,j$.
  * (15) P. Erdos and A. Renyi, Publ. Math. 6, 290 (1959).
  * (16) D. J. Watts and S. H. Strogatz, Nature 393, 440-442 (1998).
  * (17) A.-L. Barabási and R. Albert, Science 286, 509 (1999).
  * (18) Andrea Lancichinetti, Santo Fortunato, and Filippo Radicchi, Phys. Rev. E 78, 046110 (2008).


# Implication of cosmological upper bound on the validity of golden ratio
neutrino mixings under radiative corrections

Y Monitar Singh<EMAIL_ADDRESS>Department of Physics,
Manipur University, Imphal-795003, India  M Shubhakanta Singh
<EMAIL_ADDRESS>Department of Physics, Manipur University,
Imphal-795003, India  N. Nimai Singh<EMAIL_ADDRESS>Department of Physics,
Manipur University, Imphal-795003, India  Research Institute of Science and
Technology, Imphal-795003, India

###### Abstract

We study the implication of the most recent cosmological upper bound on the
sum of three neutrino masses, on the validity of the golden ratio (GR)
neutrino mixings defined at high energy seesaw scale, considering the
possibility for generating low energy values of neutrino oscillation
parameters through radiative corrections in the minimal supersymmetric
standard model (MSSM). The present study is consistent with the most stringent
and latest Planck data on cosmological upper bound, $\sum|m_{i}|<0.12$ eV. For
the radiative generation of sin$\theta_{13}$ from an exact form of golden
ratio (GR) neutrino mixing matrix defined at high seesaw energy scale, we take
opposite CP parity mass eigenvalues ($m_{1},-m_{2},m_{3}$) with a non-zero
real value of $m_{3}$, and a larger value of $\tan\beta>60$ in order to
include large effects of radiative corrections in the calculation. The present
analysis including the CP violating Dirac phase and SUSY threshold
corrections, shows the validity of golden ratio neutrino mixings defined at
high seesaw energy scale in the normal hierarchical (NH) model. The numerical
analysis with the variations of four parameters viz. $M_{R}$, $m_{s}$,
$\tan\beta$ and $\bar{\eta_{b}}$, shows that the best result for the validity
is obtained at $M_{R}=10^{15}$ GeV, $m_{s}=1$ TeV, $\tan\beta=68$ and
$\bar{\eta_{b}}=0.01$. However, the analysis based on inverted hierarchical
(IH) model does not conform with this latest Planck data on cosmological bound
but it still conforms with earlier Planck cosmological upper bound
$\sum|m_{i}|<0.23$ eV, thus indicating possible preference of NH over IH
models.

Keywords : Radiative corrections, supersymmetric standard model,
renormalisation group equations, mixing matrix, golden ratio mixing.

##  I Introduction

The values of neutrino oscillation parameters have been continuously updated
with the advancement in the technology of neutrino oscillation experiments
esteban2020fate ; de20212020 ; mc  and these updated experimental data are
also required to compare with the theoretically predicted values. The latest
Planck data on the cosmological upper bound on the sum of the three absolute
mass eigenvalues given by $\sum|m_{i}|<0.12$ eV aghanim2020planck , may be
seriously considered while comparing with other neutrino oscillation
parameters although there are also a lot of constraints associated with such
cosmological probetanseri2022updated . The theoretical predictions of these
neutrino oscillation parameters are in general defined at very high energy
seesaw scale, and the experimental data on the other hand are defined at low
energy scale of the order of $10^{2}$ GeV. In order to make a bridge between
these two energy scales, we need a set of renormalisation group equations
(RGEs) for quantum radiative corrections das2005numerical ;
king2000renormalisation . We can use two different approaches for running the
RGEs from high-energy scale to low-energy scale. In the first approach, the
running of RGEs is carried out through the neutrino mass matrix $m_{LL}$ as a
whole, and at every energy scale one can extract neutrino masses and mixing
angles through the diagonalisation of the neutrino mass matrix calculated at
that particular energy scale chankowski1993renormalization ;
babu1993renormalization ; as ; parida1998low ; king2001inverted . In the
second approach, the running of RGEs can be carried out directly in terms of
neutrino mass eigenvalues and three mixing angles with phases
chankowski2000fixed ; brignole2003low ; antusch2003running . In both cases,
the RGEs of all the neutrino parameters and the RGEs of various coupling
constants are solved simultaneously and both approaches give almost consistent
results das2005numerical . For the present analysis, we shall use second
approach which is more convenient to handle in the numerical analysis of RGEs
of neutrino oscillation parameters.

Various discrete symmetry groups like $S_{4}$, $A_{4}$, $A_{5}$ etc. which are
defined at very high energy scale, can lead to various leptonic mixing
matrices such as bi-maximal (BM), tri-bimaximal (TBM) and Golden ratio
(GR)Kajiyama_2007 . All these specific leptonic mixing matrices have their own
respective leptonic mixing angles, and two of the mixing angles ($\theta_{23}$
and $\theta_{12}$) are in good agreement with the respective non-zero neutrino
mixing angles at low energy scale. In all the above three leptonic mixing
matrices, the three leptonic neutrino mixing angles are defined at very high
energy scale, with reactor neutrino mixing angle ($\theta_{13}$) equals to
zero. The radiative magnification of reactor neutrino mixing angle
($\theta_{13}$) is studied with various leptonic mixing matrices such as BM,
TBM and GR zhang2016viability ; pramanick2021radiative .

GR neutrino mixing pattern has certain advantages over the other two neutrino
mixing patterns BM and TBM in the evolution of mixing angles under radiative
corrections as solar mixing angle ($\theta_{12}$) is always found to increase
with deccrease in energy scale. The generation of right order non-zero value
of reactor neutrino mixing angle ($\theta_{13}$) at low energy scale,
consistent with the latest cosmological upper bound on the sum of three
absolute neutrino mass eigenvalues, $\sum|m_{i}|<0.12$ eV, is mainly addressed
in the present study. Two cases of neutrino mass hierarchical models namely
normal hierarchy (NH) and inverted hierarchy (IH) are considered when we take
input mass eigenvalues at very high energy seesaw scale.

A brief description of an exact form of golden ratio mixing matrix ($U_{GR}$)
is given by king2017unified ,

$U_{GR}=\left(\begin{array}[]{ccc}\frac{\phi}{\sqrt{2+\phi}}&\frac{1}{\sqrt{2+\phi}}&0\\\
-\frac{1}{\sqrt{4+2\phi}}&\frac{\phi}{\sqrt{4+2\phi}}&\frac{1}{\sqrt{2}}\\\
\frac{1}{\sqrt{4+2\phi}}&-\frac{\phi}{\sqrt{4+2\phi}}&\frac{1}{\sqrt{2}}\end{array}\right)$
(1)

where $\phi$ has following properties :

$\phi=\phi^{2}-1=1+\frac{1}{\phi}=\frac{1+\sqrt{5}}{2}\approx 1.62$

and

$\frac{1}{\phi^{2}}=\frac{1-\frac{1}{\sqrt{5}}}{1+\frac{1}{\sqrt{5}}}\approx
0.382.$

It also predicts $\sin\theta_{13}=0$, $\sin\theta_{23}=\frac{1}{\sqrt{2}}$ and
$\tan\theta_{12}=\frac{1}{\phi}$, leading to :

$\theta_{12}=\tan^{-1}\left(\frac{1}{\phi}\right)=31.72^{0}$

and golden ratio is sometimes enforced by $A_{5}$ everett2009icosahedral . The
$U_{GR}$ is a special case for the $\mu-\tau$ symmetric mass matrix,

$M_{\nu}=\left(\begin{array}[]{ccc}D&A&\pm A\\\ A&B&C\\\ \pm
A&C&B\end{array}\right)$

where

$\tan 2\theta_{12}=\frac{2\sqrt{2}A}{B\pm C-D}.$

For case $B\pm C-D=A$, the mixing matrix goes to the tri-bimaximal mixing
matrix ($U_{TBM}$) with $\tan 2\theta_{12}=2\sqrt{2}$, and for case $B\pm
C-D=\sqrt{2}A$, the mixing matrix goes to $U_{GR}$ with $\tan 2\theta_{12}=2$
fukuyama2017twenty ; koide2004universal . When D=0, the structure of the mass
matrix predicts $\tan^{2}\theta_{12}=\frac{m_{1}}{m_{2}}\approx 0.382$, and
$m_{1}$ and $m_{2}$ are two neutrino mass eigenvalues. To check the validity
of GR neutrino mixings at high energy scale, we consider a large value of
tan$\beta>60$ in order to include large effects of radiative corrections in
the calculation of neutrino masses, mixing angles and to satisfy the latest
cosmological upper bound $\sum|m_{i}|<0.12$ eV lorenz2021reconstruction ;
choudhury2020updated  in both normal and inverted hierarchical mass models.

The paper is organised as follows. In section 2, we briefly outline the main
points on renormalisation group analysis for neutrino oscillation parameters
with phases. In section 3, we present numerical analysis of RGEs for GR
neutrino mixing matrix. In section 4, we give results and discussion. In
section 5, we give summary and conclusion.

##  II Renormalisation group analysis for neutrino oscillation parameters with
phases

We briefly present the main formalism for the evolution of neutrino
oscillation parameters ma2002all ; antusch2005charged ; antusch2002neutrino
from high energy seesaw scale to low energy scale through the RGEs with CP
violating phase in the MSSM, including appropriate SUSY threshold corrections.
The neutrino masses can be described by lowest-dimensional neutrino mass
operator compatible with the gauge symmetries of the SM. This operator reads
in the SM antusch2003running

$\mathcal{L}_{K}^{SM}=-\frac{1}{4}K_{gf}\bar{l}_{Lc}^{{C}^{g}}\varepsilon^{cd}\phi_{d}l_{Lb}^{f}\varepsilon^{ba}\phi_{a}+h.c.,$
(2)

and in its minimal supersymmetric extension, the MSSM

$\mathcal{L}_{K}^{MSSM}=-\frac{1}{4}K_{gf}\mathbb{l}_{c}^{g}\varepsilon^{cd}\mathbb{h}_{d}^{(2)}\mathbb{l}_{b}^{f}\varepsilon^{ba}\mathbb{h}_{a}^{(2)}|_{\theta\theta}+h.c.,$
(3)

where $l_{L}^{C}$ is the charge conjugate of a lepton doublet. $\varepsilon$
is the totally antisymmetric tensor in 2 dimensions, and a, b, c, d $\in$
{1,2} are $SU(2)_{L}$ indices. The double-stroke letters $\mathbb{l}$ and
$\mathbb{h}$ denote lepton doublets and the up-type Higgs superfield in the
MSSM. The coefficients $K_{gf}$ are of mass dimension -1 and related to the
Majorana neutrino mass matrix as $M_{\nu}=K\langle H\rangle^{2}$, where
$\langle H\rangle=174\text{GeV}$ is the vacuum expectation value of Higgs
field ($v_{0}$).

The most plausible explanation for neutrino mass is given by see-saw mechanism
davidson1987universal . The neutrino mass matrix $m_{LL}$(t) which is
generally obtained from see-saw mechanism, is expressible in terms of K(t),
the coefficient of the dimension five neutrino mass operator in the scale-
dependent manner, $t=\ln(\mu/1\text{GeV})$,

$m_{LL}(t)=v_{u}^{2}(t)K(t)$ (4)

where the vacuum expectation value (VEV) is $v_{u}=v_{o}\sin\beta$ and
$v_{o}$=174 GeV in the minimal supersymmetric standard model (MSSM). After
diagonalization of K(t), the above eq.(4) can be written in terms of mass
eigenvalues as follows chankowski2000fixed

$m_{i}(t)=v_{u}^{2}(t)K_{i}(t);i=1,2,3.$ (5)

This expression can be simplified as

$\frac{d}{dt}m_{i}(t)=m_{i}(t)\left(\frac{1}{K_{i}(t)}\frac{d}{dt}K_{i}(t)+\frac{2}{v_{u}(t)}\frac{d}{dt}v_{u}(t)\right).$
(6)

Now, considering the phases in neutrino mixing matrix, we parameterize the
Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix as,

$\displaystyle U_{PMNS}=$
$\displaystyle\left(\begin{array}[]{ccc}c_{13}c_{12}&c_{13}s_{12}&s_{13}e^{-i\delta}\\\
-c_{23}s_{12}-c_{12}s_{13}s_{23}e^{-i\delta}&c_{12}c_{23}-s_{12}s_{13}s_{23}e^{-i\delta}&c_{13}s_{23}\\\
s_{12}s_{23}-c_{12}s_{13}c_{23}e^{-i\delta}&-c_{12}s_{23}-c_{23}s_{13}s_{12}e^{-i\delta}&c_{13}c_{23}\end{array}\right)$
(10) $\displaystyle\times$
$\displaystyle\left(\begin{array}[]{ccc}e^{i\alpha_{1}}&0&0\\\
0&e^{i\alpha_{2}}&0\\\ 0&0&1\end{array}\right)$ (14)

where $s_{ij}=\sin\theta_{ij}$, $c_{ij}=\cos\theta_{ij}$, $\delta$= Dirac
phase, $\alpha_{1}$=first Majorana phase, $\alpha_{2}$=second Majorana phase.
Here, the three mixing angles are defined as
$\tan\theta_{12}=\frac{|U_{e2}|}{|U_{e1}|}$, $\tan\theta_{23}=\frac{|U_{\mu
3}|}{|U_{\tau 3}|}$ and $\sin\theta_{13}=|U_{e3}|$.

The RGEs for $K_{i}(t)$ and $v_{u}(t)$ in the basis where charged lepton mass
matrix is diagonal, for one-loop order in MSSM, in the energy range from
$M_{R}$ to $M_{SUSY}$, are given by babu1993renormalization ; casas2000nearly
; arason1992renormalization

$\frac{1}{K_{i}(t)}\frac{d}{dt}K_{i}(t)=\frac{1}{16\pi^{2}}\sum_{f=e,\mu,\tau}\left[-\frac{6}{5}g_{1}^{2}-6g_{2}^{2}+6h_{t}^{2}+2h_{f}^{2}U_{fi}^{2}\right]$
(15)

and

$\frac{1}{v_{u}(t)}\frac{d}{dt}v_{u}(t)=\frac{1}{16\pi^{2}}\left[\frac{3}{20}g_{1}^{2}+\frac{3}{4}g_{2}^{2}-3h_{t}^{2}\right].$
(16)

The RGEs for $K_{i}$ and $v_{0}$ in the basis where the charged lepton mass
matrix is diagonal, for one-loop order in SM, in the energy range from
$M_{SUSY}$ to $M_{Z}$, are given by babu1993renormalization ; casas2000nearly
; arason1992renormalization

$\frac{1}{K_{i}(t)}\frac{d}{dt}K_{i}(t)=\frac{1}{16\pi^{2}}\sum_{f=e,\mu,\tau}\left[-3g_{2}^{2}+2\lambda+6h_{t}^{2}+h_{\tau}^{2}+h_{f}^{2}U_{fi}^{2}\right]$
(17)

and

$\frac{1}{v_{0}(t)}\frac{d}{dt}v_{0}(t)=\frac{1}{16\pi^{2}}\left[\frac{9}{20}g_{1}^{2}+\frac{9}{4}g_{2}^{2}-3h_{t}^{2}-3h_{b}^{2}-h_{\tau}^{2}\right],$
(18)

where $g_{1},g_{2}$ are gauge couplings, and $h_{t},h_{b},h_{\tau}$ and
$\lambda$ are top-quark, bottom -quark, tau -lepton Yukawa couplings and SM
quartic Higgs coupling respectively. As VEV can affect mass terms in the RGEs,
we have two possible set of RGEs of neutrino masses where one is scale-
dependent VEV and other is scale-independent VEV. The RGEs of neutrino mass
eigenvalues for both scale dependent VEV and scale independent VEV can be
written as agarwalla2007neutrino ; antusch2003running

$\frac{d}{dt}m_{i}=-2F_{\tau}\left(P_{i}+Q_{i}\right)m_{i}-F_{u}m_{i},(i=1,2,3).$
(19)

where,

$P_{1}=s_{12}^{2}s_{23}^{2},P_{2}=c_{12}^{2}s_{23}^{2},P_{3}=c_{13}^{2}c_{23}^{2},$
$Q_{1}=-\frac{1}{2}s_{13}\sin 2\theta_{12}\sin
2\theta_{23}\cos\delta+s_{13}^{2}c_{12}^{2}c_{23}^{2},$
$Q_{2}=\frac{1}{2}s_{13}\sin 2\theta_{12}\sin
2\theta_{23}\cos\delta+s_{13}^{2}s_{12}^{2}c_{23}^{2},$ $Q_{3}=0.$

For scale-dependent VEV in the case of MSSM with $\mu\geq m_{s}$,

$F_{\tau}=-\frac{h_{\tau}^{2}}{16\pi^{2}cos^{2}\beta},\;F_{u}=\frac{1}{16\pi^{2}}\left(\frac{9}{10}g_{1}^{2}+\frac{9}{2}g_{2}^{2}\right)$

but, for SM case with $\mu\geq m_{s}$,

$F_{\tau}=\frac{3h_{\tau}^{2}}{32\pi^{2}},\;F_{u}=\frac{1}{16\pi^{2}}\left(-\frac{9}{10}g_{1}^{2}-\frac{3}{2}g_{2}^{2}+6h_{b}^{2}-2\lambda\right).$

For scale-independent VEV in the case of MSSM with $\mu\geq m_{s}$,

$F_{\tau}=-\frac{h_{\tau}^{2}}{16\pi^{2}\cos^{2}\beta},\;F_{u}=\frac{1}{16\pi^{2}}\left(\frac{6}{5}g_{1}^{2}+6g_{2}^{2}-6\frac{h_{t}^{2}}{\sin\beta^{2}}\right)$

but, for SM case with $\mu\geq m_{s}$,

$F_{\tau}=\frac{3h_{\tau}^{2}}{32\pi^{2}},\;F_{u}=\frac{1}{16\pi^{2}}\left(3g_{2}^{2}-2\lambda-6h_{t}^{2}\right).$

For the present analysis, we adopt usual sign convention $|m_{2}|>|m_{1}|$ and
we shall use RGEs of scale dependent VEV, which are different from RGEs of
scale independent VEV in the expression of $F_{u}$ involved in the equations.
The general case of fermion masses which decrease with the increase in energy
scale is consistent with that of neutrino masses in the running of RGEs with
VEV das2005numerical ; singh2001effects . The RGEs of gauge and Yukawa
couplings with and without SUSY, are given in Appendix-A. The corresponding
RGEs for three mixing angles and three phases are given in Appendix-B.

##  III Numerical analysis of RGEs for GR neutrino mixing matrix

For a complete numerical analysis of the RGEs given in the above section, we
follow here two consecutive steps: (i) bottom-up running parida1998low  in the
first step and then (ii) top-down runningking2001inverted  in the next step.

In the first step (i), the running of the RGEs for the third family Yukawa
couplings ($h_{t},h_{b},h_{\tau}$) and three gauge couplings
($g_{1},g_{2},g_{3}$) is carried out from top quark mass scale,
$t_{0}=\ln(m_{t}/\text{1}{GeV})$ at low energy end to high energy scale
$M_{R}$ via SUSY breaking scale $m_{s}$ where ($m_{t}<m_{s}<M_{R}$).

At the transition point from SM to MSSM, the appropriate matching conditions
without threshold corrections are given as follows singh2018stability ,

$g_{i}(SUSY)=g_{i}(SM),$ (20) $h_{t}(SUSY)=\frac{h_{t}(SM)}{\sin\beta},$ (21)
$h_{b}(SUSY)=\frac{h_{b}(SM)}{\cos\beta},$ (22)
$h_{\tau}(SUSY)=\frac{h_{\tau}(SM)}{\cos\beta},$ (23)

where $\tan\beta$=$\frac{v_{u}}{v_{d}}$ such that $v_{u}=v_{o}\sin\beta$,
$v_{d}=v_{o}\cos\beta$ and $v_{o}=174$ GeV is the VEV of the Higgs field
alekhin2012top .

For large value of $\tan\beta$, there should be SUSY threshold corrections
which would lead to the modification of down-type quark and charged-lepton
Yukawa coupling constants at the matching condition of SUSY breaking scale
($m_{s}$) as follows antusch2013running ; zhang2016viability ;
antusch2008quark ,

$h_{t}(SUSY)\simeq\frac{h_{t}(SM)}{\sin\bar{\beta}},$ (24)
$h_{b}(SUSY)\simeq\left(\frac{1}{1+\bar{\eta_{b}}}\right)\frac{h_{b}(SM)}{\cos\bar{\beta}},$
(25) $h_{\tau}(SUSY)\simeq\frac{h_{\tau}(SM)}{\cos\bar{\beta}},$ (26)

where $\bar{\eta_{b}}$ is a free parameter that describes the SUSY threshold
corrections, cos$\bar{\beta}=(1+\eta_{l}^{{}^{\prime}})\cos\beta$ in the
redefinition of $\beta\rightarrow\bar{\beta}$ and $\eta_{l}^{{}^{\prime}}$ is
a leptonic SUSY threshold correction parameter which is typically very small.
Neglecting the effect of the leptonic threshold correction parameters in our
parametrisation, it would simply mean that tan$\bar{\beta}$=tan$\beta$.

The latest experimental input values for physical fermion masses, gauge
couplings and Weinberg mixing angle at electroweak scale ($m_{Z}$)pdg  are
given in Table 1.

Mass in GeV | Coupling constant | Weinberg mixing angle
---|---|---
$m_{z}$($m_{Z}$)=91.1876 |  $\alpha_{em}^{-1}$($m_{Z}$)=127.952 |  $sin^{2}\theta_{W}(m_{Z})$=0.23121
$m_{t}$($m_{t}$)=172.76 |  $\alpha_{s}(m_{Z}$)=0.1179 |
$m_{b}$($m_{b}$)=4.18 |  |
$m_{\tau}$($m_{\tau}$)=1.7768 |  |
Table 1: Low energy experimental values of fermion masses, gauge coupling
constants and Weinberg mixing angle.

The three gauge couplings, $\alpha_{1}(m_{Z})=0.016943$ ,
$\alpha_{2}(m_{Z})=0.033802$ and $\alpha_{3}(m_{Z})=0.1179$ at low energy
scale ($m{{}_{Z}}$), are calculated by using latest PDG data given in Table 1,
and SM matching relations

$\frac{1}{\alpha_{em}(m_{Z})}=\frac{3}{5}\frac{1}{\alpha_{1}(m_{Z})}+\frac{1}{\alpha_{2}(m_{Z})};$
(27) $sin^{2}\theta_{W}(m_{Z})=\frac{\alpha_{em}(m_{Z})}{\alpha_{2}(m_{Z})}.$
(28)

In terms of the normalized coupling constant $(g_{i})$, $\alpha_{i}$ can be
expressed as $g_{i}=\sqrt{4\pi\alpha_{i}}$, where $i=1,2,3$ and it represents
electromagnetic, weak and strong coupling constants respectively. We adopt the
standard procedure to get the values of gauge couplings at top-quark mass
scale from the experimental measurements at $m_{Z}$, using one-loop RGEs for
simplicity, assuming the existence of one-light Higgs doublet and five quark
flavours below $m_{t}$ scaleparida1998low ; singh2001effects .

The evolution equation of gauge coupling constants of one loop for energy
range $m_{Z}\leq\mu\leq m_{t}$ in SM is given by

$\frac{1}{\alpha_{i}(\mu)}=\frac{1}{\alpha_{i}(m_{Z})}-\frac{b_{i}}{2\pi}ln\left(\frac{\mu}{m_{Z}}\right)$
(29)

where

$b_{i}=\left(\frac{53}{10},-\frac{1}{2},-4.0\right).$

Similarly, the Yukawa couplings are also evaluated at top-quark mass scale
using QCD-QED rescaling factors ($\eta_{i}$) in the standard
fashionsingh2001effects , which are given by following relations.

$h_{t}(m_{t})=\frac{m_{t}(m_{t})}{v_{0}};$ (30)

$h_{b}(m_{t})=\frac{m_{b}(m_{b})}{v_{0}\eta_{b}};$ (31)
$h_{\tau}(m_{t})=\frac{m_{\tau}(m_{\tau})}{v_{0}\eta_{\tau}}.$ (32)

The value of QCD-QED rescaling factors ($\eta_{i}$) and vacuum expectation
($v_{0}$) of Higgs field are given by $\eta_{b}=1.53$, $\eta_{\tau}=1.015$ and
$v_{0}=174$ GeV respectivelydeshpande1994predictive ; nimai1998third .

In the second step (ii), the runnings of three neutrino mass eigenvalues
($m_{1},m_{2},m_{3}$), three neutrino mixing angles ($s_{23},s_{13},s_{12}$)
and three phases ($\delta,\alpha_{1},\alpha_{2}$) are carried out, together
with the running of gauge and Yukuwa couplings, from high energy seesaw scale
$t_{R}$(=ln ($M_{R}$/1GeV)) to low energy scale
$t_{o}(=\ln(m_{t}/\text{1GeV}))$ via SUSY breaking scale
$t_{s}(=\ln(m_{s}/\text{1GeV}))$. In this case, we use the values of gauge and
Yukawa couplings evaluated earlier at the scale $t_{R}$ from the first stage
running of RGEs in (i). In principle, one can evaluate neutrino masses, mixing
angles and phasesmurayama2002theory  at every point in the energy scale.

The present numerical analysis has six unknown arbitrary input values at high
energy seesaw scale consisting of three neutrino masses and three phases.
These input values would be suitably chosen so that we can get the desired low
energy values of neutrino parameters and latest Planck cosmological upper
bound on the sum of three neutrino masses. When we choose a set of mass
eigenvalues, there are two possible cases of mass hierarchy, namely (i) normal
hierarchy ($m_{3}>>m_{2}>m_{1}$) and (ii) inverted hierarchy
($m_{2}>m_{1}>>m_{3}$) as we generally set to $|m_{2}|>|m_{1}|$. The three
neutrino mixing angles ($s_{23},s_{13},s_{12}$) used at high energy seesaw
scale ($M_{R}$), are given by the golden ratio mixing matrix, which are
constant input values in all different high energy scales, while the three
neutrino mass eigenvalues ($m_{1},m_{2},m_{3}$) and phases
($\delta,\alpha_{1},\alpha_{2}$) are suitably chosen input values which may
give the desired values of neutrino oscillation parameters
($s_{23},s_{13},s_{12},m_{1},m_{2},m_{3},\delta,\alpha_{1},\alpha_{2}$) at low
energy scale after taking radiative corrections. The main concern in our work
is to satisfy the latest upper cosmological bound on the sum of absolute
neutrino masses, $\sum|m_{i}|<0.12$ eV lorenz2021reconstruction ;
choudhury2020updated  with the generation of reactor angle, $|U_{e3}|$ at low
energy scale.

##  IV Results and discussion

For top-down running of RGEs from high to low energy scale, Table 2 represents
the high scale input parameters of gauge and Yukawa coupling constants which
are already evaluated in the bottom-up approach, for running the neutrino
oscillation parameters. In the numerical analysis of neutrino oscillation
parameters along with phases including SUSY threshold corrections, there are
four free parameters namely, $M_{R}$, $\tan\beta$, $\bar{\eta_{b}}$ and
$m_{s}$ which may affect the values of coupling constants. While taking a set
of coupling constants at high energy scale ($M_{R}$), four possible cases are
considered, where three of the four parameters are set to be fixed while other
one is taken as variable. The values of the parameters are suitably chosen
within the certain limit of ranges for checking the data of the output results
for all cases. In addition to a particular set of coupling constants at high
energy scale, we have nine neutrino oscillation parameters namely, three mass
eigenvalues ($m_{1},-m_{2},m_{3}$), three neutrino mixing angles
($\theta_{23},\theta_{12},\theta_{13}$) and three phases
($\delta,\alpha_{1},\alpha_{2}$). The negative sign in mass eigenvalues are
possible due to the absorption of two Majorana phases in the mass eigenvalues
as diag($m_{1}e^{i\bar{\alpha_{1}}},m_{2}e^{i\bar{\alpha_{2}}},m_{3}$) where
we consider $\bar{\alpha_{1}}=\alpha_{1}$ and
$\bar{\alpha_{2}}=\alpha_{2}+\pi$. This negative sign in the mass eigenvalues
may help to prevent from the possible singularity that may arise in the
evolution of RGEs which has such $(m_{i}-m_{j})$ term in the denominator.
Since we are considering GR neutrino mixing matrix, the three neutrino mixing
angles at the high energy scale, are given by $s_{23}=0.70710,s_{13}=0$ and
$s_{12}=0.52573$ respectively, and these input values are the same in all
cases. Since the reactor mixing angle ($\theta_{13}$) is exactly zero, we take
it to have an extremely small non-zero value which would be able to solve the
asymptotic function in the RGEs of Dirac phase. Now, the unknown arbitrary
input values at high energy scale, are reduced to only six parameters i.e.,
three mass eigenvalues and three phases. These six arbitrary input values
defined at high energy scale, should be suitably chosen so that the output
results are compatible with the low energy experimental neutrino oscillation
datamc ; acero2022improved , including latest cosmological upper bound.

We have considered both normal and inverted hierarchical mass models for the
numerical analysis. In the case of normal hierarchical mass model, all the low
energy neutrino parameters are found to lie within $3\sigma$ range of NuFIT
datamc  with $\sum|m_{i}|<0.12$ eV as shown in Tables 3-6. We also check the
case of inverted hierarchical mass model which fails to give the low energy
neutrino oscillation parameters and $\sum|m_{i}|<0.12$ eV within the
experimental bounds. We also study the radiative generation of $U_{e3}$ with
initial conditions, $\Delta m_{21}^{2}$=0 at high energy scale, and a non-zero
value of $m_{3}$, but it fails to give low energy experimental values of
neutrino oscillation parameters. These results are not presented in the
present work. We observe that in both cases of normal and inverted
hierarchical models, all neutrino mass eigenvalues are slightly increased in
magnitude with the decrease in energy scale, whereas the atmospheric mixing
angle ($s_{23}$) and solar mixing angle ($s_{12}$) are slightly deviated from
the mixing angles at high energy seesaw scale i.e $\theta_{23}>45^{0}$ for NH
and $\theta_{23}<45^{0}$ for IH.

Our detailed numerical analysis shows that a larger value of $\tan\beta>$60
and high energy scale ($M_{R}$) are preferred in order to satisfy the latest
cosmological upper bound on the sum of three absolute neutrino mass
eigenvalues, $\sum|m_{i}|<0.12$ eV. This result requires an additional SUSY
threshold free parameter $\bar{\eta_{b}}$ in the range from -0.6 to +0.6
antusch2013running ; zhang2016viability  arising from the threshold
corrections of heavy SUSY particlesmuramatsu2020susy ; mohapatra2005threshold
; gupta2015renormalization ; hollik2014lifting . It is found that all the
neutrino oscillation parameters are consistent with low energy experimental
data for $\bar{\eta_{b}}$=0.01 at high scale $M_{R}=10^{15}$ GeV, large value
of $tan\beta=68$, and $m_{s}=1$ TeV. It is also observed that the negative
values of $\bar{\eta_{b}}$ for the larger values of $\tan\beta$, are not
feasible in calculating the values of coupling constants at high energy scale
in the normal hierarchical model. As a result, we discard the negative values
of $\bar{\eta_{b}}$ in our numerical analysis.

The four possible ways of numerical analysis for $\sum|m_{i}|<0.12$ eV based
on high energy scale ($M_{R}$), $\tan\beta$, SUSY breaking scale ($m_{s}$) and
SUSY threshold parameter ($\bar{\eta_{b}}$), are studied as follows :

(a) Case-I : Taking fixed input values $M_{R}=10^{15}$ GeV, $\tan\beta$=68,
$m_{s}$=1TeV, we vary with $\bar{\eta_{b}}=0.01,0.2,0.4,0.6$ and this case is
allowed only when $\bar{\eta_{b}}$=0.01. Results are presented in Table 3 and
Fig.1(a).

(b) Case-II : Taking fixed values $\tan\beta$=68, $m_{s}$=1TeV,
$\bar{\eta_{b}}$=0.01, we vary with
$M_{R}=10^{12}\text{GeV},M_{R}=10^{13}\text{GeV},M_{R}=10^{14}\text{GeV},M_{R}=10^{15}\text{GeV}$
and this case is allowed only when $M_{R}=10^{15}$ GeV. Results are presented
in Table 4 and Fig.1(b).

(c) Case-III : Taking fixed values $M_{R}=10^{15}$ GeV, $m_{s}$=1TeV,
$\bar{\eta_{b}}$=0.01, we vary with $\tan\beta=38,48,58,68$ and this case is
allowed only when $\tan\beta=68$. Results are presented in Table 5 and
Fig.1(c).

(d) Case-IV : Taking fixed values $M_{R}=10^{15}$ GeV, $\tan\beta$=68,
$\bar{\eta_{b}}$=0.01, we vary with
$m_{s}=1\text{TeV},5\text{TeV},10\text{TeV},14\text{TeV}$ and this case is
allowed only when $m_{s}$=1TeV. Results are presented in Table 6 and Fig.1(d).

The required values of coupling constants for various cases are given in Table
2. The main numerical results of our analysis on neutrino oscillation
parameters with three phases and SUSY threshold corrections, are given in
Tables 3-6. The numerical values in the upper halves of the tables represent
the high energy scale input values of neutrino oscillation parameters and the
numerical values in the lower halves of the Tables 3-6 represent the neutrino
oscillation parameters at low energy scale. From numerical analysis of Table
3, it indicates that a small value of $\bar{\eta_{b}}=0.01$ can accommodate
the latest cosmological upper bound on the sum of three absolute neutrino mass
eigenvalues $\sum|m_{i}|<0.12$ eV and hence, this particular value of
$\bar{\eta_{b}}=0.01$ is fixed for the remaining three possible ways. From the
numerical analysis of Table 4-6, we also observe that the cosmological upper
bound on the sum of three absolute neutrino mass eigenvalues
$\sum|m_{i}|<0.12$ eV and the desired values of neutrino oscillation
parameters at low energy scale can be achieved for the cases at
$M_{R}=10^{15}$ GeV , $\tan\beta=68$, $\bar{\eta_{b}}=0.01$ and $m_{s}=1$ TeV
in the variation of these parameters. It is also observed that the input
values of CP violating Dirac phase and two Majorana phases at high energy
scale, have significant effects in the evolution of neutrino mixing angles at
low energy scale. The best suitable high energy input values for these phase
parameters are respectively found to be $\delta=175^{0}$, $\alpha_{1}=2^{0}$
and $\alpha_{2}=0.5^{0}$ for all possible cases in our analysis. The variation
of $\sum|m_{i}|$ with $\bar{\eta_{b}}$, $M_{R}$, $\tan\beta$ and $m_{s}$ for
(a) Case-I, (b) Case-II, (c) Case-III and (d) Case-IV are respectively shown
in Fig: 1(a, b, c and d). We have also extended our analysis for other cases,
where we fix $M_{R}=10^{15}$ GeV, $\tan\beta=68$, $m_{s}$=5 TeV, 10 TeV and 14
TeV respectively with various values of $\bar{\eta_{b}}=0.01,0.2,0.4,0.6$. The
variations of $\sum|m_{i}|$ with $\bar{\eta_{b}}$ for each case is shown in
Fig.2. Only case with $m_{s}=1$TeV falls within the acceptable region but all
cases with higher $m_{s}$ are also acceptable if $\sum{|m_{i}|}<0.23$ eV
ade2016planck . Hence, the case for inverted hierarchical model is not
presented in this work as it fails to give latest Planck cosmological bound on
the sum of three absolute neutrino mass eigenvalues,$\sum{|m_{i}|}<0.12$ eV.

Other similar work in the literature zhang2016viability  emphasises the fact
that if the Planck 2015 cosmological bound $\sum|m_{i}|<0.23$ eVade2016planck
is taken into account, none of the three mixing patterns (BM, TBM, GR) can be
identified as lepton mixing matrix below the seesaw threshold under radiative
corrections. However, our present work shows the validity of the mixing
pattern based on GR, which is consistent with the latest Planck 2021
cosmological bound $\sum|m_{i}|<0.12$ eVlorenz2021reconstruction ;
choudhury2020updated  for low SUSY breaking scale, $m_{s}=1$TeV in the normal
hierarchy. As explained before, the numerical analysis in the present work is
carried out with specific input parameters viz. $\tan\beta=68$,
$M_{R}=10^{15}$GeV, $m_{s}=1$TeV and $\bar{\eta_{b}}=0.01$. For higher values
of $1\text{TeV}\leq m_{s}\leq 14\text{TeV}$ with other input parameters
$\tan\beta=68$, $M_{R}=10^{15}$GeV and $0.01\leq\bar{\eta_{b}}\leq 0.6$, the
validity of GR is still acceptable if the cosmological bound is relaxed up to
old 2015 Planck bound, $\sum|m_{i}|<0.23$ ade2016planck  as shown in Fig.2.

We consider the three phase parameters - two Majorana phases ($\alpha_{1}$,
$\alpha_{2}$), one Dirac CP-violating phase ($\delta$) along with one free
parameter $\bar{\eta_{b}}$ to determine the characteristics of SUSY threshold
corrections at the matching scale. Further, the work in Ref.zhang2016viability
establishes an important analytical correlation
$\Delta\theta_{23}\geq\theta_{13}\tan\theta_{12}^{0}$, where $\theta_{12}^{0}$
is the high scale input value, which induces a severe tension with the
observed $\theta_{23}$ and leads to exclusion of both GR and TBM at $3\sigma$
level if the cosmological upper bound on the sum of the three absolute masses
is taken into account. Although the above correlation is not explicitly shown
in the present work, the result of our numerical analysis still agrees with
it. However, our analysis shows the validity of GR under the most stringent
latest Planck cosmological bound, $\sum m_{i}<0.12$eV at high energy seesaw
scale $M_{R}=10^{15}$ GeV with larger value of $\tan\beta=68$ whose values are
beyond the range of inputs assigned in Ref.zhang2016viability .

(a) Case-I |  $\bar{\eta_{b}}$=0.01 |  $\bar{\eta_{b}}$=0.2 |  $\bar{\eta_{b}}$=0.4 |  $\bar{\eta_{b}}$=0.6
---|---|---|---|---
$g_{1}$ | 0.6686 | 0.67031 | 0.67068 | 0.67085
$g_{2}$ | 0.7 | 0.70264 | 0.70325 | 0.70354
$g_{3}$ | 0.72562 | 0.72722 | 0.72761 | 0.72779
$h_{t}$ | 0.96331 | 0.75129 | 0.70179 | 0.67866
$h_{b}$ | 2.12322 | 0.652 | 0.445 | 0.34961
$h_{\tau}$ | 2.63357 | 1.09202 | 0.91413 | 0.84265
$\lambda$ | 0.49032 | 0.49032 | 0.49032 | 0.49032
(b) Case-II |  ${M_{R}}=10^{12}$ GeV |  ${M_{R}}=10^{13}$ GeV |  ${M_{R}}=10^{14}$ GeV |  ${M_{R}}=10^{15}$ GeV
$g_{1}$ | 0.59652 | 0.61805 | 0.64197 | 0.6686
$g_{2}$ | 0.68535 | 0.69036 | 0.69532 | 0.7
$g_{3}$ | 0.78147 | 0.76158 | 0.74305 | 0.72562
$h_{t}$ | 0.8677 | 0.88399 | 0.9115 | 0.96331
$h_{b}$ | 1.19924 | 1.33557 | 1.57371 | 2.12322
$h_{\tau}$ | 1.26012 | 1.45802 | 1.80121 | 2.63357
$\lambda$ | 0.49032 | 0.49032 | 0.49032 | 0.49032
(c) Case-III |  $\tan\beta$=38 |  $\tan\beta$=48 |  $\tan\beta$=58 |  $\tan\beta$=68
$g_{1}$ | 0.67149 | 0.6712 | 0.67066 | 0.6686
$g_{2}$ | 0.70408 | 0.70365 | 0.70286 | 0.7
$g_{3}$ | 0.72792 | 0.72767 | 0.72722 | 0.72562
$h_{t}$ | 0.66369 | 0.69434 | 0.75274 | 0.96331
$h_{b}$ | 0.26635 | 0.39168 | 0.62803 | 2.12322
$h_{\tau}$ | 0.32940 | 0.47995 | 0.75944 | 2.63357
$\lambda$ | 0.49032 | 0.49032 | 0.49032 | 0.49032
(d) Case-IV |  $m_{s}$=1 TeV |  $m_{s}$=5 TeV |  $m_{s}$=10 TeV |  $m_{s}$=14 TeV
$g_{1}$ | 0.6686 | 0.66194 | 0.65883 | 0.65749
$g_{2}$ | 0.7 | 0.68679 | 0.68089 | 0.67839
$g_{3}$ | 0.72562 | 0.71030 | 0.70373 | 0.70097
$h_{t}$ | 0.96331 | 0.79242 | 0.75688 | 0.74448
$h_{b}$ | 2.12322 | 1.12484 | 0.99582 | 0.95473
$h_{\tau}$ | 2.63357 | 1.46905 | 1.32587 | 1.28061
$\lambda$ | 0.49032 | 0.47493 | 0.46891 | 0.46747
Table 2: High energy scale input values of gauge coupling ($g_{1},g_{2},g_{3}$), Yukawa coupling ($h_{t},h_{b},h_{\tau}$) and quartic coupling ($\lambda$) constants at various SUSY threshold free parameter ($\bar{\eta_{b}}$)((a) Case-I), high energy seesaw scale ($M_{R}$) ((b) Case-II), SUSY matching condition parameter ($\tan\beta$) ((c) Case-III) and SUSY breaking scale ($m_{s}$)((d) Case-IV) for four possible cases Parameter |  $\bar{\eta_{b}}$=0.01 |  $\bar{\eta_{b}}$=0.2 |  $\bar{\eta_{b}}$=0.4 |  $\bar{\eta_{b}}$=0.6
---|---|---|---|---
$m_{1}^{0}$[eV] | 0.02269 | 0.03408 | 0.03768 | 0.03945
$m_{2}^{0}$[eV] | -0.02580 | -0.03633 | -0.03974 | -0.04144
$m_{3}^{0}$[eV] | 0.04881 | 0.05281 | 0.05481 | 0.05582
$s_{23}^{0}$ | 0.70710 | 0.70710 | 0.70710 | 0.70710
$s_{13}^{0}$ | 0.0001 | 0.0001 | 0.0001 | 0.0001
$s_{12}^{0}$ | 0.52573 | 0.52573 | 0.52573 | 0.52573
$\delta^{0}$[$/^{0}$] | 175 | 175 | 175 | 175
$\alpha_{1}^{0}$[$/^{0}$] | 2 | 2 | 2 | 2
$\alpha_{2}^{0}$[$/^{0}$] | 0.5 | 0.5 | 0.5 | 0.5
$\Delta m_{21}^{2}[10^{-5}eV^{2}]$ | 15.08 | 15.84 | 15.94 | 16.09
$\Delta m_{31}^{2}[10^{-3}eV^{2}]$ | 1.86 | 1.62 | 1.58 | 1.55
$\sum|m_{i}|[eV]$ | 0.0973 | 0.12322 | 0.13223 | 0.13671
$m_{1}$[eV] | 0.03038 | 0.04760 | 0.05302 | 0.05569
$m_{2}$[eV] | -0.03158 | -0.04839 | -0.05371 | -0.05636
$m_{3}$[eV] | 0.05798 | 0.06883 | 0.07268 | 0.07457
$s_{23}$ | 0.78946 | 0.77451 | 0.77285 | 0.77220
$s_{13}$ | 0.14289 | 0.14319 | 0.14393 | 0.14424
$s_{12}$ | 0.53576 | 0.54040 | 0.54342 | 0.54336
$\delta$[$/^{0}$] | 201.19 | 206.15 | 208.5 | 208.21
$\alpha_{1}$[$/^{0}$] | 12.5 | 16.63 | 18.41 | 18.41
$\alpha_{2}$[$/^{0}$] | 4.61 | 6.26 | 6.99 | 6.99
$\Delta m_{21}^{2}[10^{-5}eV^{2}]$ | 7.46 | 7.54 | 7.3 | 7.46
$\Delta m_{31}^{2}[10^{-3}eV^{2}]$ | 2.43 | 2.47 | 2.47 | 2.45
$\sum|m_{i}|[eV]$ | 0.11994 | 0.16482 | 0.17941 | 0.18662
Table 3: (a) Case-I: Effect of variations with SUSY threshold parameter $\bar{\eta_{b}}$ on radiative generation of low energy neutrino parameters in MSSM while running from $M_{R}=10^{15}$ GeV to low scale $m_{t}=172.76$ GeV through $m_{s}$=1TeV for a particular value of $\tan\beta$=68. Only case for $\bar{\eta_{b}}$=0.01 is allowed  Parameter |  $M_{R}=10^{12}$GeV |  $M_{R}=10^{13}$GeV |  $M_{R}=10^{14}$GeV |  $M_{R}=10^{15}$GeV
---|---|---|---|---
$m_{1}^{0}$[eV] | 0.04235 | 0.03605 | 0.02985 | 0.02269
$m_{2}^{0}$[eV] | -0.04436 | -0.03854 | -0.03254 | -0.02580
$m_{3}^{0}$[eV] | 0.06112 | 0.05652 | 0.05252 | 0.04881
$s_{23}^{0}$ | 0.70710 | 0.70710 | 0.70710 | 0.70710
$s_{13}^{0}$ | 0.0001 | 0.0001 | 0.0001 | 0.0001
$s_{12}^{0}$ | 0.52573 | 0.52573 | 0.52573 | 0.52573
$\delta^{0}$[$/^{0}$] | 175 | 175 | 175 | 175
$\alpha_{1}^{0}$[$/^{0}$] | 2 | 2 | 2 | 2
$\alpha_{2}^{0}$[$/^{0}$] | 0.5 | 0.5 | 0.5 | 0.5
$\Delta m_{21}^{2}[10^{-5}eV^{2}]$ | 19.92 | 18.57 | 16.78 | 15.08
$\Delta m_{31}^{2}[10^{-3}eV^{2}]$ | 1.94 | 1.89 | 1.86 | 1.86
$\sum|m_{i}|[eV]$ | 0.14811 | 0.13111 | 0.11491 | 0.09730
$m_{1}$[eV] | 0.05343 | 0.04667 | 0.03951 | 0.03038
$m_{2}$[eV] | -0.05414 | -0.04750 | -0.04043 | -0.03158
$m_{3}$[eV] | 0.07276 | 0.06813 | 0.06357 | 0.05798
$s_{23}$ | 0.77231 | 0.77532 | 0.78007 | 0.78946
$s_{13}$ | 0.14320 | 0.14411 | 0.14450 | 0.14289
$s_{12}$ | 0.54104 | 0.53932 | 0.53923 | 0.53576
$\delta$[$/^{0}$] | 208.72 | 206.89 | 206.6 | 201.19
$\alpha_{1}$[$/^{0}$] | 17.11 | 15.96 | 15.66 | 12.5
$\alpha_{2}$[$/^{0}$] | 6.45 | 5.98 | 5.87 | 4.61
$\Delta m_{21}^{2}[10^{-5}eV^{2}]$ | 7.59 | 7.74 | 7.28 | 7.46
$\Delta m_{31}^{2}[10^{-3}eV^{2}]$ | 2.43 | 2.46 | 2.47 | 2.43
$\sum|m_{i}|[eV]$ | 0.18033 | 0.16230 | 0.14351 | 0.11994
Table 4: (b) Case-II: Effect of variation with high energy seesaw scale $M_{R}$ on radiative generation of low energy neutrino parameters in MSSM while running from $M_{R}$ GeV to low scale $m_{t}=172.76$ GeV for a particular values of $\tan\beta$=68, $\bar{\eta_{b}}=0.01$ and $m_{s}$=1TeV. Only case for $M_{R}=10^{15}$ GeV is allowed. Parameter |  $\tan\beta$=68 |  $\tan\beta$=58 |  $\tan\beta$=48 |  $\tan\beta$=38
---|---|---|---|---
$m_{1}^{0}$[eV] | 0.02269 | 0.04752 | 0.06752 | 0.09082
$m_{2}^{0}$[eV] | -0.02580 | -0.04930 | -0.06882 | -0.09177
$m_{3}^{0}$[eV] | 0.04881 | 0.06181 | 0.07781 | 0.09821
$s_{23}^{0}$ | 0.70710 | 0.70710 | 0.70710 | 0.70710
$s_{13}^{0}$ | 0.0001 | 0.0001 | 0.0001 | 0.0001
$s_{12}^{0}$ | 0.52573 | 0.52573 | 0.52573 | 0.52573
$\delta^{0}$[$/^{0}$] | 175 | 175 | 175 | 175
$\alpha_{1}^{0}$[$/^{0}$] | 2 | 2 | 2 | 2
$\alpha_{2}^{0}$[$/^{0}$] | 0.5 | 0.5 | 0.5 | 0.5
$\Delta m_{21}^{2}[10^{-5}eV^{2}]$ | 15.08 | 17.23 | 17.72 | 17.34
$\Delta m_{31}^{2}[10^{-3}eV^{2}]$ | 1.86 | 1.56 | 1.49 | 1.39
$\sum|m_{i}|[eV]$ | 0.0973 | 0.15863 | 0.21415 | 0.2808
$m_{1}$[eV] | 0.03038 | 0.06764 | 0.09749 | 0.13210
$m_{2}$[eV] | -0.03158 | -0.06820 | -0.09788 | -0.13239
$m_{3}$[eV] | 0.05798 | 0.08432 | 0.10981 | 0.14108
$s_{23}$ | 0.78946 | 0.77197 | 0.76960 | 0.76817
$s_{13}$ | 0.14289 | 0.14857 | 0.14755 | 0.14530
$s_{12}$ | 0.53576 | 0.54926 | 0.55597 | 0.55416
$\delta$[$/^{0}$] | 201.19 | 213.25 | 216.8 | 214.85
$\alpha_{1}$[$/^{0}$] | 12.5 | 21.35 | 24.29 | 23.59
$\alpha_{2}$[$/^{0}$] | 4.61 | 8.23 | 9.51 | 9.2
$\Delta m_{21}^{2}[10^{-5}eV^{2}]$ | 7.46 | 7.59 | 7.52 | 7.63
$\Delta m_{31}^{2}[10^{-3}eV^{2}]$ | 2.43 | 2.53 | 2.55 | 2.45
$\sum|m_{i}|[eV]$ | 0.11994 | 0.22016 | 0.30518 | 0.40557
Table 5: (c) Case-III: Effect of variation with the values of tan$\beta$ on radiative generation of low energy neutrino parameters in MSSM while running from $M_{R}=10^{15}$ GeV to low scale $m_{t}=172.76$ GeV for a particular value of $m_{s}$=1TeV and $\bar{\eta_{b}}=0.01$. Only case for tan$\beta$=68 is allowed. Parameter |  $m_{s}$=1TeV |  $m_{s}$=5TeV |  $m_{s}$=10TeV |  $m_{s}$=14TeV
---|---|---|---|---
$m_{1}^{0}$[eV] | 0.02269 | 0.03018 | 0.03194 | 0.03304
$m_{2}^{0}$[eV] | -0.02580 | -0.03271 | -0.03436 | -0.03546
$m_{3}^{0}$[eV] | 0.04881 | 0.05131 | 0.05222 | 0.05322
$s_{23}^{0}$ | 0.70710 | 0.70710 | 0.70710 | 0.70710
$s_{13}^{0}$ | 0.0001 | 0.0001 | 0.0001 | 0.0001
$s_{12}^{0}$ | 0.52573 | 0.52573 | 0.52573 | 0.52573
$\delta^{0}$[$/^{0}$] | 175 | 175 | 175 | 175
$\alpha_{1}^{0}$[$/^{0}$] | 2 | 2 | 2 | 2
$\alpha_{2}^{0}$[$/^{0}$] | 0.5 | 0.5 | 0.5 | 0.5
$\Delta m_{21}^{2}[10^{-5}eV^{2}]$ | 15.08 | 15.91 | 16.04 | 16.57
$\Delta m_{31}^{2}[10^{-3}eV^{2}]$ | 1.86 | 1.72 | 1.70 | 1.74
$\sum|m_{i}|[eV]$ | 0.0973 | 0.1142 | 0.11852 | 0.12172
$m_{1}$[eV] | 0.03038 | 0.04109 | 0.04346 | 0.04492
$m_{2}$[eV] | -0.03158 | -0.04199 | -0.0443 | -0.04576
$m_{3}$[eV] | 0.05798 | 0.06424 | 0.06575 | 0.06711
$s_{23}$ | 0.78946 | 0.77816 | 0.77639 | 0.77600
$s_{13}$ | 0.14289 | 0.14360 | 0.14308 | 0.14347
$s_{12}$ | 0.53576 | 0.53853 | 0.53949 | 0.53964
$\delta$[$/^{0}$] | 201.19 | 204.94 | 205.84 | 205.92
$\alpha_{1}$[$/^{0}$] | 12.5 | 15.29 | 15.99 | 16.11
$\alpha_{2}$[$/^{0}$] | 4.61 | 5.71 | 6 | 6.04
$\Delta m_{21}^{2}[10^{-5}eV^{2}]$ | 7.46 | 7.48 | 7.35 | 7.58
$\Delta m_{31}^{2}[10^{-3}eV^{2}]$ | 2.43 | 2.43 | 2.43 | 2.48
$\sum|m_{i}|[eV]$ | 0.11994 | 0.14732 | 0.15351 | 0.15779
Table 6: (d) Case-IV: Effect of the variation with SUSY breaking scale $m_{s}$
on radiative generation of low energy neutrino parameters in MSSM while
running from $M_{R}=10^{15}$ GeV to low scale $m_{t}=172.76$ GeV for a
particular value of $\tan\beta$=68 and $\bar{\eta_{b}}=0.01$. Only case for
$m_{s}$=1TeV is allowed

Figure 1: Graphical presentation of the results given in Tables 3-6. (a)
Case-I: Variation of $\sum|m_{i}|[\text{eV}]$ with $\bar{\eta_{b}}$ for
$M_{R}=10^{15}\text{GeV},\tan\beta=68,m_{s}=1\text{TeV}$, (b) Case-II:
Variation of $\sum|m_{i}|[\text{eV}]$ with $M_{R}$ for
$\bar{\eta_{b}}=0.01,\tan\beta=68,m_{s}=1\text{TeV}$, (c) Case-III: Variation
of $\sum|m_{i}|[\text{eV}]$ with $\tan\beta$ for
$\bar{\eta_{b}}=0.01,M_{R}=10^{15}\text{GeV},m_{s}=1\text{TeV}$, (d) Case-IV:
Variation of $\sum|m_{i}|[\text{eV}]$ with $m_{s}$ for
$\bar{\eta_{b}}=0.01,\tan\beta=68,M_{R}=10^{15}\text{GeV}$

.

Figure 2: Graphical presentation for variations of the sum of three absolute
neutrino mass eigenvalues ($\sum|m_{i}|[\text{eV}]$) with the various values
of SUSY threshold parameter ($\bar{\eta_{b}}=0.01,0.2,0.4,0.6$) for different
cases of SUSY breaking scale $m_{s}$=1 TeV, 5 TeV, 10 TeV and 14 TeV. Values
of $M_{R}=10^{15}$GeV and $\tan\beta=68$ are taken. Higher upper bound
$\sum|m_{i}|<0.23$ eV can accommodate at wide range of parameters,
$\bar{\eta_{b}}$=(0.01-0.6) and $m_{s}$=(1TeV-14TeV)

.

##  V Summary and conclusion

We summarise the main points of the present work related to the implication of
the latest Planck data on the cosmological upper bound on the sum of the three
absolute neutrino masses.

Our numerical analysis is based on the evolution of RGEs of neutrino
oscillation parameters and three phases, including the effect of scale-
dependent VEV and SUSY threshold corrections. We have first found out the most
suitable value of SUSY threshold parameter $\bar{\eta_{b}}$ in the range from
-0.6 to +0.6 antusch2013running ; zhang2016viability , which is compatible
with the low energy neutrino oscillation parameters and the most stringent
cosmological upper bound, $\sum|m_{i}|<0.12$ eV lorenz2021reconstruction ;
choudhury2020updated . It has been observed that the negative values of
$\bar{\eta_{b}}$ are not feasible for large values of $\tan\beta$ and $M_{R}$
in the normal hierarchical model. The best fitted value of $\bar{\eta_{b}}$
that satisfies the cosmological upper bound $\sum|m_{i}|<0.12$ eV is found to
be 0.01 when $M_{R}=10^{15}$GeV, $\tan\beta$=68 and $m_{s}$=1 TeV.

The detailed numerical analysis shows that all neutrino mass eigenvalues as
well as mixing angles are increased with the decrease in energy scale and it
gives certain advantages on GR mixing matrix over the other two mixing
matrices such as BM and TBM in normal hierarchical mass model. It is also
observed that the other neutrino mixing patterns such as BM and TBM mixings,
do not satisfy the above cosmological upper bound in both normal hierarchical
and inverted hierarchical models. The low energy neutrino oscillation
parameters along with the latest Planck cosmological upper bound on sum of
three absolute mass eigenvalues $\sum|m_{i}|<0.12$ eV can be achieved through
radiative corrections under RGEs using GR neutrino mixing matrix defined at
high energy seesaw scale for the input values $M_{R}=10^{15}$ GeV,
$\tan\beta$=68, $m_{s}$=1 TeV and $\bar{\eta_{b}}=0.01$. However, it is
observed that a larger upper bound of $\sum|m_{i}|<0.23$ eV ade2016planck  can
accommodate at wide range of parameters: $\bar{\eta_{b}}$=(0.01-0.6) and
$m_{s}$=(1TeV-14TeV) as evident from Fig.2.

The analysis of inverted hierarchical neutrino mass model shows that the model
does not accommodate the latest cosmological upper bound but it still conforms
with earlier 2015 Planck bound $\sum|m_{i}|<0.23$ eV ade2016planck . Further
analysis for TBM case chen2019cp ; chen2019c , considering the effect of CP
violating phases and SUSY threshold corrections, will be reported in future
communication.

To conclude, the present investigation indicates the sensitivity of the value
of $\sum|m_{i}|$ on the origin of neutrino masses and mixing angles. It is
relevant in the context of the information related to the absolute neutrino
masses that has been continuously updating with recent Planck data on the
cosmological upper bound on the sum of three absolute neutrino masses
$\sum|m_{i}|<0.12$ eV. Neutrino mass model if any, is bound to be consistent
with these upper bounds on absolute neutrino masses. While the existence of
supersymmetric particles has been continuously ruling out in LHC, the
supersymmetric breaking scale ($m_{s}$) still remains as an unknown parameter.
We assume that the $m_{s}$ scale may lie somewhere in between 1 TeV and 14 TeV
within the scope of LHC, and the present work is thus confined to the
implication of SUSY breaking scale. It is a continuation of our previous
investigationwilina2022deviations ; devi2022effects ; singh2018stability  on
neutrino masses and mixings with varying SUSY breaking scale in the running of
RGEs in both normal and inverted hierarchical neutrino mass models.

The focus of the present work is the question of the validity of GR neutrino
mixing at high energy scale, with the variation of $m_{s}$ scale and other
input parameters $\tan\beta$, $\bar{\eta_{b}}$ and $M_{R}$ scale. It has
profound implications to apply on other aspects of RGEs analysis such as low
energy magnification of neutrino mixings in quark-lepton unification
hypothesis at high energy scale in SO(10) modelagarwalla2007neutrino ;
srivastava2016predictions ; abdussalam2021majorana , radiative generation of
reactor mixing angle and solar neutrino mass squared difference at low scale,
and the question of radiative stability of neutrino mass models to
discriminate between NH and IH models. These earlier good results may now be
readdressed for further analysis at low energy scale, consistent with latest
Planck data on cosmological upper bound on the sum of three absolute mass
eigenvalues.

## Acknowledgements

One of the authors (YMS) would like to thank Manipur University for granting
Fellowship for Ph.D course.

## Appendix A

The one-loop RGEs for Yukawa couplings in the MSSM in the range of mass scales
$m_{s}\leq\mu\leq M_{R}$ are given by rajkhowa2006possible

$\frac{d}{dt}h_{t}=\frac{h_{t}}{16\pi^{2}}\left(6h_{t}^{2}+h_{b}^{2}-\sum_{i=1}^{3}c_{i}g_{i}^{2}\right)$
(33)

,

$\frac{d}{dt}h_{b}=\frac{h_{b}}{16\pi^{2}}\left(6h_{b}^{2}+h_{\tau}^{2}+h_{t}^{2}-\sum_{i=1}^{3}c_{i}^{{}^{\prime}}g_{i}^{2}\right),$
(34)
$\frac{d}{dt}h_{\tau}=\frac{h_{\tau}}{16\pi^{2}}\left(4h_{\tau}^{2}+3h_{b}^{2}-\sum_{i=1}^{3}c_{i}^{{}^{\prime\prime}}g_{i}^{2}\right)$
(35)

where, for SUSY case,

$c_{i}=\left(\begin{array}[]{r}\frac{13}{15},3,\frac{16}{3}\end{array}\right),c_{i}^{{}^{\prime}}=\left(\begin{array}[]{r}\frac{7}{15},3,\frac{16}{3}\end{array}\right),c_{i}^{{}^{\prime\prime}}=\left(\begin{array}[]{r}\frac{9}{5},3,0\end{array}\right).$

The Yukawa RGEs, for non-SUSY(SM) in the range of mass scales
$m_{t}\leq\mu\leq m_{s}$,

$\frac{d}{dt}h_{t}=\frac{h_{t}}{16\pi^{2}}\left(\frac{9}{2}h_{t}^{2}+\frac{3}{2}h_{b}^{2}+h_{\tau}^{2}-\sum_{i=1}^{3}c_{i}g_{i}^{2}\right),$
(36)
$\frac{d}{dt}h_{b}=\frac{h_{b}}{16\pi^{2}}\left(\frac{9}{2}h_{b}^{2}+h_{\tau}^{2}+\frac{3}{2}h_{t}^{2}-\sum_{i=1}^{3}c_{i}^{{}^{\prime}}g_{i}^{2}\right),$
(37)
$\frac{d}{dt}h_{\tau}=\frac{h_{\tau}}{16\pi^{2}}\left(\frac{5}{2}h_{\tau}^{2}+3h_{b}^{2}+3h_{t}^{2}-\sum_{i=1}^{3}c_{i}^{{}^{\prime\prime}}g_{i}^{2}\right)$
(38)

where, for non-SUSY(SM) case,

$c_{i}=\left(\begin{array}[]{r}0.85,2.25,8.00\end{array}\right),c_{i}^{{}^{\prime}}=\left(\begin{array}[]{r}0.25,2.25,8.00\end{array}\right),c_{i}^{{}^{\prime\prime}}=\left(\begin{array}[]{r}2.25,2.25,0.0\end{array}\right).$

And one loop RGE for quartic coupling in SM is given by

$\frac{d\lambda}{dt}=\frac{1}{16\pi^{2}}[\frac{9}{4}(\frac{3}{25}g_{1}^{4}+\frac{2}{5}g_{1}^{2}g_{2}^{2}+g_{2}^{4})$
$-(\frac{9}{5}g_{1}^{2}+9g_{2}^{2})\lambda+4Y_{2}(S)\lambda-4H(S)+12\lambda_{2}]$
(39)

where,

$Y_{2}(S)=3h_{t}^{2}+3h_{b}^{2}+h_{\tau}^{2},$
$H(S)=3h_{t}^{4}+3h_{b}^{4}+h_{\tau}^{4},$ $\lambda=\frac{m_{h}^{2}}{v^{2}}$

$m_{h}$=Higgs mass, $v_{0}$=vacuum expectation value.

The two-loop RGEs for the gauge couplings are similarly expressed in the range
of mass scales $m_{s}\leq\mu\leq M_{R}$ as singh2018stability

$\frac{d}{dt}g_{i}=\frac{g_{i}}{16\pi^{2}}\left[b_{i}g_{i}^{2}+\frac{1}{16\pi^{2}}\left(\sum_{j=1}^{3}b_{ij}g_{i}^{3}g_{j}^{2}-\sum_{j=t,b,\tau}a_{ij}g_{i}^{3}h_{j}^{2}\right)\right]$
(40)

where, for SUSY case,

$b_{i}=\left(\begin{array}[]{r}6.6,1,-3\end{array}\right),b_{ij}=\left(\begin{array}[]{ccc}7.96&5.40&17.60\\\
1.80&25.00&24.00\\\
2.20&9.00&14.00\end{array}\right),a_{ij}=\left(\begin{array}[]{ccc}5.2&2.8&3.6\\\
6.0&6.0&2.0\\\ 4.0&4.0&0.0\end{array}\right)$

and, for the non-SUSY(SM) in the range of mass scales $m_{t}\leq\mu\leq
m_{s}$,

$b_{i}=\left(\begin{array}[]{r}4.1,-3.167,-7\end{array}\right),b_{ij}=\left(\begin{array}[]{ccc}3.98&2.70&8.8\\\
0.90&5.83&12.0\\\
1.10&4.50&-26.0\end{array}\right),a_{ij}=\left(\begin{array}[]{ccc}0.85&0.5&0.5\\\
1.50&1.5&0.5\\\ 2.0&2.0&0.0\end{array}\right).$

## Appendix B

The RGEs of neutrino mixing angles and CP violating Dirac and Majorana phases
are same for both scale-dependent VEV and scale independent VEV and are given
by following equationsagarwalla2007neutrino ; antusch2003running  ,

$\displaystyle\frac{ds_{12}}{dt}=\frac{F_{\tau}c_{12}\sin
2\theta_{12}s_{23}^{2}}{2(m_{2}^{2}-m_{1}^{2})}[m_{1}^{2}+m_{2}^{2}+2m_{1}m_{2}\cos(2\alpha_{2}-2\alpha_{1})],$
(41) $\displaystyle\frac{ds_{23}}{dt}=$ $\displaystyle\frac{F_{\tau}c_{23}\sin
2\theta_{23}}{2(m_{3}^{2}-m_{2}^{2})}[c_{12}^{2}(m_{3}^{2}+m_{2}^{2}+2m_{3}m_{2}\cos
2\alpha_{2})$ $\displaystyle+s_{12}^{2}(m_{3}^{2}+m_{1}^{2}+2m_{3}m_{1}\cos
2\alpha_{1})/(1+R)],$ (42) $\displaystyle\frac{ds_{13}}{dt}=$
$\displaystyle-\frac{F_{\tau}c_{13}\sin 2\theta_{12}\sin
2\theta_{23}m_{3}}{2(m_{3}^{2}-m_{1}^{2}}[m_{1}\cos(2\alpha_{1}-\delta)$
$\displaystyle-(1+R)m_{2}\cos(2\alpha_{2}-\delta)-Rm_{3}\cos\delta],$ (43)
$\displaystyle\frac{d\delta}{dt}=$ $\displaystyle-\frac{F_{\tau}m_{3}\sin
2\theta_{12}\sin
2\theta_{23}}{2\theta_{13}(m_{3}^{2}-m_{1}^{2})}\times[m_{1}\sin(2\alpha_{1}-\delta)-(1+R)m_{2}\sin(2\alpha_{2}-\delta)$
$\displaystyle+Rm_{3}\sin\delta]-2F_{\tau}[\frac{m_{1}m_{2}s_{23}^{2}\sin(2\alpha_{1}-2\alpha_{2})}{(m_{2}^{2}-m_{1}^{2})}+m_{3}s_{12}^{2}(\frac{m_{1}\cos
2\theta_{23}\sin 2\alpha_{1}}{(m_{3}^{2}-m_{1}^{2})}$
$\displaystyle+\frac{m_{2}c_{23}^{2}\sin(2\delta-2\alpha_{2})}{(m_{3}^{2}-m_{2}^{2})})+m_{3}c_{12}^{2}(\frac{m_{1}c_{23}^{2}\sin(2\delta-2\alpha_{1})}{(m_{3}^{2}-m_{1}^{2})}$
$\displaystyle+\frac{m_{2}\cos 2\theta_{23}\sin
2\alpha_{2}}{(m_{3}^{2}-m_{2}^{2})})],$ (44)
$\displaystyle\frac{d\alpha_{1}}{dt}=$ $\displaystyle-2F_{\tau}[m_{3}\cos
2\theta_{23}\frac{m_{1}s_{12}^{2}\sin\alpha_{1}+(1+R)m_{2}c_{12}^{2}\sin
2\alpha_{2}}{m_{3}^{2}-m_{1}^{2}}$
$\displaystyle+\frac{m_{1}m_{2}c_{12}^{2}s_{23}^{2}\sin(2\alpha_{1}-\alpha_{2})}{m_{2}^{2}-m_{1}^{2}}],$
(45) $\displaystyle\frac{d\alpha_{2}}{dt}=$ $\displaystyle-2F_{\tau}[m_{3}\cos
2\theta_{23}\frac{m_{1}s_{12}^{2}\sin\alpha_{1}+(1+R)m_{2}c_{12}^{2}\sin
2\alpha_{2}}{m_{3}^{2}-m_{1}^{2}}$
$\displaystyle+\frac{m_{1}m_{2}s_{12}^{2}s_{23}^{2}\sin(2\alpha_{1}-\alpha_{2})}{m_{2}^{2}-m_{1}^{2}}],$
(46)

where,

$R=\frac{(m_{2}^{2}-m_{1}^{2})}{(m_{3}^{2}-m_{2}^{2})}$

## References

  * (1) I. Esteban, M. C. González-García, M. Maltoni, T. Schwetz, and A. Zhou, “The fate of hints: updated global analysis of three-flavor neutrino oscillations,” Journal of High Energy Physics, vol. 2020, no. 9, pp. 1–22, 2020.
  * (2) P. F. de Salas, D. Forero, S. Gariazzo, P. Martínez-Miravé, O. Mena, C. Ternes, M. Tórtola, and J. Valle, “2020 global reassessment of the neutrino oscillation picture,” Journal of High Energy Physics, vol. 2021, no. 2, pp. 1–36, 2021.
  * (3) M. C. Gonzalez-Garcia, M. Maltoni, and T. Schwetz, “Nufit: Three-flavour global analyses of neutrino oscillation experiments,” Universe, vol. 7, no. 12, p. 459, 2021.
  * (4) N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. Banday, R. Barreiro, N. Bartolo, S. Basak, et al., “Planck 2018 results-vi. cosmological parameters,” Astronomy & Astrophysics, vol. 641, p. A6, 2020.
  * (5) I. Tanseri, S. Hagstotz, S. Vagnozzi, E. Giusarma, and K. Freese, “Updated neutrino mass constraints from galaxy clustering and cmb lensing-galaxy cross-correlation measurements,” Journal of High Energy Astrophysics, vol. 36, pp. 1–26, 2022.
  * (6) M. K. Das, M. Patgiri, and N. Nimai Singh, “Numerical consistency check between two approaches to radiative corrections for neutrino masses and mixings,” Pramana, vol. 65, no. 6, pp. 995–1013, 2005.
  * (7) S. King and N. N. Singh, “Renormalisation group analysis of single right-handed neutrino dominance,” Nuclear Physics B, vol. 591, no. 1-2, pp. 3–25, 2000.
  * (8) P. H. Chankowski and Z. Płuciennik, “Renormalization group equations for seesaw neutrino masses,” Physics Letters B, vol. 316, no. 2-3, pp. 312–317, 1993.
  * (9) K. Babu, C. N. Leung, and J. Pantaleone, “Renormalization of the neutrino mass operator,” Physics Letters B, vol. 319, no. 1-3, pp. 191–198, 1993.
  * (10) S. Antusch, M. Drees, J. Kersten, M. Lindner, and M. Ratz, “Neutrino mass operator renormalization in two higgs doublet models and the mssm,” Physics Letters B, vol. 525, no. 1-2, pp. 130–134, 2002.
  * (11) M. Parida and N. N. Singh, “Low-energy formulas for neutrino masses with a tan $\beta$-dependent hierarchy,” Physical Review D, vol. 59, no. 3, p. 032002, 1998.
  * (12) S. King and N. N. Singh, “Inverted hierarchy models of neutrino masses,” Nuclear Physics B, vol. 596, no. 1-2, pp. 81–98, 2001.
  * (13) P. H. Chankowski, W. Krolikowski, and S. Pokorski, “Fixed points in the evolution of neutrino mixings,” Physics Letters B, vol. 473, no. 1-2, pp. 109–117, 2000.
  * (14) A. Brignole, J. Casas, J. R. Espinosa, and I. Navarro, “Low-scale supersymmetry breaking: Effective description, electroweak breaking and phenomenology,” Nuclear Physics B, vol. 666, no. 1-2, pp. 105–143, 2003\.
  * (15) S. Antusch, J. Kersten, M. Lindner, and M. Ratz, “Running neutrino masses, mixings and cp phases: Analytical results and phenomenological consequences,” Nuclear Physics B, vol. 674, no. 1-2, pp. 401–433, 2003\.
  * (16) Y. Kajiyama, M. Raidal, and A. Strumia, “Golden ratio prediction for solar neutrino mixing,” Physical Review D, vol. 76, dec 2007.
  * (17) J. Zhang and S. Zhou, “Viability of exact tri-bimaximal, golden-ratio and bimaximal mixing patterns and renormalization-group running effects,” Journal of High Energy Physics, vol. 2016, no. 9, pp. 1–27, 2016.
  * (18) S. Pramanick, “Radiative generation of realistic neutrino mixing with a4,” Nuclear Physics B, vol. 963, p. 115282, 2021.
  * (19) S. King, “Unified models of neutrinos, flavour and cp violation,” Progress in Particle and Nuclear Physics, vol. 94, pp. 217–256, 2017.
  * (20) L. L. Everett and A. J. Stuart, “Icosahedral (a 5) family symmetry and the golden ratio prediction for solar neutrino mixing,” Physical Review D, vol. 79, no. 8, p. 085005, 2009.
  * (21) T. Fukuyama, “Twenty years after the discovery of $\mu$–$\tau$ symmetry,” Progress of Theoretical and Experimental Physics, vol. 2017, no. 3, 2017\.
  * (22) Y. Koide, “Universal texture of quark and lepton mass matrices with an extended flavor 2↔ 3 symmetry,” Physical Review D, vol. 69, no. 9, p. 093001, 2004.
  * (23) C. S. Lorenz, L. Funcke, M. Löffler, and E. Calabrese, “Reconstruction of the neutrino mass as a function of redshift,” Physical Review D, vol. 104, no. 12, p. 123518, 2021.
  * (24) S. R. Choudhury and S. Hannestad, “Updated results on neutrino mass and mass hierarchy from cosmology with planck 2018 likelihoods,” Journal of Cosmology and Astroparticle Physics, vol. 2020, no. 07, p. 037, 2020.
  * (25) E. Ma, “All-purpose neutrino mass matrix,” Physical Review D, vol. 66, no. 11, p. 117301, 2002.
  * (26) S. Antusch and S. F. King, “Charged lepton corrections to neutrino mixing angles and cp phases revisited,” Physics Letters B, vol. 631, no. 1-2, pp. 42–47, 2005.
  * (27) S. Antusch, J. Kersten, M. Lindner, and M. Ratz, “Neutrino mass matrix running for non-degenerate see-saw scales,” Physics Letters B, vol. 538, no. 1-2, pp. 87–95, 2002.
  * (28) A. Davidson and K. C. Wali, “Universal seesaw mechanism?,” Physical Review Letters, vol. 59, no. 4, p. 393, 1987.
  * (29) J. Casas, J. R. Espinosa, A. Ibarra, and I. Navarro, “Nearly degenerate neutrinos, supersymmetry and radiative corrections,” Nuclear Physics B, vol. 569, no. 1-3, pp. 82–106, 2000.
  * (30) H. Arason, D. Castano, B. Kesthelyi, S. Mikaelian, E. Piard, P. Ramond, and B. Wright, “Renormalization-group study of the standard model and its extensions: the standard model,” Physical Review D, vol. 46, no. 9, p. 3945, 1992.
  * (31) S. K. Agarwalla, M. Parida, R. Mohapatra, and G. Rajasekaran, “Neutrino mixings and leptonic c p violation from ckm matrix and majorana phases,” Physical Review D, vol. 75, no. 3, p. 033007, 2007.
  * (32) N. N. Singh, “Effects of the scale-dependent vacuum expectation values in the renormalisation group analysis of neutrino masses,” The European Physical Journal C-Particles and Fields, vol. 19, no. 1, pp. 137–141, 2001.
  * (33) K. S. Singh, S. Roy, and N. N. Singh, “Stability of neutrino parameters and self-complementarity relation with varying susy breaking scale,” Physical Review D, vol. 97, no. 5, p. 055038, 2018.
  * (34) S. Alekhin, A. Djouadi, and S. Moch, “The top quark and higgs boson masses and the stability of the electroweak vacuum,” Physics Letters B, vol. 716, no. 1, pp. 214–219, 2012.
  * (35) S. Antusch and V. Maurer, “Running quark and lepton parameters at various scales,” Journal of High Energy Physics, vol. 2013, no. 11, pp. 1–17, 2013\.
  * (36) S. Antusch and M. Spinrath, “Quark and lepton masses at the gut scale including supersymmetric threshold corrections,” Physical Review D, vol. 78, no. 7, p. 075020, 2008.
  * (37) P. D. Group, P. A. Zyla, R. M. Barnett, and e. a. Beringer, “Review of Particle Physics,” Progress of Theoretical and Experimental Physics, vol. 2020, no. 8, 2020.  083C01.
  * (38) N. Deshpande and E. Keith, “Predictive fermion mass matrix ansätze in nonsupersymmetric so (10) grand unification,” Physical Review D, vol. 50, no. 5, p. 3513, 1994.
  * (39) N. Nimai Singh and S. Biramani Singh, “Third generation yukawa couplings unification in supersymmetric so (10) model,” The European Physical Journal C-Particles and Fields, vol. 5, no. 2, pp. 363–367, 1998.
  * (40) H. Murayama, “Theory of neutrino masses and mixings,” International Journal of Modern Physics A, vol. 17, no. 24, pp. 3403–3420, 2002.
  * (41) M. Acero, P. Adamson, L. Aliaga, N. Anfimov, A. Antoshkin, E. Arrieta-Diaz, L. Asquith, A. Aurisano, A. Back, C. Backhouse, et al., “Improved measurement of neutrino oscillation parameters by the nova experiment,” Physical Review D, vol. 106, no. 3, p. 032004, 2022.
  * (42) Y. Muramatsu and Y. Shigekami, “Susy threshold corrections to quark and lepton mixing inspired by so (10) gut models,” Journal of High Energy Physics, vol. 2020, no. 2, pp. 1–24, 2020.
  * (43) R. N. Mohapatra, M. Parida, and G. Rajasekaran, “Threshold effects on quasidegenerate neutrinos with high-scale mixing unification,” Physical Review D, vol. 71, no. 5, p. 057301, 2005.
  * (44) S. Gupta, S. K. Kang, and C. Kim, “Renormalization group evolution of neutrino parameters in presence of seesaw threshold effects and majorana phases,” Nuclear Physics B, vol. 893, pp. 89–106, 2015.
  * (45) W. G. Hollik, “Lifting degenerate neutrino masses, threshold corrections and maximal mixing,” arXiv preprint arXiv:1412.5117, 2014.
  * (46) P. A. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. Banday, R. Barreiro, J. Bartlett, N. Bartolo, et al., “Planck 2015 results-xiii. cosmological parameters,” Astronomy & Astrophysics, vol. 594, p. A13, 2016.
  * (47) P. Chen, S. C. Chuliá, G.-J. Ding, R. Srivastava, and J. W. Valle, “Cp symmetries as guiding posts: revamping tri-bi-maximal mixing. part i,” Journal of High Energy Physics, vol. 2019, no. 3, pp. 1–27, 2019.
  * (48) P. Chen, S. C. Chuliá, G.-J. Ding, R. Srivastava, and J. W. Valle, “C p symmetries as guiding posts: Revamping tribimaximal mixing. ii.,” Physical Review D, vol. 100, no. 5, p. 053001, 2019.
  * (49) P. Wilina, M. S. Singh, and N. N. Singh, “Deviations from tribimaximal and golden ratio mixings under radiative corrections of neutrino masses and mixings,” arXiv preprint arXiv:2205.01936, 2022.
  * (50) K. Devi, K. S. Singh, N. N. Singh, et al., “Effects of variations of susy breaking scale on neutrino parameters at low energy scale under radiative corrections,” Advances in High Energy Physics, vol. 2022, 2022\.
  * (51) R. Srivastava, “Predictions from high scale mixing unification hypothesis,” Pramana, vol. 86, no. 2, pp. 425–436, 2016.
  * (52) S. AbdusSalam, M. Abyaneh, F. Ghelichkhani, and M. Noormandipour, “Majorana phases in high-scale mixing unification hypotheses,” International Journal of Modern Physics A, vol. 36, no. 13, p. 2150077, 2021.
  * (53) M. Rajkhowa, M. K. Das, and N. N. Singh, “Possible radiative origin of $\delta$ 2 21 and sin$\theta$ 13 revisited,” Indian Journal of Physics, vol. 80, pp. 829–835, 2006.