# Non-Markovian SIR epidemic spreading model

Lasko Basnarkov1,2<EMAIL_ADDRESS>Igor Tomovski2 Trifce
Sandev2,3,4 Ljupco Kocarev1,2 1SS. Cyril and Methodius University, Faculty of
Computer Science and Engineering, P.O. Box 393, 1000 Skopje, Macedonia
2Macedonian Academy of Sciences and Arts, P.O. Box 428, 1000 Skopje, Macedonia
3Institute of Physics & Astronomy Karl-Liebknecht-Str. 24/25, University of
Potsdam, D-14476 Potsdam-Golm, Germany  4Institute of Physics, Faculty of
Natural Sciences and Mathematics, Ss Cyril and Methodius University,
Arhimedova 3, 1000 Skopje, Macedonia

###### Abstract

We introduce non-Markovian SIR epidemic spreading model inspired by the
characteristics of the COVID-19, by considering discrete- and continuous-time
versions. The incubation period, delayed infectiousness and the distribution
of the recovery period are modeled with general functions. By taking
corresponding choice of these functions, it is shown that the model reduces to
the classical Markovian case. The epidemic threshold is analytically
determined for arbitrary functions of infectivity and recovery and verified
numerically. The relevance of the model is shown by modeling the first wave of
the epidemic in Italy, in the spring, 2020.

††preprint: APS/123-QED

##  I Introduction

The ongoing pandemics of COVID-19, has claimed millions of human lives, caused
stagnation of the global economy and excessive load on the healthcare systems
throughout the world and changed the normal life. Mathematical models of
epidemic spreading are important tools for predicting the effects that the
pandemics can have on each segment of the society. They provide support for
policy-makers to make adequate decisions in order to partially mitigate the
consequences by planning various social distancing measures, preparation of
healthcare facilities and appropriate adaptation of the economy.

The spectrum of mathematical models applied for the COVID-19 pandemic ranges
from the simplest SIR to rather complex SIDARTHE Roda  _et al._ (2020); Zhao
and Chen (2020); Calafiore  _et al._ (2020); Giordano  _et al._ (2020); Gatto
_et al._ (2020), which are used for assessment of different aspects of the
epidemics. One of the major features of these models is their Markovian
nature, which considers transitions from one state to another to be
independent on the past. As an example, when Markovian property is assumed to
hold, an individual that has just become infected can proceed to recovered
state with the same probability as another one which has been infected for
longer period. This Markovian assumption, encapsulated in constant transition
probabilities, or rates, makes the models easier to study analytically. The
outcomes of these studies with Markovian approach offer some, and in certain
instances satisfactory, assessment of the spreading dynamics. However, growing
body of evidence, particularly for the COVID-19, suggests existence of
incubation period and certain infectivity patterns, with possibility for
spreading the pathogen before onset of the symptoms, to which correspond
functions that are rather distinct from the exponential distribution which the
Markovian models rely on Qin  _et al._ (2020a, b). Although adding one or more
compartments for the Exposed, Asymptomatic, Presymptomatic, or Quarantined
persons or considering various kinds of delay Liu  _et al._ (2020); Dell’Anna
(2020); Rong  _et al._ (2020) address such observations to certain extent,
they cannot systematically incorporate the observed distributions of the
incubation period and the healing process.

The non-Markovian setting is inherent in the pioneering works in the
mathematical epidemiology by Ross Ross (1916); Ross and Hudson (1917a, b),
Kermack and McKendrick Kermack and McKendrick (1927), and in the related field
of population dynamics by Bröck Böckh (1877) and Lotka Lotka (1919). However,
the more special and mathematically more tractable, Markovian approach has
largely dominated in subsequent studies. In the recent time the non-Markovian
framework has started to gain more attention in various settings. In one
attempt Boguñá  _et al._ (2014) is proposed Gillespie algorithm as an adequate
tool for numerical analysis of non-Markovian spreading models. The effects of
the form of distribution of infection and curing (recovery) times on SIS
epidemic model occurring on complex networks in continuous time has been
analyzed in several studies Starnini  _et al._ (2017); Van Mieghem and van de
Bovenkamp (2013); Van Mieghem and Liu (2019); Liu and Van Mieghem (2018); Feng
_et al._ (2019); Krylova and Earn (2013). With the introduction of SI*V* model
Nowzari  _et al._ (2015) it was suggested that non-Markovian spreading models
have capacity to be extended to cover a wide variety of spreading sub-models
and variants. Nontrivial distribution of infectious period in an integro-
differential SIR model was considered in Riaño (2020). In a recent study, non-
Markovian SIS model on complex networks, with arbitrary function for
infectivity and recovery was proposed Tomovski  _et al._ (2021), in which
control theory was successfully applied for determination of epidemic
threshold. Another, novel key contributions in the theory of non-Markovian
epidemic spreading models can be considered Pang and Pardoux (2020); Forien
_et al._ (2021). In those works, with extensive theoretical work on models
with integro-differential equations were obtained analytical results about the
equilibria and the basic reproduction numbers. Our study adds determination of
the epidemic threshold on base on the stability analysis for general
distributions of infectivity and healing in a SIR model. By similar approach
as in Tomovski  _et al._ (2021) we show how these functions determine the
epidemic threshold. The relevance of the model, besides by numerical
simulations, is verified by fitting to the observations of the first wave of
the epidemic in Italy, in the spring, 2020.

The paper is organized as follows. After providing initial setting of the
model in Section II, we introduce the discrete-time and continuous-time models
in Sections III and IV, respectively, where we also derive the epidemic
threshold relationships. The reduction to Markovian case of the model is
presented in Section V, while numerical simulations and discussions are given
in Section VI. The paper concludes with Section VII.

##  II Preliminaries

We consider SIR model that has three compartments: Susceptible - S, Infected -
I and Recovered - R, with the usual transition $S\to I\to R$. Let the
corresponding variables $S$, $I$ and $R$ denote the fractions of the
population that are in the given state, and under assumption without births
and deaths, one has the normalization condition $S(t)+I(t)+R(t)=1$ at each
moment $t$. To capture the nontrivial dependence of the healing period and the
different contagiousness of the infected individual in different stages of the
disease we introduce two functions. The infectivity function $\beta(\tau)$
captures the rate, or probability at which individuals that became infected
before time $\tau$ are spreading the disease to the susceptible ones. Thus, by
simply taking $\beta(\tau)=0$ for $\tau<T_{0}$, one is able to introduce
incubation period with length $T_{0}$. Another important function is the
healing function $\gamma(\tau)$ that denotes the probability with which
individual can heal at moment $\tau$ after contracting the disease. To account
for asymptomatic transmitters and existence of certain time window when
presence of pathogen can be confirmed, one can introduce a reporting function
$\rho(\tau)$. It is associated to the probability that the presence of the
pathogen can be confirmed at moment $\tau$ after contraction with it. The
asymptomatic cases are conveniently handled by normalizing the reporting
function to value smaller than unity. We pursue by considering discrete- and
continuous-time models separately, and provide more details about these
functions.

##  III Discrete-time version

In this section we consider evolution in discrete time $t$ and denote the
fraction of individuals that have become infected within the continuous-time
interval $[t-1,t]$ with $I_{d}(t)$, where for simplicity the unit interval is
taken to be 1. This can be relevant for situations like those when cases are
considered on daily basis. In such scenario, we have discrete-time healing
function $\gamma(\tau)$ and infectivity one $\beta(\tau)$, on which we put the
constraint $\beta(0)=0$. The probability that the individual will heal within
$\tau$ time units is $\Gamma(\tau)=\sum_{\nu=0}^{\tau}\gamma(\tau)$. We
further assume finite duration $T$ of the disease, what implies $\Gamma(T)=1$
and for practical reasons introduce its complement
$\overline{\Gamma}(\tau)=1-\Gamma(\tau)$, to denote the probability that
individual has not healed yet for $\tau$ time units. The function
$\gamma(\tau)$ also has a meaning of fraction of individuals that have
contracted the disease within the same unit time interval, to become healed
later within another unit interval $[\tau-1,\tau]$. Similar reasoning holds
for the cumulative functions $\Gamma(\tau)$ and $\overline{\Gamma}(\tau)$. On
base on the classical SIR model, the proposed model of evolution of the
compartments is given with the system

$\displaystyle S(t+1)$ $\displaystyle=$ $\displaystyle
S(t)\left[1-\sum_{\tau=0}^{T-1}\beta(\tau)\overline{\Gamma}(\tau)I_{d}(t-\tau)\right]$
$\displaystyle I_{d}(t+1)$ $\displaystyle=$ $\displaystyle
S(t)\sum_{\tau=0}^{T-1}\beta(\tau)\overline{\Gamma}(\tau)I_{d}(t-\tau)$
$\displaystyle R(t+1)$ $\displaystyle=$ $\displaystyle
R(t)+\sum_{\tau=0}^{T-1}\gamma(\tau)I_{d}(t+1-\tau).$ (1)

One can note that the infected individuals that have contracted the pathogen
up to $T$ periods before the current moment $t$, and which are not healed yet,
can contribute to spreading of the disease, with appropriate intensity
captured in the function $\beta(\tau)$. We note that in order to determine the
infected fraction at given moment, one should sum those infected in the past,
but did not heal up to the given moment

$I(t)=\sum_{\tau=0}^{T-1}I_{d}(t-\tau)\overline{\Gamma}(\tau).$ (2)

To make the problem completely defined one has to specify the initial
conditions for $I_{d}(t)$. We assume that they are given for
$\tau=T-1,T-2,\dots,0$. In general this model cannot be solved analytically
and should be studied by application of numerical simulations.

To get insight of the conditions when epidemic can emerge, one can determine
the stability of the disease free state $S^{*}=1,I^{*}=I_{d}^{*}=R^{*}=0$,
that is an equilibrium point of the system. Its local stability is established
by linearizing the dynamical equations (1) in its neighborhood. By making the
linearization in vicinity of $S^{*}=1,I^{*}=R^{*}=0$, one can observe the
dynamical evolution of the perturbations $\delta S=S-S^{*},\delta
I_{d}=I_{d}-I_{d}^{*},\delta R=R-R^{*}$. Under linearization, the
perturbations are related with

$\displaystyle\delta S(t+1)$ $\displaystyle=$ $\displaystyle\delta
S(t)-\sum_{\tau=0}^{T-1}\beta(\tau)\overline{\Gamma}(\tau)\delta
I_{d}(t-\tau),$ $\displaystyle\delta I_{d}(t+1)$ $\displaystyle=$
$\displaystyle\sum_{\tau=0}^{T-1}\beta(\tau)\overline{\Gamma}(\tau)\delta
I_{d}(t-\tau),$ $\displaystyle\delta R(t+1)$ $\displaystyle=$
$\displaystyle\delta R(t)+\sum_{\tau=0}^{T-1}\gamma(\tau)\delta
I_{d}(t+1-\tau).$ (3)

Let us focus on the infected fraction and make $Z$-transform on the second
equation in (3). To do so, multiply first both sides of that equation by
$z^{-t}$ and sum to obtain

$\sum_{t=0}^{\infty}\delta
I_{d}(t+1)z^{-t}=\sum_{t=0}^{\infty}\sum_{\tau=0}^{T-1}\beta(\tau)\overline{\Gamma}(\tau)\delta
I_{d}(t-\tau)z^{-t}.$ (4)

By using the $Z$-transform of the fraction of the population that become
infected at unit interval $I_{d}(t)$, given as
$\mathcal{I}(z)=\sum_{t=0}^{\infty}I_{d}(t)z^{-t}$, the left hand side of (4)
will become

$\displaystyle\sum_{t=0}^{\infty}\delta I_{d}(t+1)z^{-t}$ $\displaystyle=$
$\displaystyle z\sum_{t=0}^{\infty}\delta I_{d}(t+1)z^{-(t+1)}$ (5)
$\displaystyle=$ $\displaystyle z\left[\mathcal{I}(z)-\delta I_{d}(0)\right].$

Accordingly, the right-hand side of (4) can be rearranged as

$\displaystyle\sum_{t=0}^{\infty}\sum_{\tau=0}^{T-1}\beta(\tau)\overline{\Gamma}(\tau)\delta
I_{d}(t-\tau)z^{-t}=$ $\displaystyle=$
$\displaystyle\sum_{\tau=0}^{T-1}\beta(\tau)\overline{\Gamma}(\tau)z^{-\tau}\sum_{t=0}^{\infty}\delta
I_{d}(t-\tau)z^{-(t-\tau)}.$ (6)

By using substitution $\nu=t-\tau$, the last sum for $\tau\leq-1$ can be
expressed as

$\sum_{\nu=-\tau}^{\infty}\delta
I_{d}(\nu)z^{-\nu}=\sum_{\nu=-\tau}^{-1}\delta
I_{d}(\nu)z^{-\nu}+\mathcal{I}(z)=\mathcal{I}_{0}(\tau,z)+\mathcal{I}(z),$ (7)

where we have introduced a function $\mathcal{I}_{0}(\tau,z)$ that corresponds
to the initial conditions. Now, combining the relationships (5) – (7) one has

$z\left[\mathcal{I}(z)-\delta
I_{d}(0)\right]=\sum_{\tau=0}^{T-1}\beta(\tau)\overline{\Gamma}(\tau)\left[\mathcal{I}_{0}(\tau,z)+\mathcal{I}(z)\right]z^{-\tau}.$
(8)

To shorten the notation, one can introduce the following two complex functions

$\displaystyle\mathcal{E}(z)$ $\displaystyle=$
$\displaystyle\sum_{\tau=0}^{T-1}\beta(\tau)\overline{\Gamma}(\tau)z^{-\tau},$
$\displaystyle\mathcal{E}_{0}(z)$ $\displaystyle=$
$\displaystyle\sum_{\tau=0}^{T-1}\beta(\tau)\overline{\Gamma}(\tau)\mathcal{I}_{0}(\tau,z)z^{-\tau}.$
(9)

The first one is simply the $Z$-transform $\mathcal{E}(z)$ of what might be
called epidemic function $E(\tau)=\beta(\tau)\overline{\Gamma}(\tau)$, that is
a combination of the infecting and healing functions because
$\sum_{\tau=0}^{T-1}\beta(\tau)\overline{\Gamma}(\tau)z^{-\tau}=\sum_{\tau=0}^{\infty}\beta(\tau)\overline{\Gamma}(\tau)z^{-\tau}$.
The second complex function $\mathcal{E}_{0}(z)$ is related to the initial
conditions. Now, one has the following relationship

$z\left[\mathcal{I}(z)-\delta
I_{d}(0)\right]=\mathcal{I}(z)\mathcal{E}(z)+\mathcal{E}_{0}(z),$ (10)

from where

$\mathcal{I}(z)=\frac{z\delta I_{d}(0)+\mathcal{E}_{0}(z)}{z-\mathcal{E}(z)}.$
(11)

From a result in theory of discrete linear time-invariant systems, a sequence
(the impulse response of such system) is decaying if the poles of its
$Z$-transform are within the unit circle Oppenheim Alan and Willsky Alan
(2013). Thus, when the poles of the function $\mathcal{I}(z)$ of the complex
function (11), or the roots of the polynomial $z-\mathcal{E}(z)$ lie within
the unit circle, the perturbation dies out at infinity. So, the epidemic
threshold can be obtained by taking $z=1$ in the denominator in (11), that
results in

$\sum_{\tau=0}^{T-1}\beta(\tau)\overline{\Gamma}(\tau)=1,$ (12)

which obviously depends on the functional forms of the infectivity and healing
functions.

We should finally note that any initial infection would not shift back the
population to the disease-free state $S=1,I=R=0$, but to some endemic
$S_{e}^{*},I_{e}^{*}=0,R_{e}^{*}=1-S^{*}$. However, if the conditions are not
favoring epidemic both equilibria will be rather close $S_{e}^{*}\approx 1$.

##  IV Continuous-time version

We will pursue similarly to the discrete-time approach, where the fractions of
individuals within given compartment and the functions modeling the
infectivity, healing and reporting are defined for continuous time $t$ and we
use the same notation. Thus, $S(t)$ is the fraction of susceptible individuals
at given moment $t$ and $R(t)$ corresponds to the recovered and again assume
finite healing period $T$. The fraction of infected individuals is
conveniently modeled with the rate of infection, or the fraction of newly
infected individuals $I_{d}(t)$ within the infinitesimal interval $(t-dt,t)$.
The total fraction of infected persons is given with the integral

$I(t)=\int_{0}^{T}I_{d}(t-\tau)\overline{\Gamma}(\tau)d\tau,$ (13)

which accounts for those that had become infected in the past and have not
healed yet. Now, the dynamical evolution of the respective fractions is given
with

$\displaystyle\dot{S}$ $\displaystyle=$
$\displaystyle-S(t)\int_{0}^{T}\beta(\tau)\overline{\Gamma}(\tau)I_{d}(t-\tau)d\tau$
$\displaystyle I_{d}(t)$ $\displaystyle=$ $\displaystyle
S(t)\int_{0}^{T}\beta(\tau)\overline{\Gamma}(\tau)I_{d}(t-\tau)d\tau$
$\displaystyle\dot{R}$ $\displaystyle=$
$\displaystyle\int_{0}^{T}\gamma(\tau)I_{d}(t-\tau)d\tau.$ (14)

One should note that in their original approach, the general version of the
model by Kermack and McKendrick assumes dependence of the infectivity on the
age of infection just as the last relationships (14) suggests Kermack and
McKendrick (1927); Brauer (2017). In order to determine whether the initial
perturbation will grow to epidemics, one could focus on the second equation in
the vicinity of the disease-free state $S^{*}=1,R^{*}=I^{*}=0$. Then, the
perturbation of newly infected individuals will evolve as

$\delta I_{d}(t)=\int_{0}^{T}\beta(\tau)\overline{\Gamma}(\tau)\delta
I_{d}(t-\tau)d\tau,$ (15)

where it is assumed that in vicinity of the disease-free state $S(t)\approx
1$. Now, make Laplace transform of the perturbation of the rate of infection,
$\mathcal{I}(s)=\int_{0}^{\infty}\delta I_{d}(t)e^{-st}dt$ and use it in the
last equation (15). To do that, we will follow the same approach as in the
discrete-time version. Multiply both sides with $e^{-st}$ and integrate. The
left hand side will result in the Laplace transform of $\delta I_{d}(t)$,
while the right hand one will be

$\displaystyle A$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\int_{0}^{T}\beta(\tau)\overline{\Gamma}(\tau)\delta
I_{d}(t-\tau)e^{-st}d\tau dt$ (16) $\displaystyle=$
$\displaystyle\int_{0}^{T}\beta(\tau)\overline{\Gamma}(\tau)e^{-s\tau}\int_{0}^{\infty}\delta
I_{d}(t-\tau)e^{-s(t-\tau)}dt$ $\displaystyle=$
$\displaystyle\int_{0}^{T}\beta(\tau)\overline{\Gamma}(\tau)e^{-s\tau}\int_{-\tau}^{\infty}\delta
I_{d}(\nu)e^{-s\nu}d\nu$

The last integral can be expressed with

$\displaystyle\int_{-\tau}^{\infty}I_{d}(\nu)e^{-s\nu}d\nu$ $\displaystyle=$
$\displaystyle\int_{-\tau}^{0}I_{d}(\nu)e^{-s\nu}d\nu+\mathcal{I}(s)$ (17)
$\displaystyle=$ $\displaystyle\mathcal{I}_{0}(\tau,s)+\mathcal{I}(s).$

Now, one has

$A=\int_{0}^{T}\beta(\tau)\overline{\Gamma}(\tau)e^{-s\tau}\left[\mathcal{I}_{0}(\tau,s)+\mathcal{I}(s)\right]d\tau.$
(18)

Similarly to the discrete-time case we can introduce the Laplace transform of
the epidemic function $E(\tau)=\beta(\tau)\overline{\Gamma}(\tau)$ and its
initial conditions contribution

$\displaystyle\mathcal{E}(s)$ $\displaystyle=$
$\displaystyle\int_{0}^{T}\beta(\tau)\overline{\Gamma}(\tau)e^{-s\tau}d\tau,$
(19) $\displaystyle\mathcal{E}_{0}(s)$ $\displaystyle=$
$\displaystyle\int_{0}^{T}\beta(\tau)\overline{\Gamma}(\tau)\mathcal{I}_{0}(\tau,s)e^{-s\tau}d\tau.$

Finally, one obtains

$\mathcal{I}(s)=\mathcal{I}(s)\mathcal{E}(s)+\mathcal{E}_{0}(s),$ (20)

from where the Laplace transform of the perturbation of the infection rate is

$\mathcal{I}(s)=\frac{\mathcal{E}_{0}(s)}{1-\mathcal{E}(s)}.$ (21)

From the results of control theory, a continuous-time linear time-invariant
system is stable if the poles of its transfer function, or Laplace transform
of its impulse response have negative real part Oppenheim Alan and Willsky
Alan (2013). Thus, the perturbations $\delta I_{d}(t)$ will decay if the poles
of its Laplace transform $\mathcal{I}(s)$ (21), or eigenvalues of the system
(14) lie within negative half-plane $Re\\{s\\}<0$. Then, the epidemic
threshold can be obtained with $s=0$ which leads to

$\int_{0}^{T}\beta(\tau)\overline{\Gamma}(\tau)d\tau=1,$ (22)

that represents the relationship, which corresponds to the discrete-time case
(12).

##  V Markovian SIR model

In order to obtain the classical Markovian SIR model, from the non-Markovian
case (1), one should consider taking $T\rightarrow\infty$,
$\beta(\tau)=\beta$, $\gamma(0)=0$, $\gamma(\tau)=\gamma(1-\gamma)^{\tau-1}$
where $\beta$ and $\gamma$ are constants. This further yields $\Gamma(0)=0$,
$\Gamma(\tau)=1-(1-\gamma)^{\tau}$ and
$\overline{\Gamma}(\tau)=(1-\gamma)^{\tau}=\gamma(\tau+1)/\gamma$. First, one
could observe that by using constant infectivity $\beta(\tau)=\beta$ in the
first relationship of the model (1) and using (2) one will obtain the
classical form for evolution of the susceptible population

$S(t+1)=S(t)\left[1-\beta I(t)\right].$ (23)

Next, by implementing the condition $\gamma(0)=0$, and the relationship
$\overline{\Gamma}(\tau)=\gamma(\tau+1)/\gamma$ one can drop the first term in
the sum in the recovered population in (1), and further obtain

$\displaystyle\sum_{\tau=0}^{T-2}\gamma(\tau+1)I_{d}(t-\tau)$
$\displaystyle=\gamma\sum_{\tau=0}^{T-1}\overline{\Gamma}(\tau)I_{d}(t-\tau)-\gamma\overline{\Gamma}(T-1)I_{d}(t-T+1)$
$\displaystyle=\gamma I(t)-\gamma(1-\gamma)^{T-1}I_{d}(t-T+1),$ (24)

from where, for $T\rightarrow\infty$, the recovered population evolves as

$R(t+1)=R(t)+\gamma I(t).$ (25)

Finally, from the conservation relationship $I(t)+S(t)+R(t)=1$, one can find
that the infected fraction is given as

$I(t+1)=\beta S(t)I(t)+(1-\gamma)I(t).$ (26)

The relationships (23), (25) and (26) represent the classical SIR model in
discrete time.

As an example, in the figure 1 we make a comparison between numerical
solutions of the discrete classical SIR model and the classical SIR -
equivalent model obtained from the non-Markovian form.

Figure 1: Comparison between the discrete classical SIR model and the
classical SIR - equivalent model obtained from the non-Markovian form, for
$\beta=0.2$, $\gamma=0.03$. It is used rather large finite duration of the
healing $T=150$, as a proxy for $T\to\infty$.

Similarly to the discrete-time version, to verify that the proposed continuous
model is generalization of the classical, Markovian SIR model, one should
consider two characteristics of the latter: 1. The infection rate is
independent on the moment when the disease was contracted $\beta(\tau)=\beta$;
and 2. The duration of infectivity is infinite and exponentially distributed
which implies that the healing function is $\gamma(\tau)=\lambda
e^{-\lambda\tau}$. We note that the respective cumulative distribution is
$\Gamma(\tau)=1-e^{-\lambda\tau}$, and accordingly
$\overline{\Gamma}(\tau)=e^{-\lambda\tau}$. By using the functional form of
the healing function, the total infectious population will be

$I(t)=\int_{0}^{\infty}\overline{\Gamma}(\tau)I_{d}(t-\tau)d\tau=\int_{0}^{\infty}e^{-\lambda\tau}I_{d}(t-\tau)d\tau.$
(27)

Similarly, by using $\beta(\tau)=\beta$, for the dynamics of the susceptible
fraction one has

$\dot{S}=-S(t)\beta\int_{0}^{\infty}e^{-\lambda\tau}I_{d}(t-\tau)d\tau=-\beta
SI,$ (28)

that represents the corresponding relationship in the classical SIR model.
Furthermore, by applying the functional form for the healing function, the
dynamics of the recovered population will be as follows

$\displaystyle\dot{R}=\int_{0}^{\infty}\lambda
e^{-\lambda\tau}I_{d}(t-\tau)d\tau=\lambda I,$ (29)

that is the respective relationship in the classical SIR model. Finally, by
using the conservation principle $S(t)+I(t)+R(t)=1$, the total infectious
fraction will evolve as

$\dot{I}=-\dot{S}-\dot{R}=\beta SI-\lambda I,$ (30)

that is the remaining familiar relationship from the classical case. As a
final note, we just mention that using respective forms for the infectivity
and recovery functions for the Markovian case in the epidemic threshold
relationships (12) and (22), one will obtain the familiar threshold
$\beta_{\mathrm{t}h}=\gamma$.

##  VI Numerical experiments and discussion

Our numerical experiments with the proposed model were based on solution of
the integro-differential equations for the continuous-time case. We have used
the Euler method with step $\Delta t=0.01$. Although in the model can be used
arbitrary functions of infection and recovery, we have chosen to use those
that can been found in the literature as appropriate for the COVID-19
pandemic. As suggested in Qin  _et al._ (2020a) the infectivity function
$\beta(\tau)$ is conveniently represented with Weibull probability density
function, with parameters $\alpha=2.04$ and $\lambda=0.103$, which is further
truncated to 35 days and normalized. The daily recovering probabilities were
modeled with log-normal probability density function
$L(\tau;\mu;\sigma)=1/(\tau\sigma\sqrt{2\pi})\exp(-\left(\ln\tau-\mu\right)^{2}/(2\sigma^{2}))$,
with parameters $\mu=\ln(\mu_{X}^{2}/(\sqrt{\mu_{X}^{2}+\sigma_{X}^{2}})$,
$\sigma^{2}=\ln\left(1+\sigma_{X}^{2}/\mu_{X}^{2}\right)$ chosen to match a
mean value of $\mu_{X}=21$ and standard deviation $\sigma_{X}=6$. The
distribution is then normalized to 61 days, and time-shifted for 4 days in
order to exclude immediate recovery. This results in the healing function
$\gamma(\tau)$ with mean recovery time of $25\pm 6$ days, in the following
fashion

$\displaystyle\gamma(\tau)=\begin{cases}\frac{L(\tau-4;\mu;\sigma)}{\int_{0}^{61}L(\tau;\mu;\sigma)d\tau},&4\leq\tau\leq
65,\\\ 0,&\text{otherwise}.\end{cases}$ (31)

This construct was based on the results from Sreevalsan-Nair  _et al._ (2020);
Faes  _et al._ (2020), assuming that: 1. Onset of symptoms (on average) occurs
after four days (the time shift); 2. It takes another 7-10 days from onset of
symptoms to diagnosis confirmation and hospitalization; 3. Another 10-11 days,
on average, are needed from hospitalization to recovery. The period of $T=65$
days is considered in order to include even most extreme cases in which
hospitalization exceeded 40 days.

Furthermore, we have chosen to scale the infectivity function with a parameter
$\beta$ given in terms of the epidemic threshold $\beta_{\mathrm{th}}$. The
threshold value was obtained from the condition (22)

$\beta_{\mathrm{th}}\int_{0}^{T}\beta(\tau)\overline{\Gamma}(\tau)d\tau=1.$
(32)

To verify the value of the epidemic threshold we have varied the infectivity
parameter in vicinity of the critical value obtained from (32) and run the
continuous-time model for total time equal to 5000. The final values of the
susceptible and recovered fraction are plotted as function of the infectivity
parameter in the figure 2. As one can see, once $\beta$ is larger than its
critical value, the epidemic emerges.

Figure 2: Fractions of susceptible (red stars) and recovered (blue dots)
individuals at the end of the epidemic as a function of the scaling of the
infectivity function $\beta$ given in terms of its threshold value
$\beta_{\mathrm{th}}$.

In order to verify how well the approach can be used to model the COVID-19
pandemic we have chosen to use value of the infectivity parameter $\beta$ that
nearly matches the growth patterns of the epidemic in the countries before
countermeasures were applied. As was obtained in a detailed study Pellis  _et
al._ (2020), the epidemic doubling time in many countries is approximately
three days. For that reason, we have opted to use the value
$\beta=4.85\beta_{\mathrm{th}}$ that produces such growth. We have numerically
verified that in the initial stage of the epidemic, the newly confirmed daily
cases and the total number of infected individuals grow with the same rate,
and have the same doubling time. Also, by running the model with
$\beta=4.85\beta_{\mathrm{th}}$ for very long time, we have obtained that at
the end less than 1% of the population will remain susceptible! This result
means that, if the doubling time is three days in case of free spreading of
the virus, then prevention of the epidemic would need nearly everyone should
be either vaccinated or had healed from the virus. This is particular
challenge of the model that should be addressed carefully.

We have finally attempted to check how well the model can explain the
observations. To do so, we have used the COVID-19 data from Our World in Data,
for Italy. Our focus was put on the first wave of the pandemic, since in its
beginning no preventive measures were used. We have chosen to study the
epidemic in Italy, where the wave was the strongest. The window of data under
study starts from February 21, 2020, that is the date from which every day
were reported new cases. The countrywide lockdown started on March 10, 2020,
that corresponds to day 19 in this study. We have used value of infectivity
$\beta\approx 3.7\beta_{\mathrm{th}}$ that provided good fit to the observed
data for the period from February 21 until March 9. This value was used until
the start of the lockdown, when it was set to certain value smaller than the
threshold. The initial condition was set to $I_{d}(0)=10^{-7}$, that for Italy
would mean about 6 persons infected at the starting day of simulation. We have
chosen to apply detection of the infected individuals on based on a function
that has identical form as the infectivity one, but which is delayed for
certain number of days. This corresponds to situation that only those with
symptoms are tested, and their appearance is delayed few days after the onset
of infectivity. Also, there is certain delay that corresponds to the whole
process from onset of symptoms, to visit to hospital to obtaining positive
result. We note that the testing function was normalized to 0.8 that
corresponds to assuming existence of 20% asymptomatic cases Buitrago-Garcia
_et al._ (2020). To reach a good fit to the observations we had to take the
start of the simulation, that is the day when the initial seed of infection
was set, to be approximately 60 days before the day 1, when comparison with
the real data starts. Its exact value was obtained by fitting the logarithms
of the daily detected cases from the simulation to the respective ones from
the data. More precisely, we have looked for a shift $s$, that will result in
minimal squared error as follows

$\epsilon=\arg\min_{s}\left\\{\frac{1}{19}\sum_{k=1}^{19}\left[\ln(I_{d}^{\mathrm{data}}(k))-\ln(I_{d}(k+s))\right]\right\\}.$
(33)

We report in the top panel of figure 3 two simulated scenarios compared to the
observations. In the first case we took testing function that is delayed after
the infectivity one for two days, that actually becomes nonzero at the
possible onset of the symptoms He  _et al._ (2020), while the other case
corresponds to delay of five days. The latter scenario provides much better
fit to the observations, particularly in the period after the lockdown has
started, and even further in the period after the peak, as one can notice in
the figure 3. We have tried with all integer values of the delay from two to
ten (not shown) and five days correspond to the best fit. We remind that the
lockdown corresponds to day 19 in the plot, while the peaks of the daily
reported cases are delayed: at day 27 and day 30, for scenario one, and two,
respectively. The peak at the latter case, appears at March 21, the day when
largest number of new cases were registered. This fit to the peak and beyond
of the model simulation with the observation, makes a good basis for the
relevance of the proposed framework. In the bottom panel in figure 3 we show
how modification of the value of infectivity parameter $\beta$ during the
lockdown phase influences the daily cases.

Figure 3: Daily confirmed cases in the first epidemic wave in Italy in spring
2020 (in blue squares), compared to numerical simulations of the model. Top
panel: Confirmation function is delayed for two days after onset of
infectivity (green circles) and five days (red stars) and
$\beta=0.75\beta_{\mathrm{th}}$. Bottom panel: Confirmation function is
delayed for five days, while the infectivity parameter is:
$\beta=0.5\beta_{\mathrm{th}}$ (green circles),
$\beta=0.75\beta_{\mathrm{th}}$ (red stars), and $\beta=\beta_{\mathrm{th}}$
(magenta crosses)

Although providing natural framework for incorporation of observed
distributions of the infectiousness of the infected individuals and the
typical development of the disease, the proposed model has drawbacks as well.
First, before using it, one needs to specify the functions modeling the
infectiousness, healing and discovering the infected individuals. Their
determination is a serious issue by their own and needs careful study. As more
complex one, the tuning of the model would need in general more data than the
classical Markovian counterparts. Also, its full specification needs providing
initial conditions that represent a high-dimensional vector, or an interval of
values. How all these factors shape the outcome of the model, and how much is
it robust to perturbations of any kind is unknown. We believe that their
understanding could provide the epidemiologists with valuable information for
better understanding of the possible outcomes of epidemics with pronounced
non-Markovian nature.

##  VII Conclusions

The proposed general non-Markovian epidemic spreading model captures the
typical patterns of the disease in person infected with SARS-CoV-2: delayed
onset of symptoms and infectivity and impossibility of immediate cure of those
that will become sick. We have studied both discrete- and continuous-time
versions and derived analytically the relationships for determination of the
epidemic threshold. The model reduces to the classical SIR model with the
corresponding choice of the functions of infection and healing. The
theoretical analysis was supported by numerical confirmation of the epidemic
threshold values. The good fit of the model to the real data shows its
promising potential for application for modeling the spread of other
infectious diseases. By introducing other appropriate functions one could
possibly generalize this model to versions that include other compartments
that correspond to hospitalized, quarantined, or deceased persons.

Although the epidemic threshold as key quantity was determined, we did not
calculated the basic reproduction number $R_{0}$, that represents another
important quantity. Furthermore, the relationship between the scaling of the
infectivity function $\beta/\beta_{\mathrm{th}}$ from one side and $R_{0}$ and
the doubling time, from another should be explored as well. With this regard,
we think that it is even more important to determine the herd immunity level
needed to prevent the epidemic. Finally, analysis of epidemic spreading by
nontrivial contact patterns, modeled with complex networks, and by
incorporating the proposed approach could provide further insight in the
evolution of the epidemics. These issues could provide better understanding of
the non-Markovian setting in modeling the epidemic spreading.

##  VIII Acknowledgement

This research was partially supported by the Faculty of Computer Science and
Engineering, at the Ss. Cyril and Methodius University in Skopje, Macedonia.
The Authors acknowledge support by the German Science Foundation (DFG, Grant
number ME 1535/12-1).

## References

  * Roda  _et al._ (2020) W. C. Roda, M. B. Varughese, D. Han, and M. Y. Li, Infect. Dis. Model. 5, 271 (2020).
  * Zhao and Chen (2020) S. Zhao and H. Chen, Quant. Biol. , 1 (2020).
  * Calafiore  _et al._ (2020) G. C. Calafiore, C. Novara, and C. Possieri, Annu. Rev. Control (2020).
  * Giordano  _et al._ (2020) G. Giordano, F. Blanchini, R. Bruno, P. Colaneri, A. Di Filippo, A. Di Matteo, and M. Colaneri, Nat. Med. 26, 855 (2020).
  * Gatto  _et al._ (2020) M. Gatto, E. Bertuzzo, L. Mari, S. Miccoli, L. Carraro, R. Casagrandi, and A. Rinaldo, Proceedings of the National Academy of Sciences 117, 10484 (2020).
  * Qin  _et al._ (2020a) J. Qin, C. You, Q. Lin, T. Hu, S. Yu, and X.-H. Zhou, medRxiv , 2020.03.06.20032417 (2020a).
  * Qin  _et al._ (2020b) J. Qin, C. You, Q. Lin, T. Hu, S. Yu, and X.-H. Zhou, Sci. Adv. 6, eabc1202 (2020b).
  * Liu  _et al._ (2020) Z. Liu, P. Magal, O. Seydi, and G. Webb, Infectious Disease Modelling 5, 323 (2020).
  * Dell’Anna (2020) L. Dell’Anna, Sci. Rep. 10, 1 (2020).
  * Rong  _et al._ (2020) X. Rong, L. Yang, H. Chu, and M. Fan, Math. Biosci. Eng. 17, 2725 (2020).
  * Ross (1916) R. Ross, Proceedings of the Royal Society of London. Series A, Containing papers of a mathematical and physical character 92, 204 (1916).
  * Ross and Hudson (1917a) R. Ross and H. P. Hudson, Proceedings of the Royal Society of London. Series A, Containing papers of a mathematical and physical character 93, 212 (1917a).
  * Ross and Hudson (1917b) R. Ross and H. P. Hudson, Proceedings of the Royal Society of London. Series A, Containing papers of a mathematical and physical character 93, 225 (1917b).
  * Kermack and McKendrick (1927) W. O. Kermack and A. G. McKendrick, Proc. R. Soc. Lond. Series A, Mathematical and Physical Sciences 115, 700 (1927).
  * Böckh (1877) R. Böckh,  _Statistisches Jahrbuch der Stadt Berlin_ (Leonh. Simion, 1877).
  * Lotka (1919) A. J. Lotka, Journal of the Washington Academy of Sciences 9, 73 (1919).
  * Boguñá  _et al._ (2014) M. Boguñá, L. F. Lafuerza, R. Toral, and M. A. Serrano, Phys. Rev. E 90, 042108 (2014).
  * Starnini  _et al._ (2017) M. Starnini, J. P. Gleeson, and M. Boguñá, Phys. Rev. Lett. 118, 128301 (2017).
  * Van Mieghem and van de Bovenkamp (2013) P. Van Mieghem and R. van de Bovenkamp, Phys. Rev. Lett. 110, 108701 (2013).
  * Van Mieghem and Liu (2019) P. Van Mieghem and Q. Liu, Phys. Rev. E 100, 022317 (2019).
  * Liu and Van Mieghem (2018) Q. Liu and P. Van Mieghem, Phys. Rev. E 97, 022309 (2018).
  * Feng  _et al._ (2019) M. Feng, S.-M. Cai, M. Tang, and Y.-C. Lai, Nat. Commun. 10, 3748 (2019).
  * Krylova and Earn (2013) O. Krylova and D. J. Earn, J. R. Soc. Interface 10, 20130098 (2013).
  * Nowzari  _et al._ (2015) C. Nowzari, M. Ogura, V. M. Preciado, and G. J. Pappas, in  _2015 54th IEEE Conference on Decision and Control (CDC)_ (2015) pp. 5073–5078.
  * Riaño (2020) G. Riaño, medRxiv (2020), https://doi.org/10.1101/2020.05.15.20103465.
  * Tomovski  _et al._ (2021) I. Tomovski, L. Basnarkov, and A. Abazi, IEEE Trans. Netw. Sci. Eng. (2021), DOI: 10.1109/TNSE.2021.3125191.
  * Pang and Pardoux (2020) G. Pang and É. Pardoux, arXiv preprint arXiv:2003.03249 (2020).
  * Forien  _et al._ (2021) R. Forien, G. Pang, and É. Pardoux, Royal Society open science 8, 202327 (2021).
  * Oppenheim Alan and Willsky Alan (2013) V. Oppenheim Alan and S. Willsky Alan, “Signals and systems,” (2013).
  * Brauer (2017) F. Brauer, Infectious Disease Modelling 2, 113 (2017).
  * Sreevalsan-Nair  _et al._ (2020) J. Sreevalsan-Nair, R. R. Vangimalla, and P. R. Ghogale, medRxiv , 2020.04.17.20069724 (2020).
  * Faes  _et al._ (2020) C. Faes, S. Abrams, D. Van Beckhoven, G. Meyfroidt, E. Vlieghe, N. Hens, _et al._ , Int. J. Environ. Res. Public Health 17, 7560 (2020).
  * Pellis  _et al._ (2020) L. Pellis, F. Scarabel, H. B. Stage, C. E. Overton, L. H. Chappell, K. A. Lythgoe, E. Fearon, E. Bennett, J. Curran-Sebastian, R. Das, _et al._ , arXiv preprint arXiv:2004.00117 (2020).
  * Buitrago-Garcia  _et al._ (2020) D. Buitrago-Garcia, D. Egli-Gany, M. J. Counotte, S. Hossmann, H. Imeri, A. M. Ipekci, G. Salanti, and N. Low, PLoS Med. 17, e1003346 (2020).
  * He  _et al._ (2020) X. He, E. H. Lau, P. Wu, X. Deng, J. Wang, X. Hao, Y. C. Lau, J. Y. Wong, Y. Guan, X. Tan, _et al._ , Nat. Med. 26, 672 (2020).


# Investigating Two-zero Textures of Inverse Neutrino Mass Matrix under the
Lamp Post of LMA and LMA-D Solutions and Symmetry Realization

Labh Singh, Monal Kashav and Surender Verma
sainilabh5@gmail.commonalkashav<EMAIL_ADDRESS>Corresponding
Author

(Department of Physics and Astronomical Science
Central University of Himachal Pradesh
Dharamshala, India 176215 )

###### Abstract

In this work we have investigated the phenomenological consequences of two-
zero textures of inverse neutrino mass matrix ($M_{\nu}^{-1}$) in light of the
large mixing angle (LMA) and large mixing angle-dark (LMA-D) solutions, later
of which originates if neutrinos exhibit non-standard interactions with
matter. Out of fifteen possibilities, only seven two-zero textures of
$M_{\nu}^{-1}$ are found to be phenomenologically allowed under LMA and/or
LMA-D descriptions. In particular, five textures are in consonance with both
LMA and LMA-D solutions and are necessarily $CP$ violating while remaining two
textures are found to be consistent with LMA solution only. The textures with
vanishing (1, 1) element of $M_{\nu}^{-1}$ are, in general, disallowed. All
the textures allowed under LMA and LMA-D solutions follow the same neutrino
mass hierarchy. Furthermore, textures with vanishing (2, 3) element of
$M_{\nu}^{-1}$ are found to be either disallowed or are consistent with LMA
description only. We have, also, obtained the implication of the model for
$0\nu\beta\beta$ decay amplitude $|M_{ee}|$. For most of the textures the
calculated $3\sigma$ lower bound on $|M_{ee}|$ is $\mathcal{O}(10^{-2})$,
which is within the sensitivity reach of $0\nu\beta\beta$ decay experiments.
We have, also, proposed a flavor model based on discrete non-Abelian flavor
group $A_{4}$ wherein such textures of $M_{\nu}^{-1}$ can be realized within
Type-I seesaw setting.

Keywords: Neutrino mass matrix; Phenomenology; Majorana Neutrino; Texture
zeros; $A_{4}$ symmetry.
PACS Nos.: 14.60.Pq

##  1 Introduction

The experimental evidences accumulated in last two decades have convincingly
established that not only neutrino have non-zero mass but its dynamical origin
is beyond our current understanding of the standard model (SM). Despite the
best efforts to precisely decipher the structure of neutrino mixing matrix, it
still has certain unknowns, for example, $CP$ violating phases, octant of
$\theta_{23}$ and neutrino mass hierarchy, to name a few. The theoretical
frameworks to explain neutrino mass and its manifestations like neutrino
oscillations have been developed assuming standard (charged (CC) and neutral
current (NC)) interactions between neutrino and matter. These frameworks
culminated in large mixing angle (LMA) solution to the solar neutrino problem
(SNP) and is well established by the neutrino oscillation experiments. In fact
LMA has been independently confirmed as the solution to solar neutrino problem
(SNP) in solar and KamLAND reactor experiments[1]. It has been shown that for
positive solar mass-squared difference, $\sin^{2}\theta_{12}$ cannot be
greater than 0.5. However, at the subdominant level, there still have
possibilities for additional contributions to neutrino oscillations such as
non-standard interactions (NSI) of neutrino with matter fields[2, 3, 4]. The
future hyper-technological experiments will have the access to these
unexplored regions. In presence of NSI, $\theta_{12}$ can be in the second
octant. Thus, solar neutrino problem may have another degenerate solution in
which $\sin^{2}\theta_{12}\approx 0.7$. This solution is termed as large
mixing angle-dark (LMA-D) solution. In general, the degeneracy between LMA and
LMA-D solutions cannot be alleviated by oscillation experiments due to
generalized mass hierarchy degeneracy in presence of NSI. However, a combined
measurements from neutrino oscillation and scattering experiments may have
imperative implication with regard to lifting of these degeneracies[5, 6, 7].
Although non-standard interactions of neutrino with matter are severely
constrained, the latest global fit, incorporating neutrino oscillation and
COHERENT data, still allows LMA-D solution at $3\sigma$ confidence level[8].
The only difference in LMA and LMA-D is in the octant of $\theta_{12}$,
however, the solar mass-squared difference remains the same.

The progress in understanding the origin of neutrino mass centrally involves
explaining the observed pattern of neutrino mixing which is encoded in the
neutrino mass matrix obtained after electroweak symmetry breaking. Assuming
neutrino to be Majorana particle the mass matrix contain nine free parameters
viz. three neutrino mass eigenvalues, three mixing angles and three CP
violating phases. The two mass-squared differences and three mixing angles
have been measured in neutrino oscillation experiments with high degree of
precision. Seesaw mechanism is a natural and most effective way to explain the
smallness of neutrino mass. Within Type-I seesaw[9, 10] paradigm, the low
energy effective neutrino mass matrix ($M_{\nu}$) is generated from Dirac
neutrino mass matrix ($M_{D}$) and heavy right-handed Majorana neutrino mass
matrix ($M_{R}$) using the relation: $M_{\nu}\approx
M_{D}M_{R}^{-1}M_{D}^{T}$. The existence of near degeneracy in LMA and LMA-D
solutions have imperative implications for models of neutrino mass and
associated phenomenology. Recently, the LMA and LMA-D phenomenology of
Majorana neutrino mass matrix has been studied assuming (i) zero textures of
the neutrino mass matrix, $M_{\nu}$ [11, 12, 13] (ii) in presence of one-
sterile neutrino[14]. Texture-zeros in the effective low energy Majorana
neutrino mass matrix may have seesaw origin in which they can be realized from
the zeros in $M_{D}$ and $M_{R}$[15]. In literature, there have been
phenomenological studies with texture one-zero[16, 17, 18, 19] and two-zeros
[20, 21, 22, 23, 24] while three-zeros or more, in neutrino mass matrix, are
ruled out by current neutrino oscillation data[25]. In Type-I seesaw, an
interesting scenario may emerge if we work in $M_{D}$-diagonal basis. In this
basis, the zero(s) in $M_{\nu}^{-1}$ is same as the the zero(s) in $M_{R}$
i.e. $M_{\nu}^{-1}\approx M_{D}^{-1}M_{R}M_{D}^{-1}$ [26]. In Ref. [27], the
authors have investigated the phenomenological consequences of one-zero
texture in $M_{\nu}^{-1}$ with in the context of trimaximal mixing. The LMA
phenomenology of two-zero textures of $M_{\nu}^{-1}$ has been investigated in
[28]. It is to be noted that the texture zeros in $M_{\nu}$ and $M_{\nu}^{-1}$
are, in general, independent and may have distinguishing phenomenology.
Motivated by the capabilities of the future neutrino oscillation experiments
in resolving these subdominant effects[29] and its important model building
perspective, we investigate the phenomenological consequences of two-zero
textures of $M_{\nu}^{-1}$. Also, we have proposed a flavor model based on
non-Abelian discrete group $A_{4}$ and Type-I seesaw, where such zeros can be
realized in the inverse neutrino mass matrix. There are fifteen possibilities
to have two zeros in $M_{\nu}^{-1}$. We investigate the LMA and LMA-D
phenomenology of these fifteen textures.

The paper is organized as follows. In Sec. II, we briefly introduce the
formalism of two-zero texture of $M_{\nu}^{-1}$ and the details of the
numerical analysis. Sec. III is devoted to the investigation and discussion of
LMA/LMA-D phenomenology of the seven allowed textures of $M_{\nu}^{-1}$. A
flavor model based on non-Abelian $A_{4}$ symmetry is discussed in Sec. IV.
Finally, we brief our conclusions in Sec. V.

##  2 Formalism of Two-zero Textures of Inverse Neutrino Mass matrix

In charged lepton basis, the neutrino mass matrix is given by

$M_{\nu}=VM^{diag}_{\nu}V^{T},$ (1)

where $M^{diag}_{\nu}$= $diag(m_{1},m_{2},m_{3})$ is the neutrino mass
eigenvalue matrix, $V=U.P$ is complex unitary neutrino mixing matrix. The
matrix $U$ is Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix which in the PDG
representation can be written as [30]

$U=\begin{pmatrix}U_{11}&U_{12}&U_{13}\\\ U_{21}&U_{22}&U_{23}\\\
U_{31}&U_{32}&U_{33}\\\
\end{pmatrix}=\begin{pmatrix}c_{12}c_{13}&s_{12}c_{13}&s_{13}e^{-i\delta}\\\
-s_{12}c_{23}-c_{12}s_{23}s_{13}e^{i\delta}&c_{12}c_{23}-s_{12}s_{23}s_{13}e^{i\delta}&s_{23}c_{13}\\\
s_{12}s_{23}-c_{12}c_{23}s_{13}e^{i\delta}&-c_{12}s_{23}-s_{12}c_{23}s_{13}e^{i\delta}&c_{23}c_{13}\end{pmatrix},$
(2)

where $c_{ij}=cos\theta_{ij}$ and $s_{ij}=sin\theta_{ij}$ and $\delta$ is
Dirac-type $CP$ violating phase. Also, $P$ is the phase matrix given by

$P=\begin{pmatrix}1&0&0\\\ 0&e^{i\alpha}&0\\\ 0&0&e^{i(\beta+\delta)}\\\
\end{pmatrix},$

where $\alpha$, $\beta$ are Majorana-type $CP$ violating phases. The inverse
neutrino mass matrix can be derived from Eqn. (1) as

$M^{-1}_{\nu}=(VM^{diag}_{\nu}V^{T})^{-1}.$ (3)

Using Eqns. (1) and (2), six independent elements of ($M^{-1}_{\nu}$) can be
written as

$\left.\begin{aligned}
(M^{-1}_{\nu})_{11}&=\dfrac{1}{m_{1}m_{2}m_{3}}\left[c^{2}_{13}e^{-2i\alpha}m_{3}(c^{2}_{12}e^{2i\alpha}m_{2}+m_{1}s^{2}_{12})+e^{-2i\beta}s^{2}_{13}m_{1}m_{2}\right],\quad\\\
(M^{-1}_{\nu})_{12}&=\dfrac{1}{m_{1}m_{2}m_{3}}\left[e^{-i(2(\alpha+\beta)+\delta)}(c_{13}e^{2i\beta}m_{1}m_{3}s_{12}(c_{12}c_{23}e^{i\delta}-s_{12}s_{13}s_{23})\right.\\\
\left.\right.&\left.+c_{13}e^{2i\alpha}m_{2}(-c_{12}c_{23}e^{i(2\beta+\delta)}m_{3}s_{12})+(-c^{2}_{12}e^{2i\beta}m_{3}+m_{1})s_{13}s_{23})\right],\quad\\\
(M^{-1}_{\nu})_{13}&=\dfrac{1}{m_{1}m_{2}m_{3}}\left[e^{-i(2(\alpha+\beta)+\delta)}(c_{13}c_{12}(e^{2i\alpha}m_{1}m_{2}-e^{2i\alpha}m_{3}(c^{2}_{12}e^{2i\alpha}m_{2}+m_{1}s^{2}_{12}))s_{13}\right.\\\
\left.\right.&\left.-c_{12}c_{13}e^{i(2\beta+\delta)}(m_{1}-e^{2i\alpha}m_{2})m_{3}s_{12}s_{23})\right],\quad\\\
(M^{-1}_{\nu})_{22}&=\dfrac{1}{m_{1}m_{2}m_{3}}\left[e^{-2i(\alpha+\delta)}m_{1}m_{3}(c_{12}c_{23}e^{i\delta}-s_{12}s_{13s_{23}})^{2}+m_{2}(c^{2}_{23}m_{3}s^{2}_{12}+2c_{12}c_{23}m_{3}s_{12}s_{13}s_{23}\right.\\\
\left.\right.&\left.+e^{-2i(\beta+\delta)}(c^{2}_{13}m_{1}+c^{2}_{12}e^{i\beta}m_{3}s^{2}_{13})s^{2}_{23})\right],\quad\\\
(M^{-1}_{\nu})_{23}&=\dfrac{1}{m_{1}m_{2}m_{3}}\left[e^{-2i(\alpha+\beta)}m_{1}m_{3}(c_{23}s_{12}s_{13}+c_{12}e^{i\delta}s_{23})(c_{12}c_{23}e^{i\delta}-s_{12}s_{23}s_{13})\right.\\\
\left.\right.&\left.+e^{-2i(\beta+\delta)}m_{2}(c_{23}(c^{2}_{13}m_{1}-e^{2i(\beta+\delta)}m_{3}s^{2}_{12})+c^{2}_{12}c_{23}e^{i\beta}m_{3}s^{2}_{13}s_{23}\right.\\\
\left.\right.&\left.+c_{12}e^{i(2\beta+\delta)}m_{3}s_{12}s_{23}(c^{2}_{23}-s^{2}_{23}))\right],\quad\\\
(M^{-1}_{\nu})_{33}&=\dfrac{1}{m_{1}m_{2}m_{3}}\left[e^{-2i(\alpha+\beta+\delta)}c^{2}_{23}(c^{2}_{13}e^{2i\alpha}m_{1}m_{2}+e^{2i\beta}m_{3}(c^{2}_{12}e^{2i\alpha}m_{2}+m_{1}s^{2}_{12})s^{2}_{13})\right.\\\
\left.\right.&\left.+e^{i(2\alpha+\delta)}2c_{12}c_{23}(m_{1}-e^{2i\alpha}m_{2})m_{3}s_{12}s_{13}s_{23}+e^{-2i\alpha}m_{3}(c^{2}_{12}m_{1}+e^{2i\alpha}m_{2}s^{2}_{12})s^{2}_{23}\right].\end{aligned}\right\\}$
(4)

There are fifteen possible two-zero textures of $M^{-1}_{\nu}$ which we
categorize in five classes viz. class $A,B,C,D$ and $E$ as shown in Table (1).
Symbolically, these 15 textures can be represented as

$A_{1}=\begin{pmatrix}0&0&X\\\ 0&X&X\\\ X&X&X\\\
\end{pmatrix},A_{2}=\begin{pmatrix}0&X&0\\\ X&X&X\\\ 0&X&X\\\
\end{pmatrix},A_{3}=\begin{pmatrix}0&X&X\\\ X&0&X\\\ X&X&X\\\ \end{pmatrix},$
$A_{4}=\begin{pmatrix}0&X&X\\\ X&X&0\\\ X&0&X\\\
\end{pmatrix},A_{5}=\begin{pmatrix}0&X&X\\\ X&X&X\\\ X&X&0\\\ \end{pmatrix};$

$B_{1}=\begin{pmatrix}X&0&0\\\ 0&X&X\\\ 0&X&X\\\
\end{pmatrix},B_{2}=\begin{pmatrix}X&0&X\\\ 0&0&X\\\ X&X&X\\\
\end{pmatrix},B_{3}=\begin{pmatrix}X&0&X\\\ 0&X&0\\\ X&0&X\\\
\end{pmatrix},B_{4}=\begin{pmatrix}X&0&X\\\ 0&X&X\\\ X&X&0\\\ \end{pmatrix};$

$C_{1}=\begin{pmatrix}X&X&0\\\ X&0&X\\\ 0&X&X\\\
\end{pmatrix},C_{2}=\begin{pmatrix}X&X&0\\\ X&X&0\\\ 0&0&X\\\
\end{pmatrix},C_{3}=\begin{pmatrix}X&X&0\\\ X&X&X\\\ 0&X&0\\\ \end{pmatrix};$

$D_{1}=\begin{pmatrix}X&X&X\\\ X&0&0\\\ X&0&X\\\
\end{pmatrix},D_{2}=\begin{pmatrix}X&X&X\\\ X&0&X\\\ X&X&0\\\ \end{pmatrix};$

$E_{1}=\begin{pmatrix}X&X&X\\\ X&X&0\\\ X&0&0\\\ \end{pmatrix},$

where “$X$” denotes the non-zero element of inverse neutrino mass matrix. In
the present work, we investigate the phenomenological consequences of these 15
possible two-zero textures in $M_{\nu}^{-1}$ within the paradigm of LMA and
LMA-D descriptions of neutrino oscillation phenomenon.
In general, two-zero texture of $M_{\nu}^{-1}$ result in constraining
equations

$\displaystyle(M^{-1}_{\nu})_{pq}=0;\hskip 14.22636pt(M^{-1}_{\nu})_{rs}=0$
(5)

where $p,q,r$ and $s$ can take value 1,2,3 such that $p\leq q$, $r\leq s$.
Using Eqn. (3), the constraints described by Eqn. (5) can be written as

$\lambda_{1}{U^{-1}}_{1p}{U^{-1}}_{1q}+\lambda_{2}{U^{-1}}_{2p}{U^{-1}}_{2q}+\lambda_{3}{U^{-1}}_{3p}{U^{-1}}_{3q}=0,$
(6)

and

$\lambda_{1}{U^{-1}}_{1r}{U^{-1}}_{1s}+\lambda_{2}{U^{-1}}_{2r}{U^{-1}}_{2s}+\lambda_{3}{U^{-1}}_{3r}{U^{-1}}_{3s}=0,$
(7)

where

$\displaystyle\lambda_{1}=\dfrac{1}{m_{1}},\hskip
14.22636pt\lambda_{2}=\dfrac{e^{-2i\alpha}}{m_{2}},\hskip
14.22636pt\lambda_{3}=\dfrac{e^{-2i(\beta+\delta)}}{m_{3}}.$

Class | Textures | $(M^{-1}_{\nu})_{pq}=0$ | $(M^{-1}_{\nu})_{rs}=0$
---|---|---|---
A | $A_{1}$ | $(M^{-1}_{\nu})_{11}$ | $(M^{-1}_{\nu})_{12}$
$A_{2}$ | $(M^{-1}_{\nu})_{11}$ | $(M^{-1}_{\nu})_{13}$
$A_{3}$ | $(M^{-1}_{\nu})_{11}$ | $(M^{-1}_{\nu})_{22}$
$A_{4}$ | $(M^{-1}_{\nu})_{11}$ | $(M^{-1}_{\nu})_{23}$
$A_{5}$ | $(M^{-1}_{\nu})_{11}$ | $(M^{-1}_{\nu})_{33}$
B | $B_{1}$ | $(M^{-1}_{\nu})_{12}$ | $(M^{-1}_{\nu})_{13}$
$B_{2}$ | $(M^{-1}_{\nu})_{12}$ | $(M^{-1}_{\nu})_{22}$
$B_{3}$ | $(M^{-1}_{\nu})_{12}$ | $(M^{-1}_{\nu})_{23}$
$B_{4}$ | $(M^{-1}_{\nu})_{12}$ | $(M^{-1}_{\nu})_{33}$
C | $C_{1}$ | $(M^{-1}_{\nu})_{13}$ | $(M^{-1}_{\nu})_{22}$
$C_{2}$ | $(M^{-1}_{\nu})_{13}$ | $(M^{-1}_{\nu})_{23}$
$C_{3}$ | $(M^{-1}_{\nu})_{13}$ | $(M^{-1}_{\nu})_{33}$
D | $D_{1}$ | $(M^{-1}_{\nu})_{22}$ | $(M^{-1}_{\nu})_{23}$
$D_{2}$ | $(M^{-1}_{\nu})_{22}$ | $(M^{-1}_{\nu})_{33}$
E | $E_{1}$ | $(M^{-1}_{\nu})_{23}$ | $(M^{-1}_{\nu})_{33}$
Table 1: Fifteen possible two-zero texture patterns with corresponding first
and second zero.

We solve Eqns. (6) and (7) for two mass ratios
$\left(\dfrac{m_{1}}{m_{2}}e^{-2i\alpha},\hskip
1.42271pt\dfrac{m_{1}}{m_{3}}e^{-2i(\beta+\delta)}\right)$

$\displaystyle\dfrac{m_{1}}{m_{2}}e^{-2i\alpha}=\dfrac{{U^{-1}}_{3r}{U^{-1}}_{3s}{U^{-1}}_{1p}{U^{-1}}_{1q}-{U^{-1}}_{3p}{U^{-1}}_{3q}{U^{-1}}_{1r}{U^{-1}}_{1s}}{{U^{-1}}_{3p}{U^{-1}}_{3q}{U^{-1}}_{2r}{U^{-1}}_{2s}-{U^{-1}}_{2p}{U^{-1}}_{2q}{U^{-1}}_{3r}{U^{-1}}_{3s}},$
(8)
$\displaystyle\dfrac{m_{1}}{m_{3}}e^{-2i(\beta+\delta)}=\dfrac{{U^{-1}}_{1r}{U^{-1}}_{1s}{U^{-1}}_{2p}{U^{-1}}_{2q}-{U^{-1}}_{1p}{U^{-1}}_{1q}{U^{-1}}_{2r}{U^{-1}}_{2s}}{{U^{-1}}_{3p}{U^{-1}}_{3q}{U^{-1}}_{2r}{U^{-1}}_{2s}-{U^{-1}}_{2p}{U^{-1}}_{2q}{U^{-1}}_{3r}{U^{-1}}_{3s}}.$
(9)

The absolute values of mass ratios are given by

$\frac{m_{1}}{m_{2}}=\left|\dfrac{{U^{-1}}_{3r}{U^{-1}}_{3s}{U^{-1}}_{1p}{U^{-1}}_{1q}-{U^{-1}}_{3p}{U^{-1}}_{3q}{U^{-1}}_{1r}{U^{-1}}_{1s}}{{U^{-1}}_{3p}{U^{-1}}_{3q}{U^{-1}}_{2r}{U^{-1}}_{2s}-{U^{-1}}_{2p}{U^{-1}}_{2q}{U^{-1}}_{3r}{U^{-1}}_{3s}}\right|,$
(10)
$\frac{m_{1}}{m_{3}}=\left|\dfrac{{U^{-1}}_{1r}{U^{-1}}_{1s}{U^{-1}}_{2p}{U^{-1}}_{2q}-{U^{-1}}_{1p}{U^{-1}}_{1q}{U^{-1}}_{2r}{U^{-1}}_{2s}}{{U^{-1}}_{3p}{U^{-1}}_{3q}{U^{-1}}_{2r}{U^{-1}}_{2s}-{U^{-1}}_{2p}{U^{-1}}_{2q}{U^{-1}}_{3r}{U^{-1}}_{3s}}\right|,$
(11)

Class | Texture |  Mass ratios $\left(R_{12},R_{13}\right)$
---|---|---
| $A_{1}$ | $\dfrac{m_{1}}{m_{2}}e^{-2i\alpha}=\dfrac{c_{12}\left(c_{23}e^{i\delta}s_{12}s_{13}+c_{12}s_{23}\right)}{s_{12}\left(c_{12}c_{23}e^{i\delta}s_{13}-s_{12}s_{23}\right)}$
|
$\dfrac{m_{1}}{m_{3}}e^{-i(2\beta+\delta)}=\dfrac{c_{12}c^{2}_{13}c_{23}}{s_{13}\left(-c_{12}c_{23}e^{i\delta}s_{13}+s_{12}s_{23}\right)}$
A | $A_{2}$ | $\dfrac{m_{1}}{m_{2}}e^{-2i\alpha}=-\dfrac{c_{12}\left(c_{12}c_{23}-e^{i\delta}s_{12}s_{13}s_{23}\right)}{s_{12}\left(c_{23}s_{12}+c_{12}e^{i\delta}s_{13}s_{23}\right)}$
$\dfrac{m_{1}}{m_{3}}e^{-i(2\beta+\delta)}=\dfrac{c_{12}c^{2}_{13}s_{23}}{s_{13}\left(c_{23}s_{12}+c_{12}e^{i\delta}s_{13}s_{23}\right)}$
$A_{3}$ | $\dfrac{m_{1}}{m_{2}}e^{-2i\alpha}=-\dfrac{\left(c_{23}e^{i\delta}s_{12}s_{13}+c_{12}s_{23}\right)\left(c_{23}e^{i\delta}s_{12}s_{13}+c_{12}(s^{2}_{13}-c^{2}_{13})s_{23}\right)}{c^{2}_{12}c^{2}_{23}e^{2i\delta}s^{2}_{13}-2c_{12}c_{23}e^{i\delta}s_{12}s^{2}_{13}s_{23}+s^{2}_{12}(s^{4}_{13}-c^{4}_{13})s_{23}}$
$\dfrac{m_{1}}{m_{3}}e^{-i(2\beta+\delta)}=\dfrac{c^{2}_{13}c_{23}\left(c_{23}e^{i\delta}(c^{2}_{12}-s^{2}_{12})-2c_{12}s_{12}s_{13}s_{23}\right)}{-c^{2}_{12}c^{2}_{23}e^{2i\delta}s^{2}_{13}+2c_{12}c_{23}e^{i\delta}s_{12}s^{2}_{13}s_{23}+s^{2}_{12}(s^{4}_{13}-c^{4}_{13})s_{23}}$
$A_{4}$ | $\dfrac{m_{1}}{m_{2}}e^{-2i\alpha}=\dfrac{c^{2}_{12}c^{2}_{23}(s^{4}_{13}-c^{4}_{13})-2c_{12}c_{23}e^{i\delta}s_{12}s^{3}_{13}s_{23}+e^{2i\delta}s^{2}_{12}s^{2}_{13}s^{2}_{23}}{c^{2}_{13}c^{2}_{23}s^{2}_{12}-s^{2}_{13}\left(c_{23}s_{12}s_{13}+c_{12}e^{i\delta}s_{23}\right)^{2}}$
|
$\dfrac{m_{1}}{m_{3}}e^{-i(2\beta+\delta)}=-\dfrac{c^{2}_{13}s_{23}\left(2c_{12}c_{23}s_{12}s_{13}+e^{i\delta}(c^{2}_{12}-s^{2}_{12})s_{23}\right)}{-c^{2}_{13}c^{2}_{23}s^{2}_{12}+s^{2}_{13}\left(c_{23}s_{12}s_{13}+c_{12}e^{i\delta}s_{23}\right)^{2}}$
| $A_{5}$ | $\dfrac{m_{1}}{m_{2}}e^{-2i\alpha}=\dfrac{-c_{23}e^{2i\delta}s^{2}_{12}s^{2}_{13}s_{23}+c^{2}_{12}c_{23}(s^{4}_{13}-c^{4}_{13})s_{23}+c_{12}e^{i\delta}s_{12}s^{3}_{13}(c^{2}_{13}-s^{2}_{13})}{c^{2}_{12}c_{23}e^{2i\delta}s^{2}_{13}s_{23}+c_{23}s^{2}_{12}(c^{4}_{13}-s^{4}_{13})s_{23}+c_{12}e^{i\delta}s_{12}s^{3}_{13}(c^{2}_{23}-s^{2}_{23})}$
|
$\dfrac{m_{1}}{m_{3}}e^{-i(2\beta+\delta)}=-\dfrac{-c^{2}_{13}\left(c_{23}e^{i\alpha}(c^{2}_{12}-s^{2}_{12})s_{23}+c_{12}s_{12}s_{13}(c^{2}_{23}-s^{2}_{23})\right)}{c^{2}_{12}c_{23}e^{2i\delta}s^{2}_{13}s_{23}+c_{23}s^{2}_{12}(c^{4}_{13}-s^{4}_{13})s_{23}+c_{12}e^{i\delta}s_{12}s^{3}_{13}(c^{2}_{23}-s^{2}_{23})}$
Table 2: Mass ratios for class $A$. Class | Texture |  Mass ratios $\left(R_{12},R_{13}\right)$
---|---|---
| $B_{1}$ | $\dfrac{m_{1}}{m_{2}}e^{-2i\alpha}=1$
|  | $\dfrac{m_{1}}{m_{3}}e^{-i(2\beta+\delta)}=1$
B | $B_{2}$ | $\dfrac{m_{1}}{m_{2}}e^{-2i\alpha}=-\dfrac{\left(c_{23}e^{i\delta}s_{12}+c_{12}s_{13}s_{23}\right)\left(c_{23}e^{i\delta}s_{12}s_{13}+c_{12}s_{23}\right)}{\left(c_{12}c_{23}e^{i\delta}-s_{12}s_{23}s_{13}\right)\left(c_{12}c_{23}e^{i\delta}s_{13}-s_{12}s_{23}\right)}$
$\dfrac{m_{1}}{m_{3}}e^{-i(2\beta+\delta)}=\dfrac{c_{23}\left(c_{23}e^{i\delta}s_{12}+c_{12}s_{13}s_{23}\right)}{s_{23}\left(c_{12}c_{23}e^{i\delta}s_{13}-s_{12}s_{23}\right)}$
$B_{3}$ | $\dfrac{m_{1}}{m_{2}}e^{-2i\alpha}=\dfrac{\left(c_{23}e^{i\delta}s_{12}+c_{12}s_{13}s_{23}\right)\left(c_{12}c_{23}-e^{i\delta}s_{12}s_{13}s_{23}\right)}{\left(c_{23}s_{12}+c_{12}e^{i\delta}s_{13}s_{23}\right)\left(c_{12}c_{23}e^{i\delta}-s_{12}s_{13}s_{23}\right)}$
$\dfrac{m_{1}}{m_{3}}e^{-i(2\beta+\delta)}=\dfrac{c_{23}e^{i\delta}s_{12}+c_{12}s_{13}s_{23}}{c_{23}s_{12}+c_{12}e^{i\delta}s_{13}s_{23}}$
| $B_{4}$ | $\dfrac{m_{1}}{m_{2}}e^{-2i\alpha}=-\dfrac{\left(c^{2}_{12}c^{2}_{23}s_{13}s_{23}+e^{2i\delta}s^{2}_{12}s_{13}s^{3}_{23}+c_{12}c_{23}e^{i\delta}s_{12}(c^{2}_{13}c^{2}_{23}-s^{2}_{13}s^{2}_{23})\right)}{c^{2}_{23}s^{2}_{12}s_{13}s_{23}+c^{2}_{12}e^{2i\delta}s_{13}s^{3}_{23}-c_{12}c_{23}e^{i\delta}s_{12}(c^{2}_{13}c^{2}_{23}-s^{2}_{13}s^{2}_{23})}$
|
$\dfrac{m_{1}}{m_{3}}e^{-i(2\beta+\delta)}=\dfrac{\left(c_{12}c_{23}e^{2i\delta}s_{12}s^{2}_{23}+e^{i\delta}(c^{2}_{12}-s^{2}_{12})s_{13}s^{3}_{23}+c_{12}c_{23}s_{12}s^{2}_{13}(c^{2}_{23}+2s^{2}_{23})\right)}{c^{2}_{23}s^{2}_{12}s_{13}s_{23}+c^{2}_{12}e^{2i\delta}s_{13}s^{3}_{23}-c_{12}c_{23}e^{i\delta}s_{12}(c^{2}_{13}c^{2}_{23}-s^{2}_{13}s^{2}_{23})}$
Table 3: Mass ratios for class $B$. Class | Texture |  Mass ratios $\left(R_{12},R_{13}\right)$
---|---|---
| $C_{1}$ | $\dfrac{m_{1}}{m_{2}}e^{-2i\alpha}=-\dfrac{\left(c^{3}_{23}e^{2i\delta}s^{2}_{12}s_{13}+c^{2}_{12}c_{23}s_{13}s^{2}_{23}+c_{12}e^{i\delta}s_{12}s_{23}(2c^{2}_{23}s^{2}_{13}-c^{2}_{13}s^{2}_{23})\right)}{c^{2}_{12}c^{3}_{23}s_{13}e^{2i\delta}+c_{23}s^{2}_{12}s_{13}s^{2}_{23}+c_{12}e^{i\delta}s_{12}s_{23}(-2c^{2}_{23}s^{2}_{13}+c^{2}_{13}s^{2}_{23})}$
|
$\dfrac{m_{1}}{m_{3}}e^{-i(2\beta+\delta)}=-\dfrac{\left(c^{3}_{23}e^{i\delta}(c^{2}_{12}-s^{2}_{12})s_{13}+c_{12}c^{2}_{23}e^{2i\delta}s_{12}s_{23}+c_{12}s_{12}s^{2}_{13}s_{23}(2c^{2}_{23}+s^{2}_{23})\right)}{c^{2}_{12}c^{3}_{23}s_{13}e^{2i\delta}+c_{23}s^{2}_{12}s_{13}s^{2}_{23}+c_{12}e^{i\delta}s_{12}s_{23}(-2c^{2}_{23}s^{2}_{13}+c^{2}_{13}s^{2}_{23})}$
C | $C_{2}$ | $\dfrac{m_{1}}{m_{2}}e^{-2i\alpha}=\dfrac{\left(c_{12}c_{23}s_{13}-e^{i\delta}s_{12}s_{23}\right)\left(c_{23}e^{i\delta}s_{12}s_{13}+c_{12}s_{23}\right)}{\left(c_{23}s_{12}s_{13}+c_{12}e^{i\delta}s_{23}\right)\left(c_{12}c_{23}e^{i\delta}s_{13}-s_{12}s_{23}\right)}$
$\dfrac{m_{1}}{m_{3}}e^{-i(2\beta+\delta)}=-\dfrac{\left(-c_{12}c_{23}s_{13}+e^{i\delta}s_{12}s_{23}\right)}{c_{12}c_{23}e^{i\delta}s_{13}-s_{12}s_{23}}$
$C_{3}$ | $\dfrac{m_{1}}{m_{2}}e^{-2i\alpha}=\dfrac{\left(c_{12}c_{23}s_{13}-e^{i\delta}s_{12}s_{23}\right)\left(c_{12}c_{23}-e^{i\delta}s_{12}s_{13}s_{23}\right)}{\left(c_{23}s_{12}s_{13}+c_{12}e^{i\delta}s_{23}\right)\left(c_{23}s_{12}+c_{12}e^{i\delta}s_{13}s_{23}\right)}$
|
$\dfrac{m_{1}}{m_{3}}e^{-i(2\beta+\delta)}=\dfrac{s_{23}\left(c_{12}c_{23}s_{13}-e^{i\delta}s_{12}s_{23}\right)}{c_{23}\left(c_{23}s_{12}+c_{12}e^{i\delta}s_{13}s_{23}\right)}$
| $D_{1}$ | $\dfrac{m_{1}}{m_{2}}e^{-2i\alpha}=-\dfrac{s_{12}\left(c_{23}e^{i\delta}s_{12}+c_{12}s_{13}s_{23}\right)}{c_{12}\left(c_{12}c_{23}e^{i\delta}-s_{12}s_{23}s_{13}\right)}$
D | $\dfrac{m_{1}}{m_{3}}e^{-i(2\beta+\delta)}=\dfrac{s_{13}\left(c_{23}s_{12}+e^{-i\delta}c_{12}s_{13}s_{23}\right)}{c_{12}c^{2}_{13}s_{23}}$
$D_{2}$ | $\dfrac{m_{1}}{m_{2}}e^{-2i\alpha}=-\dfrac{s_{12}\left(2c_{12}c_{23}s_{13}s_{23}+e^{i\delta}s_{12}(c^{2}_{23}-s^{2}_{23})\right)}{c_{12}\left(-2c_{23}s_{12}s_{13}s_{23}+c_{12}e^{i\delta}(c^{2}_{23}-s^{2}_{23})\right)}$
|
$\dfrac{m_{1}}{m_{3}}e^{-i(2\beta+\delta)}=-\dfrac{s_{13}\left(c^{2}_{23}s_{13}(c^{2}_{12}-s^{2}_{12})-2c_{12}c_{23}s_{12}s_{23}(e^{-i\delta}s^{2}_{13}+e^{i\delta})+s_{13}s^{2}_{23}(s^{2}_{12}-c^{2}_{12})\right)}{c_{12}c^{2}_{13}\left(c_{12}c^{2}_{23}e^{i\delta}-2c_{23}s_{12}s_{13}s_{23}-c_{12}e^{i\delta}s^{2}_{23}\right)}$
$E$ | $E_{1}$ | $\dfrac{m_{1}}{m_{2}}e^{-2i\alpha}=\dfrac{s_{12}\left(c_{12}c_{23}s_{13}-e^{i\delta}s_{12}s_{23}\right)}{c_{12}\left(c_{23}s_{12}s_{13}+c_{12}e^{i\delta}s_{23}\right)}$
|  | $\dfrac{m_{1}}{m_{3}}e^{-i(2\beta+\delta)}=\dfrac{s_{13}\left(-e^{-i\delta}c_{12}c_{23}s_{13}+s_{12}s_{23}\right)}{c_{12}c_{23}c^{2}_{13}}$
Table 4: Mass ratios for class $C,D$ and $E$.

and two Majorana phases are obtained as

$\displaystyle\alpha=-\dfrac{1}{2}Arg\left(\dfrac{{U^{-1}}_{3r}{U^{-1}}_{3s}{U^{-1}}_{1p}{U^{-1}}_{1q}-{U^{-1}}_{3p}{U^{-1}}_{3q}{U^{-1}}_{1r}{U^{-1}}_{1s}}{{U^{-1}}_{3p}{U^{-1}}_{3q}{U^{-1}}_{2r}{U^{-1}}_{2s}-{U^{-1}}_{2p}{U^{-1}}_{2q}{U^{-1}}_{3r}{U^{-1}}_{3s}}\right),$
(12)
$\displaystyle\beta=-\dfrac{1}{2}Arg\left(\dfrac{{U^{-1}}_{1r}{U^{-1}}_{1s}{U^{-1}}_{2p}{U^{-1}}_{2q}-{U^{-1}}_{1p}{U^{-1}}_{1q}{U^{-1}}_{2r}{U^{-1}}_{2s}}{{U^{-1}}_{3p}{U^{-1}}_{3q}{U^{-1}}_{2r}{U^{-1}}_{2s}-{U^{-1}}_{2p}{U^{-1}}_{2q}{U^{-1}}_{3r}{U^{-1}}_{3s}}\right)-\delta.$
(13)

It is to be noted that the mass ratios
$\left(\dfrac{m_{1}}{m_{2}}e^{-2i\alpha},\hskip
1.42271pt\dfrac{m_{1}}{m_{3}}e^{-2i(\beta+\delta)}\right)$ are different for
each texture as they depend on the position of zero in $M^{-1}_{\nu}$. For
example, in case of $A_{1}$ texture

$\displaystyle p=1,\hskip 14.22636ptq=1,\hskip 14.22636ptr=1,\hskip
14.22636pts=2,$

therefore, Eqns. (8) and (9) become

$\displaystyle\dfrac{m_{1}}{m_{2}}e^{-2i\alpha}=\dfrac{{U^{-1}}_{31}{U^{-1}}_{32}{U^{-1}}_{11}{U^{-1}}_{11}-{U^{-1}}_{31}{U^{-1}}_{31}{U^{-1}}_{11}{U^{-1}}_{12}}{{U^{-1}}_{31}{U^{-1}}_{31}{U^{-1}}_{21}{U^{-1}}_{22}-{U^{-1}}_{21}{U^{-1}}_{21}{U^{-1}}_{31}{U^{-1}}_{32}},$
$\displaystyle\dfrac{m_{1}}{m_{3}}e^{-2i(\beta+\delta)}=\dfrac{{U^{-1}}_{11}{U^{-1}}_{12}{U^{-1}}_{21}{U^{-1}}_{21}-{U^{-1}}_{11}{U^{-1}}_{11}{U^{-1}}_{21}{U^{-1}}_{22}}{{U^{-1}}_{31}{U^{-1}}_{31}{U^{-1}}_{21}{U^{-1}}_{22}-{U^{-1}}_{21}{U^{-1}}_{21}{U^{-1}}_{31}{U^{-1}}_{32}},$

while for $B_{2}$ texture

$\displaystyle p=1,\hskip 14.22636ptq=2,\hskip 14.22636ptr=2,\hskip
14.22636pts=2,$

resulting in mass ratios

$\displaystyle\dfrac{m_{1}}{m_{2}}e^{-2i\alpha}=\dfrac{{U^{-1}}_{32}{U^{-1}}_{32}{U^{-1}}_{11}{U^{-1}}_{12}-{U^{-1}}_{31}{U^{-1}}_{32}{U^{-1}}_{12}{U^{-1}}_{12}}{{U^{-1}}_{31}{U^{-1}}_{32}{U^{-1}}_{22}{U^{-1}}_{22}-{U^{-1}}_{21}{U^{-1}}_{22}{U^{-1}}_{32}{U^{-1}}_{32}},$
$\displaystyle\dfrac{m_{1}}{m_{3}}e^{-2i(\beta+\delta)}=\dfrac{{U^{-1}}_{12}{U^{-1}}_{12}{U^{-1}}_{21}{U^{-1}}_{22}-{U^{-1}}_{11}{U^{-1}}_{12}{U^{-1}}_{22}{U^{-1}}_{22}}{{U^{-1}}_{31}{U^{-1}}_{32}{U^{-1}}_{22}{U^{-1}}_{22}-{U^{-1}}_{21}{U^{-1}}_{22}{U^{-1}}_{32}{U^{-1}}_{32}}.$

We have analytically solved the mass ratios for all possible two-zero textures
which are shown in the Tables (2-4).
The two mass-squared differences $\Delta m_{21}^{2}=m_{2}^{2}-m_{1}^{2}$ and
$\left|\Delta m_{32}^{2}\right|=m_{3}^{2}-m_{2}^{2}$ alongwith mass ratios
$\dfrac{m_{1}}{m_{2}}e^{-2i\alpha}\equiv R_{12}$,
$\dfrac{m_{1}}{m_{3}}e^{-i(2\beta+\delta)}\equiv R_{13}$ yield two values of
neutrino mass $m_{1}$ given by

$m_{1}^{a}=\left|R_{12}\right|\sqrt{\dfrac{\Delta
m^{2}_{21}}{1-\left|R_{12}\right|^{2}}},\hskip
8.53581ptm_{1}^{b}=\left|R_{13}\right|\sqrt{\dfrac{\Delta m^{2}_{21}+|\Delta
m^{2}_{32}|}{1-\left|R_{13}\right|^{2}}},$ (14)

respectively. In order to ensure the consistency of the formalism the two
values ($m_{1}^{a},m_{1}^{b}$) must be equal which results in mass ratio
parameter

$\displaystyle\dfrac{\Delta m^{2}_{21}}{|\Delta
m^{2}_{32}|}=\dfrac{|R_{13}|^{2}\left(1-|R_{12}|^{2}\right)}{|R_{12}|^{2}-|R_{13}|^{2}}\equiv
R_{\nu}.$ (15)

The $3\sigma$ experimental range of parameter $R_{\nu}$ defined in Eqn. (15)
is 0.02590 $<R_{\nu}<$ 0.03656. The allowed phenomenology of the model is
obtained by restricting $R_{\nu}$ in the $3\sigma$ experimental range. The
neutrino mass eigenvalues $m_{2}$ and $m_{3}$ can be obtained using mass
square differences $\left(\Delta m_{21}^{2},\hskip 1.42271pt\Delta
m_{32}^{2}\right)$ as

$\displaystyle m_{2}=\sqrt{m_{1}^{2}+\Delta m_{21}^{2}};\hskip
2.84544ptm_{3}=\sqrt{m_{2}^{2}+\Delta m_{32}^{2}}\hskip 5.69046pt\text{for
normal hierarchy (NH) $(m_{1}<m_{2}<m_{3})$},$ (16) and (17) $\displaystyle
m_{2}=\sqrt{m_{1}^{2}+\Delta m_{21}^{2}};\hskip
2.84544ptm_{1}=\sqrt{m_{3}^{2}+\Delta m_{32}^{2}-\Delta m_{21}^{2}}\hskip
5.69046pt\text{for inverted hierarchy (IH) $(m_{3}<m_{1}<m_{2})$}.$

Also, the effective Majorana mass which governs the neutrinoless double
beta($0\nu\beta\beta$) decay process is given by

$|M_{ee}|=\left|\sum_{i}V^{2}_{ei}m_{i}\right|=\left|m_{1}c^{2}_{12}c^{2}_{13}+m_{2}s^{2}_{12}c^{2}_{13}e^{2i\alpha}+m_{3}s^{2}_{13}e^{2i\beta}\right|.$
(18)

The Jarlskog $CP$ invariant is defined as [31, 32]

$J_{CP}=s_{23}c_{23}s_{12}c_{12}s_{13}c_{13}^{2}\sin\delta.$

Parameters | best fit$\pm 1\sigma$ range(NH) | best fit$\pm 1\sigma$ range(IH)
---|---|---
$\Delta m^{2}_{21}[10^{-5}$ eV${}^{2}]$ | $7.50^{+0.22}_{-0.20}$ | $7.50^{+0.22}_{-0.20}$
$\left|\Delta m^{2}_{31}\right|[10^{-3}$ eV${}^{2}]$ | $2.55^{+0.02}_{-0.03}$ | $2.45^{+0.02}_{-0.03}$
$\theta^{\circ}_{12}$ | $34.3\pm 1.0$ | $34.3\pm 1.0$
$\theta_{23}^{\circ}$ | $49.26\pm 0.79$ | $49.46^{+0.60}_{-0.97}$
$\theta_{13}^{\circ}$ | $8.53^{+0.13}_{-0.14}$ | $8.53^{+0.12}_{-0.14}$
Table 5: Global fit values of neutrino oscillation parameters[37].

In the numerical analysis, to study LMA phenomenology of the model, we have
randomly generated the known neutrino oscillation parameters such as
$\theta_{ij}$ $(i,j=1,2,3$; $i<j)$ and $\Delta m_{ij}^{2}$ $(i>j)$ using
Gaussian distribution within allowed experimental range shown in Table (5).
However, in order to study the viability of the model under LMA-D solution,
$\theta_{12}$ is randomly generated using the uniform distribution within the
range (53.71∘-58.37∘)[33, 34, 35, 36].

##  3 LMA and LMA-D phenomenology

In this section, we have investigated the phenomenology of all possible two-
zero textures of $M_{\nu}^{-1}$ under the paradigm of LMA and LMA-D solutions
to the neutrino oscillation phenomenon.

###  3.1 Class A

Class $A$ is disallowed for both LMA and LMA-D solutions. As a representative
case, we have discussed the viability of $A_{1}$ texture in the following. The
mass ratios ($|R_{12}|,|R_{13}|$) for $A_{1}$ texture up-to first order in
$s_{13}$ can be written as

$\displaystyle|R_{12}|\equiv\dfrac{m_{1}}{m_{2}}\approx\dfrac{c_{12}^{2}}{s_{12}^{2}}+\dfrac{c_{12}c_{23}\cos{\delta}}{s_{12}^{3}s_{23}}s_{13},$
(19)
$\displaystyle|R_{13}|\equiv\dfrac{m_{1}}{m_{3}}\approx\dfrac{c_{12}^{2}c_{23}^{2}\cos{\delta}}{s_{12}^{2}s_{23}^{2}}+\dfrac{c_{12}c_{23}\left(c_{12}^{2}c_{23}^{2}-4s_{12}^{2}s_{23}^{2}+3c_{12}^{2}c_{23}^{2}\cos{2\delta}\right)}{4s_{12}^{3}s_{23}^{3}}.$
(20)

LMA Scenario: For $\delta$ in the range $0^{\circ}\leq\delta\leq 90^{\circ}$
or $270^{\circ}\leq\delta\leq 360^{\circ}$, Eqn. (19) results in
$\dfrac{m_{1}}{m_{2}}>1$. Furthermore, if $\delta$ lie in the range
$90^{\circ}<\delta<270^{\circ}$ and
$\sin{\theta_{23}}\approx\cos{\theta_{23}}$, solar mass hierarchy requires
$\sin{2{\theta_{12}}}<4\sin{\theta_{13}}$, however, from neutrino oscillation
data (Table (5)) $\sin{2{\theta_{12}}}>4\sin{\theta_{13}}$ implying that
$A_{1}$ texture is disallowed. LMA-D Scenario: The parameter $(R_{\nu})$ up-to
first order in $s_{13}$ can be written as

$\displaystyle
R_{\nu}\approx\left(-1+\dfrac{c_{12}^{4}}{s_{12}^{4}}\right)+\dfrac{2c_{12}^{3}c_{23}\cos{\delta}}{s_{12}^{5}s_{23}}s_{13}.$
(21)

Using $\theta_{12}=55^{\circ}$ and $\delta$= $80^{\circ}$ $(190^{\circ})$,
numerical value of $R_{\nu}$ is found to be $0.74$ $(0.89)$ which lies outside
the $3\sigma$ range of $R_{\nu}$. The above observations are, also, evident
from Fig. (1).

Similar analysis can be done for all remaining textures in class $A$.
Therefore, neutrino mass model, with two-zero textures in $M_{\nu}^{-1}$,
wherein one of the texture zero is at $(1,1)$ place in $M_{\nu}^{-1}$ is
disallowed.

Figure 1: Correlation between $\left|R_{12}\right|$ and $R_{\nu}$ for texture
$A_{1}$.

###  3.2 Class B

In class $B$, mass ratios for textures $B_{1}$ are given by

$\left|R_{12}\right|\equiv\dfrac{m_{1}}{m_{2}}=1,\hskip
8.53581pt\left|R_{13}\right|\equiv\dfrac{m_{1}}{m_{3}}=1,$

which result in degenerate neutrino masses (i.e. $m_{1}=m_{2}=m_{3}$) and are
inconsistent with neutrino oscillation data. Also, for texture $B_{3}$, mass
ratios are given by

$\displaystyle|R_{12}|=\left|\dfrac{\left(c_{23}e^{i\delta}s_{12}+c_{12}s_{13}s_{23}\right)\left(c_{12}c_{23}-e^{i\delta}s_{12}s_{13}s_{23}\right)}{\left(c_{23}s_{12}+c_{12}e^{i\delta}s_{13}s_{23}\right)\left(c_{12}c_{23}e^{i\delta}-s_{12}s_{13}s_{23}\right)}\right|=1,$
(22)
$\displaystyle|R_{13}|=\left|\dfrac{c_{23}e^{i\delta}s_{12}+c_{12}s_{13}s_{23}}{c_{23}s_{12}+c_{12}e^{i\delta}s_{13}s_{23}}\right|=1,$
(23)

resulting in degenerate neutrino masses. Therefore, textures $B_{1}$ and
$B_{3}$ are disallowed. In the following we have investigated the
phenomenological consequences of $B_{2}$ and $B_{4}$ textures under LMA and
LMA-D solutions.
The mass ratios ($|R_{12}|,|R_{13}|$) for texture $B_{2}$, up-to first order
in $s_{13}$, can be written as

$\displaystyle\left|R_{12}\right|\equiv\dfrac{m_{1}}{m_{2}}\approx
1+\dfrac{2\cos{\delta}}{c_{12}c_{23}s_{12}s_{23}}s_{13},$ (24)
$\displaystyle\left|R_{13}\right|\equiv\dfrac{m_{1}}{m_{3}}\approx\dfrac{c_{23}^{2}}{s_{23}^{2}}+\dfrac{c_{12}c_{23}\cos{\delta}}{s_{12}s_{23}^{3}}s_{13}.$
(25)

From Eqn. (24), in order to satisfy the solar mass hierarchy $i.e$
$|R_{12}|\equiv\frac{m_{1}}{m_{2}}<1$, $\delta$ should be in the range
$90^{\circ}<\delta<270^{\circ}$. Also, texture $B_{2}$ predicts the normal
hierarchical neutrino masses for $\theta_{23}$ above maximality
($\theta_{23}>45^{o}$) and $\cos\delta$ negative such that
$\left|R_{13}\right|\equiv\dfrac{m_{1}}{m_{3}}<1$. The above observations are,
also, depicted in the Fig. (2). Similar analysis can also be done for texture
$B_{4}$.
The texture $B_{2}$ ($B_{4}$) is allowed for both LMA and LMA-D solutions with
normal (inverted) hierarchy. The allowed parameter space for these textures
are shown in Fig. (2) as correlation plots amongst different parameters. Both
these textures are found to have identical phenomenology under LMA and LMA-D
solutions. In Fig. (2) we have shown the LMA scenario for $B_{2}$ (NH) and
$B_{4}$ (IH) textures. In the left (right) panel we have depicted the
correlation plots for $B_{2}$ ($B_{4}$) texture. The atmospheric mixing angle
$\theta_{23}$ is found to be above maximality for both the textures. The $CP$
violating phases $\alpha,\beta$ and $\delta$ are found to be sharply
constrained. The Dirac-type $CP$ violating phase $\delta$ is found to be
maximal (around $90^{o}$ and $270^{o}$) and the Jarlskog rephasing invariant
$J_{CP}\neq 0$, thus, these textures are necessarily $CP$ violating. We can,
also, appreciate the $CP$ violating nature of these textures analytically. For
example, for texture $B_{2}$ with $\delta=0^{o}$ ($CP$ conserving scenario),
we obtain $R_{\nu}$ to the first order in $s_{13}$ using Eqn. (15) and values
of mass ratios given in Table (3)

$\displaystyle
R_{\nu}\approx\frac{2c_{23}}{c_{12}s_{12}s_{23}\left(c_{23}^{2}-s_{23}^{2}\right)}s_{13}-\frac{2c_{23}s_{23}}{c_{12}s_{12}\left(c_{23}^{2}-s_{23}^{2}\right)}s_{13}.$

Using the best-fit values given in Table (5), $|R_{\nu}|\approx 1.7$ which is
outside 3$\sigma$ range, thus, $\delta=0^{o}$ is disallowed implying $B_{2}$
texture is necessarily $CP$ violating. Similarly, for texture $B_{4}$, taking
$\delta=0^{\circ}$ $R_{\nu}$ can, approximately, be written as

$\displaystyle
R_{\nu}\approx\frac{s_{23}^{5}\left(2c_{23}^{4}c_{12}^{2}-2c_{12}^{2}s_{23}^{2}c_{23}^{2}+2c_{23}^{2}s_{12}^{2}-s_{23}^{2}\right)}{c_{12}c_{23}^{5}s_{12}\left(s_{23}^{4}-c_{23}^{4}\right)}s_{13},$

which again result in $|R_{\nu}|>1$ for best-fit values of the mixing angles.
Similar analysis can be done for $CP$-violating textures of class C and D
discussed in the following sections.

We have, also, obtained the implication of the model for neutrinoless double
beta ($0\nu\beta\beta$) decay amplitude $|M_{ee}|$. It is evident from Fig.
(2) that there exist a lower bound on $|M_{ee}|$ in both the textures. For
texture $B_{2}$ ($B_{4}$), $|M_{ee}|<0.03$ eV ($0.06$ eV). The prediction for
$|M_{ee}|$ has, also, been tabulated in Table (7).

###  3.3 Class C

For texture $C_{2}$, the mass ratios $|R_{12}|$ and $|R_{13}|$ are equal to 1
resulting in degenerate neutrino masses which is in contradiction with
neutrino oscillation data. Therefore, texture $C_{2}$ is disallowed. For
textures $C_{1}$, the mass ratios $|R_{12}|$ and $|R_{13}|$, up-to first order
in $s_{13}$, are given by

$\displaystyle|R_{12}|\equiv\dfrac{m_{1}}{m_{2}}\approx
1-\frac{c_{23}s_{13}\cos{\delta}}{c_{12}s_{12}s_{23}^{3}},$ (26)
$\displaystyle|R_{13}|\equiv\dfrac{m_{1}}{m_{3}}\approx\frac{c_{23}^{2}}{s_{23}^{2}}-\frac{c_{12}c_{23}^{3}s_{13}\cos\delta}{s_{12}s_{23}^{5}}.$
(27)

It can be seen from Eqn. (26) that $\delta$ should lie in first and fourth
quadrant to have $|R_{12}|\equiv\frac{m_{1}}{m_{2}}<1$. Also, for
$\theta_{23}$ above maximality ($\theta_{23}>45^{o}$) and $\cos\delta$
positive the model predicts normal hierarchical neutrino masses. The above
observations are, also, supplemented by Fig. (3). Similar analysis can, also,
be done for texture $C_{3}$. In fact, it is evident from Fig. (3) that $C_{3}$
admit inverted hierarchical neutrino masses. Also, it is to be noted that the
phenomenology of both the textures are found to be similar under LMA and LMA-D
solutions. In Fig. (3), we have shown the allowed parameter space considering
LMA solution. The Majorana phases are sharply correlated and constrained to
very narrow ranges giving a lower bound on $0\nu\beta\beta$ decay amplitude
$|M_{ee}|$. For texture $C_{1}$ ($C_{3}$) $|M_{ee}|>0.02$ $(0.06)$ eV (Table
(7)). The Dirac-type $CP$ violating phase $\delta$ is sharply constrained
around $90^{o}$ and $270^{o}$, thus, allowing a maximal $CP$ violation. The
Jarlskog rephasing invariant $J_{CP}$ is non-zero which is, also, shown in
Fig. (3).

###  3.4 Class D

For textures $D_{1}$, the mass ratios, up-to first order in $s_{13}$, are
given by

$\displaystyle|R_{12}|\equiv\dfrac{m_{1}}{m_{2}}\approx\frac{s_{12}^{2}}{c_{12}^{2}}+\frac{s_{12}s_{23}s_{13}\cos\delta}{c_{12}^{3}c_{23}},$
(28)
$\displaystyle|R_{13}|\equiv\dfrac{m_{1}}{m_{3}}\approx\frac{\tan\theta_{12}}{\tan\theta_{23}}s_{13}.$
(29)

For texture $D_{1}$, LMA-D solution is disallowed as $|R_{12}|>1$ (Eqn. (28)).
However, for LMA solution, it is evident that $|R_{13}|<1$(Eqn. (29)) implying
normal hierarchical neutrino masses. The above analytical observations are,
also, supplemented by the correlation plots in Figs. (4) and (5). In addition,
the Majorana phases are sharply constrained and correlated in such a way
giving vanishing value of $0\nu\beta\beta$ decay amplitude $|M_{ee}|$(Table
(7)).

Texture $D_{2}$ is found to be consistent with both LMA and LMA-D descriptions
with normal hierarchical neutrino masses, as shown in Figs. (6) and (7). It is
evident from ($\theta_{23}-|M_{ee}|$) correlation plot in Fig. (7) that there
exist a $3\sigma$ lower bound on $0\nu\beta\beta$ decay amplitude
$|M_{ee}|>0.02$ eV (see Table (7)). Furthermore, in contrast to $D_{1}$,
texture $D_{2}$ is found to be necessarily $CP$ violating as depicted in
($\theta_{13}$-$J_{CP}$) (Fig. (7)) and ($\theta_{23}$-$\delta$) (Fig.
(8))correlation plots.

###  3.5 Class E

The mass ratios for texture $E_{1}$ can be obtained from Eqns. (28) and (29)
by using the transformation $c_{23}\rightarrow s_{23};s_{23}\rightarrow
c_{23}$ viz.,

$\displaystyle|R_{12}|\equiv\dfrac{m_{1}}{m_{2}}\approx\frac{s_{12}^{2}}{c_{12}^{2}}+\frac{s_{12}c_{23}s_{13}\cos\delta}{c_{12}^{3}s_{23}},$
(30)
$\displaystyle|R_{13}|\equiv\dfrac{m_{1}}{m_{3}}\approx\tan\theta_{12}\tan\theta_{23}s_{13}.$
(31)

It can be seen from Eqn. (30) that $E_{1}$ is not consistent with LMA-D
solution as $|R_{12}|>1$. The allowed parameter space for LMA with NH is shown
in Figs. (9) and (10). Also, the Majorana phases are correlated in such a way
that $0\nu\beta\beta$ decay amplitude $|M_{ee}|$ is found to be vanishing in
this case. Furthermore, the texture allows for both $CP$ conserving and
violating solutions as is evident from Fig. (10).

Class | Texture | LMA | LMA-D
---|---|---|---
A |  $A_{1}$ (NH/IH) |  $\times$/$\times$ |  $\times$/$\times$
$A_{2}$ (NH/IH) |  $\times$/$\times$ |  $\times$/$\times$
$A_{3}$ (NH/IH) |  $\times$/$\times$ |  $\times$/$\times$
$A_{4}$ (NH/IH) |  $\times$/$\times$ |  $\times$/$\times$
$A_{5}$ (NH/IH) |  $\times$/$\times$ |  $\times$/$\times$
B |  $B_{1}$ (NH/IH) |  $\times$/$\times$ |  $\times$/$\times$
$B_{2}$ (NH/IH) |  $\checkmark$/$\times$ |  $\checkmark$/$\times$
$B_{3}$ (NH/IH) |  $\times$/$\times$ |  $\times$/$\times$
$B_{4}$ (NH/IH) |  $\times$/$\checkmark$ |  $\times$/$\checkmark$
C |  $C_{1}$ (NH/IH) |  $\checkmark$/$\times$ |  $\checkmark$/$\times$
$C_{2}$ (NH/IH) |  $\times$/$\times$ |  $\times$/$\times$
$C_{3}$ (NH/IH) |  $\times$/$\checkmark$ |  $\times$/$\checkmark$
D |  $D_{1}$ (NH/IH) |  $\checkmark$/$\times$ |  $\times$/$\times$
$D_{2}$ (NH/IH) |  $\checkmark$/$\times$ |  $\checkmark$/$\times$
E |  $E_{1}$ (NH/IH) |  $\checkmark$/$\times$ |  $\times$/$\times$
Table 6: Allowed/disallowed two-zero textures of $M_{\nu}^{-1}$ under LMA and LMA-D solutions. The $\checkmark$ ($\times$) mark is used to denote allowed (disallowed) texture. Allowed Textures | Lower bound on $|M_{ee}|$ (eV)
---|---
LMA | LMA-D
NH | IH | NH | IH
$B_{2}$ | 0.02 | $\times$ | 0.04 | $\times$
$B_{4}$ | $\times$ | 0.06 | $\times$ | 0.06
$C_{1}$ | 0.02 | $\times$ | 0.02 | $\times$
$C_{3}$ | $\times$ | 0.06 | $\times$ | 0.06
$D_{1}$ | 0 | $\times$ | $\times$ | $\times$
$D_{2}$ | 0.01 | $\times$ | 0.02 | $\times$
$E_{1}$ | 0 | $\times$ | $\times$ | $\times$
Table 7: $3\sigma$ lower bound on effective Majorana mass ($|M_{ee}|$) in eV
for all allowed textures. “$\times$” symbolize the disallowed hierarchy in the
corresponding texture.

Figure 2: Correlation plots for textures $B_{2}$-NH(left panel) and
$B_{4}$-IH(right panel) under LMA scenario.

Figure 3: Correlation plots for textures $C_{1}$-NH(left panel) and
$C_{3}$-IH(right panel) under LMA scenario. Figure 4: Correlation plots
between $(\theta_{23}-\delta)$(left) and
$(\left|R_{12}\right|-R_{\nu})$(right) for $D_{1}$ texture under LMA solution
with NH. Figure 5: Correlation plots for $D_{1}$ texture under LMA solution
with NH.

Figure 6: Correlation between ($\left|R_{12}\right|-R_{\nu}$) for $D_{2}$
texture. Figure 7: Correlation plots for $D_{2}$ texture under LMA solution
with NH. Figure 8: Correlation between $(\theta_{23}-\delta)$ for $D_{2}$
texture with NH under LMA(left) and LMA-D(right). Figure 9: Correlation
between ($\left|R_{12}\right|-R_{\nu}$) for texture $E_{1}$ under LMA solution
with NH. Figure 10: Correlation plots for texture $E_{1}$ under LMA solution
with NH.

##  4 Symmetry Realization

In this section, we discuss the minimal realization of two-zero texture of
$M_{\nu}^{-1}$. Motivated by the pivotal character played by the discrete
flavor symmetry groups in explaining the observed neutrino oscillation
data[39, 38], we obtain the $A_{4}$ flavor group based symmetry realisation of
inverse neutrino mass matrix ($M_{\nu}^{-1}$) by extending the standard model
particle content in the lepton sector. $A_{4}$ is a non-Abelian discrete group
of even permutations. $A_{4}$ group is a orientation-preserving symmetry of a
regular tetrahedron. It has four irreducible representations 1, $1^{\prime}$,
$1^{\prime\prime}$ and 3 and can be generated using S and T generators
satisfying the relations

$S^{2}=T^{3}=(ST)^{3}=1.$

Here, we choose basis for $A_{4}$ group in which T generator takes the
diagonal form. The reason behind choosing this particular representation is
that it facilitates the diagonal mass matrix for charged leptons. In
T-diagonal basis, one dimensional unitary representation 1, $1^{\prime}$ and
$1^{\prime\prime}$ with generator S and T can be written as

$\displaystyle 1:\hskip 14.22636ptS=1,\hskip 2.84544ptT=1,$ $\displaystyle
1^{\prime}:\hskip 14.22636ptS=1,\hskip 2.84544ptT=\omega,$ $\displaystyle
1^{\prime\prime}:\hskip 14.22636ptS=1,\hskip 2.84544ptT=\omega^{2},$

such that $\omega$= $e^{i2\pi/3}$ whereas three-dimensional unitary
representation is given by

$T=\begin{pmatrix}1&0&0\\\ 0&\omega&0\\\ 0&0&\omega^{2}\\\
\end{pmatrix},\hskip 8.53581ptS=\dfrac{1}{3}\begin{pmatrix}-1&2&2\\\ 2&-1&2\\\
2&2&-1\\\ \end{pmatrix}.$

The multiplication rules for the representations of $A_{4}$ are as follows

$\displaystyle 1^{\prime}\otimes 1^{\prime}=1^{\prime\prime},\hskip
14.22636pt1^{\prime\prime}\otimes 1^{\prime\prime}=1^{\prime},\hskip
14.22636pt1^{\prime}\otimes 1^{\prime\prime}=1,\hskip
14.22636pt1^{\prime\prime}\otimes 1=1^{\prime\prime},$ $\displaystyle 1\otimes
1^{\prime}=1^{\prime},\hskip 14.22636pt3\otimes 1^{\prime}=3,\hskip
14.22636pt3\otimes 1^{\prime\prime}=3,\hskip 14.22636pt3\otimes 3=1\oplus
1^{\prime}\oplus 1^{\prime\prime}\oplus 3_{s}\oplus 3_{a}.$

In the $T$-diagonal basis, the Clebsch–Gordan decomposition of two triplets, a
= ($a_{1},a_{2},a_{3}$) and b = ($b_{1},b_{2},b_{3}$) is given as

$\displaystyle(a\otimes b)_{1}=a_{1}b_{1}+a_{2}b_{3}+a_{3}b_{2},$
$\displaystyle(a\otimes b)_{1^{\prime}}=a_{3}b_{3}+a_{1}b_{2}+a_{2}b_{1},$
$\displaystyle(a\otimes
b)_{1^{\prime\prime}}=a_{2}b_{2}+a_{1}b_{3}+a_{3}b_{1},$
$\displaystyle(a\otimes
b)_{3_{s}}=\dfrac{1}{3}\left(2a_{1}b_{1}-a_{2}b_{3}-a_{3}b_{2},2a_{3}b_{3}-a_{1}b_{2}-a_{2}b_{1},2a_{2}b_{2}-a_{1}b_{3}-a_{3}b_{1}\right),$
$\displaystyle(a\otimes
b)_{3_{a}}=\dfrac{1}{2}\left(a_{2}b_{3}-a_{3}b_{2},a_{1}b_{2}-a_{2}b_{1},a_{1}b_{3}-a_{3}b_{1}\right).$
(32)

Symmetry | ${D}_{eL}$ | ${D}_{\mu L}$ | ${D}_{\tau L}$ | $e_{R}$ | $\mu_{R}$ | $\tau_{R}$ | $\nu_{1R}$ | $\nu_{2R}$ | $\nu_{3R}$ | $\phi$ | $\chi$ | Obtained two-zero textures
---|---|---|---|---|---|---|---|---|---|---|---|---
$SU(2)_{L}$ | 2 | 2 | 2 | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 1 |
$A_{4}$ | 1 | $1^{\prime\prime}$ | $1^{\prime}$ | 1 | $1^{\prime}$ | $1^{\prime\prime}$ | 1 | $1^{\prime}$ | $1^{\prime\prime}$ | 1 | $1^{\prime}$ |  $B_{4}$= $\begin{pmatrix}X&0&X\\\ 0&X&X\\\ X&X&0\\\ \end{pmatrix}$
1 | $1^{\prime\prime}$ | $1^{\prime}$ | 1 | $1^{\prime}$ | $1^{\prime\prime}$ | 1 | $1^{\prime}$ | $1^{\prime\prime}$ | 1 | $1^{\prime\prime}$ |  $C_{1}$=$\begin{pmatrix}X&X&0\\\ X&0&X\\\ 0&X&X\\\ \end{pmatrix}$
Table 8: Field content and charge assignments in the model under $SU(2)_{L}$
and $A_{4}$ symmetries.

Here, we have worked in the framework of Type-I seesaw. In addition, we have
minimally extended the standard model by adding three right-handed neutrino
fields ($\nu_{iR}$; $i=1,2,3$) and one scalar field ($\chi$), having singlet
representation under $A_{4}$ symmetry as shown in the Table (7). In general,
for any Yukawa coupling to be non-zero, its Yukawa Lagrangian term must be in
singlet-invariant representation of $A_{4}$ with mass dimension four at tree
level.
Texture $B_{4}$:
Using the tensor products in Eqn. (4), the invariant Yukawa Lagrangian is
given by

$\displaystyle-\mathcal{L}=........$ $\displaystyle+y_{e}\bar{D}_{eL}\phi
e_{R}+y_{\mu}\bar{D}_{\mu L}\phi\mu_{R}+y_{\tau}\bar{D}_{\tau L}\phi\tau_{R}$
(33) $\displaystyle+y_{1}\bar{D}_{eL}\tilde{\phi}\nu_{eR}+y_{2}\bar{D}_{\mu
L}\tilde{\phi}\nu_{\mu R}+y_{3}\bar{D}_{\tau L}\tilde{\phi}\nu_{\mu R}$
$\displaystyle+\dfrac{1}{2}\left[M_{1}(\nu^{T}_{1R}C^{-1}\nu_{1R})+M_{2}(\nu^{T}_{2R}C^{-1}\nu_{3R}+\nu^{T}_{3R}C^{-1}\nu_{2R})\right]$
$\displaystyle+\dfrac{1}{2}\left[\left(y_{\chi_{1}}(\nu^{T}_{2R}C^{-1}\nu_{2R})+y_{\chi_{2}}(\nu^{T}_{1R}C^{-1}\nu_{3R}+\nu^{T}_{3R}C^{-1}\nu_{1R})\right)\chi\right],$

where $y_{k},y_{i}$ ($k=e,\mu,\tau$; $i=1,2,3$) are Yukawa coupling constants,
$M_{1,2}$ are bare mass terms for right-handed Majorana neutrinos,
$y_{\chi_{1,2}}$ denotes Yukawa coupling constant for interaction terms with
scalar field $\chi$ and $\tilde{\phi}=i\tau_{2}\phi^{*}$; $\tau_{2}$ being
Pauli matrix. The dots in the Lagrangian represents the other kinetic and
scalar potential terms. We have restricted up to Yukawa interactions
pertaining to mass terms. Spontaneous symmetry breaking (SSB) occurred with
vacuum expectation values ($vev$’s) $v$ and $w$ for the Higgs doublet and
scalar singlet field, respectively. The Yukawa Lagrangian (Eqn. (33)) leads to
the mass matrices as

$\displaystyle M_{l}=\begin{pmatrix}y_{e}v&0&0\\\ 0&y_{\mu}v&0\\\
0&0&y_{\tau}v\\\ \end{pmatrix},\hskip
28.45274ptM_{D}=\begin{pmatrix}y_{1}v&0&0\\\ 0&y_{2}v&0\\\ 0&0&y_{3}v\\\
\end{pmatrix}$ (34)

and

$\displaystyle M_{R}=\begin{pmatrix}M_{1}&0&y_{\chi_{2}}w\\\
0&y_{\chi_{1}}w&M_{2}\\\ y_{\chi_{2}}w&M_{2}&0\\\ \end{pmatrix},$ (35)

where $M_{l}$, $M_{D}$ and $M_{R}$ corresponds to charged lepton mass matrix,
Dirac mass matrix and right-handed Majorana mass matrix, respectively.
Type-I seesaw contribution to effective Majorana neutrino mass matrix is given
by

$\displaystyle M_{\nu}=M_{D}M_{R}^{-1}M_{D}^{T}.$ (36)

Also, the inverse neutrino mass matrix can be written as

$\displaystyle M_{\nu}^{-1}=M_{D}^{-T}M_{R}M_{D}^{-1}.$ (37)

In $M_{D}$-diagonal basis, the peculiar feature of implementation of type-I
seesaw for $M_{\nu}^{-1}$ is that the zero(s) in $M_{R}$ corresponds to
zero(s) in $M_{\nu}^{-1}$. Using the Eqns. (34) and (35), the $M_{\nu}^{-1}$
is given by

$\displaystyle
M^{-1}_{\nu}=\begin{pmatrix}\dfrac{M_{1}}{v^{2}y^{2}_{1}}&0&\dfrac{y_{\chi_{1}}w}{v^{2}y_{1}y_{3}}\\\
0&\dfrac{y_{\chi_{2}}w}{v^{2}y^{2}_{3}}&\dfrac{M_{2}}{v^{2}y_{2}y_{3}}\\\
\dfrac{y_{\chi_{1}}w}{v^{2}y_{1}y_{3}}&\dfrac{M_{2}}{v^{2}y_{2}y_{3}}&0\\\
\end{pmatrix},$ (38)

which symbolically can be written as

$\displaystyle M_{\nu}^{-1}=\begin{pmatrix}X&0&X\\\ 0&X&X\\\
X&X&0\end{pmatrix}$ (39)

corresponding to texture $B_{4}$.
Texture $C_{1}$:
For the realization of texture $C_{1}$, we change the irreducible
representation of scalar field $\chi$ to be $1^{\prime\prime}$. The relevant
Yukawa Lagrangian is

$\displaystyle-\mathcal{L}=........$ $\displaystyle+y_{e}\bar{D}_{eL}\phi
e_{R}+y_{\mu}\bar{D}_{\mu L}\phi\mu_{R}+y_{\tau}\bar{D}_{\tau L}\phi\tau_{R}$
(40) $\displaystyle+y_{1}\bar{D}_{eL}\tilde{\phi}\nu_{eR}+y_{2}\bar{D}_{\mu
L}\tilde{\phi}\nu_{\mu R}+y_{3}\bar{D}_{\tau L}\tilde{\phi}\nu_{\mu R}$
$\displaystyle+\dfrac{1}{2}\left[M_{1}(\nu^{T}_{1R}C^{-1}\nu_{1R})+M_{2}(\nu^{T}_{2R}C^{-1}\nu_{3R}+\nu^{T}_{3R}C^{-1}\nu_{2R})\right]$
$\displaystyle+\dfrac{1}{2}\left[\left(y_{\chi_{1}}(\nu^{T}_{3R}C^{-1}\nu_{3R})+y_{\chi_{2}}(\nu^{T}_{1R}C^{-1}\nu_{2R}+\nu^{T}_{2R}C^{-1}\nu_{1R})\right)\chi\right]$

where $y_{k},y_{i}$($k=e,\mu,\tau$; $i=1,2,3$) are Yukawa coupling constants,
$M_{1,2}$ are bare mass terms for right-handed Majorana neutrinos,
$y_{\chi_{1,2}}$ denotes Yukawa coupling constant for interaction terms with
scalar field $\chi$ and $\tilde{\phi}=i\tau_{2}\phi^{*}$; $\tau_{2}$ being
Pauli matrix.
After SSB, charged lepton mass matrix and Dirac mass matrix remains diagonal
as shown in Eqn. (34). But Majorana mass matrix gets modified and takes the
form

$\displaystyle M_{R}=\begin{pmatrix}M_{1}&y_{\chi_{2}}w&0\\\
y_{\chi_{2}}w&0&M_{2}\\\ 0&M_{2}&y_{\chi_{1}}w\\\ \end{pmatrix}.$ (41)

In $M_{D}$-diagonal basis, the inverse neutrino mass matrix is given by

$\displaystyle
M^{-1}_{\nu}=\begin{pmatrix}\dfrac{M_{1}}{v^{2}y^{2}_{1}}&\dfrac{y_{\chi_{2}}w}{v^{2}y_{1}y_{2}}&0\\\
\dfrac{y_{\chi_{2}}w}{v^{2}y_{1}y_{2}}&0&\dfrac{M_{2}}{v^{2}y_{2}y_{3}}\\\
0&\dfrac{M_{2}}{v^{2}y_{2}y_{3}}&\dfrac{y_{\chi_{1}}w}{v^{2}y^{2}_{3}}\\\
\end{pmatrix},$ (42)

which symbolically can be written as

$\displaystyle M^{-1}_{\nu}=\begin{pmatrix}X&X&0\\\ X&0&X\\\ 0&X&X\\\
\end{pmatrix},$ (43)

corresponding to texture $C_{1}$.

##  5 Conclusions

In conclusion, we have investigated the phenomenological implications of two-
zero textures in inverse neutrino mass matrix $M_{\nu}^{-1}$. In the basis
where Dirac neutrino mass matrix $M_{D}$ is diagonal, zeros in right-handed
Majorana neutrino mass matrix $M_{R}$ corresponds to zeros of $M_{\nu}^{-1}$.
We have, also, proposed symmetry realization, based on discrete flavor group
$A_{4}$, wherein such texture zeros can emerge in $M_{\nu}^{-1}$. Further, we
investigate the viability of all two-zero textures in $M_{\nu}^{-1}$ under LMA
and LMA-D solutions. We have categorized the possible textures in five classes
viz. class A, B, C, D and E. Out of fifteen possible two-zero textures of
$M_{\nu}^{-1}$, seven are found to be in consonance with LMA and/or LMA-D
scenario. The general remarks about the obtained phenomenology are as under:

  * •

The textures in class A are all disallowed as they do not reproduce the
correct neutrino phenomenology. Thus, texture with
$\left(M_{\nu}^{-1}\right)_{11}=0$ is disallowed, in general.

  * •

In class B, $B_{2}$ and $B_{4}$ textures are consistent with both LMA and
LMA-D solutions. $B_{2}$ ($B_{4}$) predicts normal (inverted) hierarchical
neutrino masses. For texture $B_{2}$, the $3\sigma$ lower bound on
$0\nu\beta\beta$ decay amplitude $|M_{ee}|$ is found to be $0.02$ eV ($0.04$
eV) under LMA and LMA-D, respectively. For $B_{4}$ texture, it is about $0.06$
eV for both LMA and LMA-D solutions.

  * •

In class C, $C_{2}$ is disallowed. $C_{1}$ and $C_{3}$ textures are consistent
with both LMA and LMA-D solutions. $C_{1}$ ($C_{3}$) predicts normal
(inverted) hierarchical neutrino masses. For texture $C_{1}$ ($C_{3}$), the
$3\sigma$ lower bound on $|M_{ee}|$ is $0.02$ eV ($0.06$) eV under both LMA
and LMA-D solutions.

  * •

For textures $B_{2}$, $B_{4}$, $C_{1}$ and $C_{3}$, the Dirac and Majorana-
type $CP$ violating phases are sharply constrained and these textures are
found to be necessarily $CP$ violating.

  * •

Textures $D_{1}$ and $E_{1}$ predict normal hierarchical neutrino masses and
are found to be consistent with LMA solution. LMA-D solution is disallowed by
these textures. Also, $|M_{ee}|$ is vanishing in these textures. In general,
we conclude that the textures for which LMA-D is disallowed, $|M_{ee}|$ is
vanishing. Similar inference is observed in Ref. [11] wherein the authors
analyzed phenomenology of Majorana neutrino textures in the light of LMA-D
solution.

  * •

Texture $D_{2}$ predicts normal hierarchical neutrino masses and is consistent
with both LMA and LMA-D phenomenology. In LMA (LMA-D) scenario, there exist a
$3\sigma$ lower bound on $|M_{ee}|>0.01$ eV ($0.02$ eV).

  * •

The generic feature of the class of model, discussed in the present work, is
the existence of neutrino mass hierarchy degeneracy in a particular texture.
For example, if a texture is allowed by LMA solution with “X” neutrino mass
hierarchy then, if LMA-D is allowed, it is allowed with the same hierarchy
“X”.

The allowed two-zero texture of $M_{\nu}^{-1}$ viz. $B_{2}$, $B_{4}$, $C_{1}$,
$C_{3}$, $D_{1}$, $D_{2}$ and $E_{1}$ has imperative predictions for
$|M_{ee}|$[40]. Except for $D_{1}$ and $E_{1}$ textures, the predicted
$3\sigma$ lower bound on $0\nu\beta\beta$ decay amplitude $|M_{ee}|$ is
$\mathcal{O}(10^{-2})$ which is within the sensitivity reach of
$0\nu\beta\beta$ decay experiments like SuperNEMO[41], KamLAND-Zen[42],
NEXT[43, 44], and nEXO[45]. For example, the non-observation of
$0\nu\beta\beta$ decay down to these high sensitivities will refute all the
textures except $D_{1}$ and $E_{1}$. Also, we have shown that the allowed
$M_{\nu}^{-1}$ textures can be accommodated in an extension of the SM with
three right-handed neutrinos and one scalar singlet field. As representative
realizations, we have obtained two such textures $B_{4}$ and $C_{1}$ within
Type-I seesaw scenario using $A_{4}$ discrete flavor symmetry.

Acknowledgments
M. K. acknowledges the financial support provided by Department of Science and
Technology(DST), Government of India vide Grant No. DST/INSPIRE
Fellowship/2018/IF180327. The authors, also, acknowledge Department of Physics
and Astronomical Science for providing necessary facility to carry out this
work.

## References

  * [1] S. T. Petcov and M. Piai, Phys. Lett. B 533, 94-106 (2002).
  * [2] L. Wolfenstein, Phys. Rev. D 17, 2369-2374 (1978).
  * [3] P. S. Bhupal Dev et al., SciPost Phys. Proc. 2, 001 (2019).
  * [4] C. Biggio, M. Blennow and E. Fernandez-Martinez, JHEP 08, 090 (2009).
  * [5] P. Coloma et al., JHEP 04, 116 (2017).
  * [6] P. B. Denton, Y. Farzan and I. M. Shoemaker, JHEP 07, 037 (2018).
  * [7] P. B. Denton and J. Gehrlein, [arXiv:2204.09060 [hep-ph]].
  * [8] I. Esteban et al., JHEP 08, 180 (2018).
  * [9] P. Minkowski, Phys. Lett. B 67, 421-428 (1977).
  * [10] R. N. Mohapatra and G. Senjanovic, Phys. Rev. D 23, 165 (1981).
  * [11] H. Borgohain and D. Borah, J. Phys. G 47, no.12, 125002 (2020).
  * [12] M. Ghosh, S. Goswami and A. Mukherjee, Nucl. Phys. B 969, 115460 (2021).
  * [13] Z. H. Zhao, X. Zhang, S. S. Jiang and C. X. Yue, Int. J. Mod. Phys. A 35, no.07, 2050039 (2020).
  * [14] K. N. Deepthi, S. Goswami, V. K.N. and T. K. Poddar, Phys. Rev. D 102, no.01, 015020 (2020).
  * [15] D. Borah, M. Ghosh, S. Gupta and S. K. Raut, Phys. Rev. D 96, no.5, 055017 (2017).
  * [16] E. I. Lashin and N. Chamoun, Phys. Rev. D 85 (2012), 113011.
  * [17] H. Fritzsch and Z. z. Xing, Nucl. Phys. B 556 (1999), 49-75 doi:10.1016/S0550-3213(99)00337-5.
  * [18] A. Merle and W. Rodejohann, Phys. Rev. D 73 (2006), 073012.
  * [19] S. Dev and S. Kumar, Mod. Phys. Lett. A 22 (2007), 1401-1410.
  * [20] H. Fritzsch, Z. z. Xing and S. Zhou, JHEP 09 (2011), 083.
  * [21] P. H. Frampton, S. L. Glashow and D. Marfatia, Phys. Lett. B 536 (2002), 79-82.
  * [22] B. R. Desai, D. P. Roy and A. R. Vaucher, Mod. Phys. Lett. A 18 (2003), 1355-1366.
  * [23] Z. z. Xing, Phys. Lett. B 530 (2002), 159-166.
  * [24] W. l. Guo and Z. z. Xing, Phys. Rev. D 67 (2003), 053002.
  * [25] M. Singh, G. Ahuja and M. Gupta, PTEP 2016 (2016) no.12, 123B08.
  * [26] L. Lavoura, Phys. Lett. B 609, 317-322 (2005).
  * [27] Z. H. Zhao, X. Zhang, S. S. Jiang and C. X. Yue, Int. J. Mod. Phys. A 35, no.07, 2050039 (2020).
  * [28] S. Verma, Nucl. Phys. B 854, 340-349 (2012).
  * [29] S. Choubey and D. Pramanik, JHEP 12, 133 (2020).
  * [30] P. A. Zyla et al. [Particle Data Group], PTEP 2020, no.08, 083C01 (2020).
  * [31] P. I. Krastev and S. T. Petcov, Phys. Lett. B 205, 84-92 (1988).
  * [32] C. Jarlskog, Phys. Rev. Lett. 55, 1039 (1985)
  * [33] O. G. Miranda, M. A. Tortola and J. W. F. Valle, AIP Conf. Proc. 917, no.01, 100-107 (2007).
  * [34] K. N. Vishnudath, S. Choubey and S. Goswami, Phys. Rev. D 99, no.09, 095038 (2019).
  * [35] F. J. Escrihuela, O. G. Miranda, M. A. Tortola and J. W. F. Valle, Phys. Rev. D 80, 105009 (2009).
  * [36] Y. Farzan and M. Tortola, Front. in Phys. 6, 10 (2018).
  * [37] P. F. de Salas et al., JHEP 02, 071 (2021).
  * [38] H. Ishimori et al., Lect. Notes Phys. 858, 1-227 (2012).
  * [39] E. Ma and G. Rajasekaran, Phys. Rev. D 64, 113012 (2001).
  * [40] K. N. Vishnudath, S. Choubey and S. Goswami, Phys. Rev. D 99, no.09, 095038 (2019).
  * [41] A. S. Barabash, J. Phys. Conf. Ser. 375, 042012 (2012).
  * [42] A. Gando et al., Phys. Rev. Lett. 117, 082503 (2016).
  * [43] F. Granena et al., arXiv:0907.4054 [hep-ex].
  * [44] J. J. Gomez-Cadenas et al., Adv. High Energy Phys. 2014, 907067 (2014).
  * [45] C. Licciardi, J. Phys. Conf. Ser. 888, 012237 (2017).