Neutrino mass model at a three-loop level from a non-holomorphic modular A₄ symmetry

The successful construction of the non-holomorphic modular symmetry framework by Qu and Ding [1] has enabled us to safely handle the beyond the standard model (BSM) without super-symmetric theories when using the framework for a flavor symmetry. In fact, the non-holomorphic symmetries have been applied to some non-supersymmetric models [2−12] in order to restrict the number of model parameters. In constructing a model, we have the advantage of applying non-supersymmetric framework to reduce the number of new fields when extra fields are required to cancel a gauge anomaly in the supersymmetric case.

Radiatively induced neutrino mass models are representative scenarios that do not require the super-symmetric framework, and new particles can be connected to the standard model particles. Sometimes, the model can possess a dark matter (DM) candidate [13] that often requires an additional symmetry to stabilize it. Thus, constructing radiative neutrino mass models (with DM) using the non-holomorphic modular symmetry can make a model more attractive by realizing more predictability.

In this study, we apply a non-holomorphic $ A_4 $ flavor symmetry to a well-known three-loop neutrino mass model [14]. The three-loop neutrino model is phenomenologically interesting, as the scale of new particles would be smaller compared to lower loop (or tree) level models owing to loop suppression. We then expect rich phenomenology such as collider and lepton flavor physics. The non-holomorphic modular symmetry framework is suitable for constructing such a three-loop model in a minimal manner; if we consider a holomorphic framework, we need to add more fields to cancel the gauge anomaly. Then, we need to determine the minimal number of free parameters to fit the observables in the lepton sector under the non-holomorphic modular $ A_4 $ symmetry. Through chi-square numerical analysis, we search for the minimum model to predict the lepton masses and mixing angles in addition to reproducing the current neutrino observables in Nufit 6.0 [15]. Then, we perform further numerical analyses to satisfy lepton flavor violations (LFVs), the muon anomalous magnetic moment, (muon $ g-2 $), lepton universality, and DM. The results show that relic density is too large within the limit of perturbation, thereby requiring a new interaction that adds one singlet scalar boson without changing predictions in the neutrino sector.

This paper is organized as follows. In Sec. II, we explain our minimum three-loop neutrino mass model and construct the renormalizable Lagrangian in the lepton sector, Higgs sector, charged-lepton sector, heavier Majorana fermion sector, and active-neutrino sector. Then, we formulate the LFVs, muon $ g-2 $, lepton universality, and relic density of the DM. In Sec. III, we perform χ square analysis and present predictions for normal and inverted hierarchies in the neutrino sector. By employing the benchmark points of the best-fit values in the lepton sector, we further demonstrate the numerical analyses for the LFVs, muon $ g-2 $, lepton universality, and relic density of the DM. We present the conclusions and discussion in Sec. IV. In Appendix A, we show the three-loop function in the neutrino sector.

In this section, we show the setup of the model based on a $ G_{\rm{SM}} \times A_4 $ symmetry, where $ G_{\rm{SM}} $ is the SM gauge symmetry and $ A_4 $ is the modular one. In the lepton sector, we introduce a singlet fermion, which is a triplet under $ A_4 $ with modular weight 0. In the scalar sector, we introduce two charged singlets distinguished by modular weights +2 and –1. The SM leptons $ \overline{L_L} $ and $ \ell_R $ are also $ A_4 $ triplets with modular weights –1 and +1, respectively. The assignments are summarized in Table 1. By assigning modular weights, we can eliminate unwanted terms such as $ \overline{N_R} L_L H $, and the neutrino masses are generated at the three-loop level, as discussed below.

The relevant Lagrangian under these symmetries is given by

where we define $ Y_3^{(0)} = [y_1,y_2,y_3] $ [1], and ''$ \cdot $'' indicates the $ {\rm i} \sigma_2 $ factor that makes the term $ S U(2)_L $ invariant. The first two terms generate the mass of the charged leptons, and parameters $ \{a_e,a_\mu,a_\tau \} $ are real without loss of generality and are rephased into $ e_R,\mu_R,\tau_R $, respectively.

In this section, we show setup of the model based on $ G_{\rm{SM}} \times A_4 $ symmetry where $ G_{\rm{SM}} $ being the SM gauge symmetry and $ A_4 $ is modular one. In lepton sector, We introduce singlet fermion which is triplet under $ A_4 $ with modular weight 0. In scalar sector, we introduce two charged singlets distinguished by modular weight +2 and -1. The SM leptons $ \overline{L_L} $ and $ \ell_R $ are also $ A_4 $ triplets with modular weight -1 and +1 respectively. The assignments are summarized in Table 1. By the assignments of modular weight, we can eliminate unwanted terms like $ \overline{N_R} L_L H $ and neutrino masses are generated at three-loop level as discussed below.

	Leptons				Bosons
	$ \overline{L_L}$	${\ell}_R$	$N_R$		$H$	$S^+_1$	$S^+_2$
$SU(2)_L$	2	1	1		2	1	1
$U(1)_Y$	$-\frac12$	1	0		$\frac12$	+1	+1
$A_4$	3	$\{ 1\}$	3		1	1	1
$-k_I$	-1	+1	0		0	+2	-1

Table 1.  Field contents and stheir charge assignments in the model under $SU(2)_L\times U(1)_Y\times A_4 $ where $-k_I $ is the number of modular weight. Here $\{1\} $ stands for the combination of $ A_4$ singlets $\{1,1',1''\} $.

The relevant Lagrangian under these symmetries is given by

where we define $ Y_3^{(0)} = [y_1,y_2,y_3] $ [1] and $ "\cdot" $ indicates $ i \sigma_2 $ factor to make the term $ SU(2)_L $ invariant. The first two terms generates the mass of charged-leptons and parameters $ \{a_e,a_\mu,a_\tau \} $ are real without loss of generality by rephasing them into $ e_R,\mu_R,\tau_R $, respectively.

The scalar potential in the model is given by

The SM Higgs field is denoted by

and $ v\approx 246 $ GeV is the vacuum expectation value (VEV) in the Higgs basis after the spontaneous symmetry breaking, z is absorbed by the neutral gauge boson of the SM Z, and $ w^+ $ is absorbed by the charged gauge boson of the SM $ W^+ $. The charged scalar masses are given by

In the numerical analysis, we consider $ m_{S_{1,2}} $ to be free parameters.

The scalar potential in the model is given by

The SM Higgs field is denoted by

and $ v\approx 246 $ GeV is vacuum expectation value (VEV) in the Higgs basis after the spontaneous symmetry breaking, z is absorbed by the neutral gauge boson of the SM Z, and $ w^+ $ is absorbed by the charged gauge boson of the SM $ W^+ $. The charged scalar masses are respectively given by

In the numerical analysis we consider $ m_{S_{1,2}} $ to be free parameters.

After the spontaneous electroweak symmetry breaking, the charged-lepton mass matrix $ M_e $ is given by

Then, the charged-lepton mass matrix is diagonalized by a bi-unitary mixing matrix as $ D_{\ell} \equiv{\rm{diag}}(m_e,m_\mu,m_\tau)= V^{\dagger}_{eL} M_e V_{eR} $. Therefore, $ \ell_{L(R)}\equiv V_{eL(R)} \ell'_{L(R)} $ where $ \ell'_{L(R)} $ is the mass eigenstate. These three parameters are used in order to fit the mass eigenvalues of charged-leptons by solving the following three relations:

For our convenience to construct the neutrino mass matrix below, we define $ \tilde D_\ell $ that is given by $ D_\ell\equiv m_\tau \tilde D_\ell $.

After the spontaneous electroweak symmetry breaking, the charged-lepton mass matrix $ M_e $ is given by

Then, the charged-lepton mass matrix is diagonalized by a bi-unitary mixing matrix as $ D_{\ell} \equiv{\rm{diag}}(m_e,m_\mu,m_\tau)= V^{\dagger}_{eL} M_e V_{eR} $. Therefore, $ \ell_{L(R)}\equiv V_{eL(R)} \ell'_{L(R)} $, where $ \ell'_{L(R)} $ is the mass eigenstate. These three parameters are used to fit the mass eigenvalues of charged-leptons by solving the following three relations:

For convenience in constructing the neutrino-mass matrix, we define $ \tilde D_\ell $, which is given by $ D_\ell\equiv m_\tau \tilde D_\ell $.

The heavier Majorana mass matrix is given by

where $ \tilde M_2\equiv M_2/M_1 $ can be real without loss of generality. $ M_N $ is diagonalized by $ D_N\equiv U_N^T M_N U_N $ ($ \tilde D_N \equiv U_N^T \tilde M_N U_N $)., therefore $ N_R\equiv U_N \psi_R $. Here, $ \psi_R $ is the mass eigenstate.

The heavier Majorana mass matrix is given by

where $ \tilde M_2\equiv M_2/M_1 $ can be real without loss of generality. $ M_N $ is diagonalized by $ D_N\equiv U_N^T M_N U_N $ ($ \tilde D_N \equiv U_N^T \tilde M_N U_N $); therefore, $ N_R\equiv U_N \psi_R $. Here, $ \psi_R $ is the mass eigenstate.

The active neutrino mass matrix is given at the three-loop level via the following Lagrangian in terms of the mass eigenstates:

where $ H\equiv h V^*_{eL} $ and $ Y\equiv V^T_{eR} y U_N $. The Yukawa matrices y and h are as follows:

where $ \tilde c(\tilde d)_\nu\equiv c(d)_\nu/b_\nu $ are complex free parameters. The neutrino mass matrix is then given by

Here, F is a loop function via three loop diagram and it depends on the mass eigenvalues of $ \{\psi_R, S_1^+, S_2^+\} $. ¹ Since the masses of $ \psi_R $ contribute to the structure of neutrino mass matrix, there would be too many free parameters to get some predictions for the neutrino sector. Thus, we consider a special situation among the mass hierarchies of $ \psi_R, S_1^+, S_2^+ $ so that F is independent of the structure of neutrino mass matrix. When we assume $ D_{N_i} \ll m_{S_1}\sim m_{S_2} $, one finds that the dominant part of the loop-function F is a constant and can explicitly be given by $ F\approx 0.062 $. In detail, one finds Appendix A. Thus, we redefine the neutrino mass matrix as follows:

The dimensionless matrix $ \tilde m_\nu $ is diagonalized by a unitary matrix $ U_\nu $ as $ U_\nu^T \tilde m_\nu U_\nu =\tilde D_\nu $, where $ \tilde D_\nu = {\rm{diag}}[\tilde D_{\nu_1}, \; \tilde D_{\nu_2}, \tilde D_{\nu_3}] $, and the Pontecorvo-Maki-Nakagawa-Sakata unitary matrix $ U_{\rm{PMNS}} $ is defined by $ V_{eL}^{\dagger} U_\nu $. Note here that the lightest neutrino mass is zero due to the two-matrix rank of the neutrino. Thus, the atmospheric mass squared difference $ \Delta m^2_{\rm{atm}} $ is as follows:

where NH(IH) represents the normal(inverted) hierarchy. The solar mass squared difference $ \Delta m^2_{\rm{sol}} $ is given by

The effective mass for neutrinoless double beta decay is given by

where the Majorana phase is defined by ${\rm{diag}}[1, {\rm e}^{{\rm i}\alpha_{21}/2} ,1]$ and we adopt the standard parametrization for the PMNS unitary matrix. The current KamLAND-Zen data [16] provide measured observables, and their upper bound is given by $ \langle m_{ee}\rangle<(28-122) $ meV at a 90% confidence level. The minimal cosmological model ΛCDM $ +\sum D_{\nu} $ provides an upper bound on $ \sum D_{\nu}\le $ 120 meV [17, 18]. Moreover, the recent combination of DESI and CMB data gives a more stringent upper bound on this bound; $ \sum D_{\nu}\le $ 72 meV [19].

The active neutrino mass matrix is given at three-loop level via the following Lagrangian in terms of mass eigenstates

where $ H\equiv h V^*_{eL} $ and $ Y\equiv V^T_{eR} y U_N $. The Yukawa matrices y and h are respectively found as

where $ \tilde c(\tilde d)_\nu\equiv c(d)_\nu/b_\nu $ are complex free parameters. The neutrino mass matrix is then given by

Here, F is a loop function via three loop diagram and it depends on the mass eigenvalues of $ \{\psi_R, S_1^+, S_2^+\} $. ¹ Since the masses of $ \psi_R $ contribute to the structure of neutrino mass matrix, there would be too many free parameters to get some predictions for the neutrino sector. Thus, we consider a special situation among the mass hierarchies of $ \psi_R, S_1^+, S_2^+ $ so that F is independent of the structure of neutrino mass matrix. When we assume $ D_{N_i} \ll m_{S_1}\sim m_{S_2} $, one finds that the dominant part of the loop-function F is a constant and can explicitly be given by $ F\approx 0.062 $. In detail, one finds Appendix A. Thus, we redefine the neutrino mass matrix as follows:

The dimensionless matrix $ \tilde m_\nu $ is diagonalized by a unitary matrix $ U_\nu $ as $ U_\nu^T \tilde m_\nu U_\nu =\tilde D_\nu $, where $ \tilde D_\nu = {\rm{diag}}[\tilde D_{\nu_1},\tilde D_{\nu_2},\tilde D_{\nu_3}] $ and the Pontecorvo-Maki-Nakagawa-Sakata unitary matrix $ U_{\rm{PMNS}} $ is defined by $ V_{eL}^{\dagger} U_\nu $. Note here that the lightest neutrino mass is zero due to two matrix rank of the neutrino. The atmospheric mass squared difference $ \Delta m^2_{\rm{atm}} $ is thus found as

where NH(IH) represents normal(inverted) hierarchy. The solar mass squared difference $ \Delta m^2_{\rm{sol}} $ is given by

The effective mass for neutrinoless double beta decay is given by

where Majorana phase is defined by $ {\rm{diag}}[1,e^{i\alpha_{21}/2} ,1] $ and we adopt the standard parametrization for the PMNS unitary matrix. A current KamLAND-Zen data [16]. provides measured observable in future and its upper bound is given by $ \langle m_{ee}\rangle<(28-122) $ meV at 90% confidence level. The minimal cosmological model ΛCDM $ +\sum D_{\nu} $ provides upper bound on $ \sum D_{\nu}\le $ 120 meV [17, 18]. Moreover, recently combination of DESI and CMB data gives more stringent upper bound on this bound; $ \sum D_{\nu}\le $ 72 meV [19].

$ \ell_\alpha \to \ell_\beta \gamma $ process: First of all, let us consider the processes $ \ell_\alpha \to \ell_\beta \gamma $ at one-loop level ². The formula for the branching ratio can generally be written as

where $ \alpha_{\rm{em}}\approx1/137 $ is the fine-structure constant, $ C_{\alpha\beta} \approx(1,0.1784, 0.1736) $ for ($ (\alpha,\beta)=(\mu,e),(\tau,e),(\tau,\mu) $), $ {\rm{G_F}}\approx1.17\times 10^{-5} $ GeV^-2 is the Fermi constant, and $ a_{L/R} $ is respectively given by

where

Once we assume that $ m_{\ell_\alpha} \gg m_{\ell_\beta} $, the formula can be simplified to

The formula for the muon $ g-2 $ can be written in terms of $ a_L $ and $ a_R $, and simplified as follows:

Notice here that this contribution to the muon $ g-2 $ is negative, yet it is negligible compared to the deviation in the experimental value $ {\cal O}(10^{-9}) $ [22].

$ \ell_\alpha \to \ell_\beta \gamma $ process: First of all, let us consider the processes $ \ell_\alpha \to \ell_\beta \gamma $ at one-loop level ². The formula for the branching ratio can generally be written as

where $ \alpha_{\rm{em}}\approx 1/137 $ is the fine-structure constant, $ C_{\alpha\beta} \approx (1,0.1784, \; 0.1736) $ for ($ (\alpha,\beta)=(\mu,e),\; (\tau,e), \; (\tau,\mu) $), $ G_{\rm F} \approx 1.17\times 10^{-5} $ GeV^-2 is the Fermi constant, and $ a_{L/R} $ is given by

where

By assuming that $ m_{\ell_\alpha} \gg m_{\ell_\beta} $, the formula can be simplified to

The formula for the muon $ g-2 $ can be written in terms of $ a_L $ and $ a_R $ and simplified as follows:

Notice here that this contribution to the muon $ g-2 $ is negative; however, it is negligible compared to the deviation in the experimental value $ {\cal O}(10^{-9}) $ [22].

Here, we just employ the results of lepton universality from precursor's works [23] whose results provide us the upper bounds on coupling H in terms of $ m_{S_1} $ and $ a_\nu $. We summarize these results in Table 3.

Process	Experiments	Bound (90% CL)
$ {\rm{Lepton/hadron\ universality}} $	$ \sum_{q=b,s,d}|V^{\rm{exp}}_{uq}|^2=0.9999\pm0.0006 $:	$ |H^{\dagger}_{e\mu}|^2<0.007\left(\frac{m_{S_1}}{a_\nu {\rm{TeV}}}\right)^2 $
$ {\rm{\mu/e\ universality}} $	$ \frac{G_\mu^{\rm{exp}}}{G_e^{\rm{exp}}}=1.0010\pm0.0009 $	$ ||H^{\dagger}_{\mu\tau}|^2-|H^{\dagger}_{e\tau}|^2|<0.024\left(\frac{m_{S_1}}{a_\nu {\rm{TeV}}}\right)^2 $
$ {\rm{\tau/\mu\ universality}} $	$ \frac{G_\tau^{\rm{exp}}}{G_\mu^{\rm{exp}}}=0.9998\pm0.0013 $	$ ||H^{\dagger}_{e\tau}|^2-|H^{\dagger}_{e\mu}|^2|<0.035\left(\frac{m_{S_1}}{a_\nu {\rm{TeV}}}\right)^2 $
$ {\rm{\tau/e\ universality}} $	$ \frac{G_\tau^{\rm{exp}}}{G_e^{\rm{exp}}}=1.0034\pm0.0015 $	$ ||H^{\dagger}_{\mu\tau}|^2-|H^{\dagger}_{e\mu}|^2|<0.04\left(\frac{m_{S_1}}{a_\nu {\rm{TeV}}}\right)^2 $

Table 3.  Summary of the lepton universality and the corresponding bounds on $ f_{\alpha\beta} $.

Process	$ (\alpha,\beta) $	Experimental bounds (90% CL)	References
$ \mu^{-} \to e^{-} \gamma $	$ (\mu,e) $	$ {BR}(\mu \to e\gamma)< 4.2 \times 10^{-13} $	[20]
$ \tau^{-} \to e^{-} \gamma $	$ (\tau,e) $	$ {BR}(\tau \to e\gamma)< 3.3 \times 10^{-8} $	[21]
$ \tau^{-} \to \mu^{-} \gamma $	$ (\tau,\mu) $	$ {BR}(\tau \to \mu\gamma)< 4.4 \times 10^{-8} $	[21]

Table 2.  Summary for the experimental bounds of the LFV processes $ \ell_\alpha \to \ell_\beta \gamma $.

Here, we employ only the results of lepton universality from a precursor work [23]; the results provide the upper bounds on coupling H in terms of $ m_{S_1} $ and $ a_\nu $. We summarize these results in Table 3.

Relic density: Our DM is identified as the lightest Majorana fermion $ N_1 $ where we denote $ N_1 $ as X hereafter and its mass is $ m_\chi $. In order to analyze it simpler, we impose the following condition, $ 1.2 m_\chi\lesssim D_{N_2}\le D_{N_3} $, in order to evade an effect of co-annihilation interactions for the relic density of DM. ³ Under the condition, the dominant contribution to the relic density arises from Y. Then, the non-relativistic cross section is expanded by relative velocity $ v_{\rm{rel}}^2 $; $ (\sigma v_{\rm{rel}})\approx a_{\rm{eff}} + b_{\rm{eff}} v^2_{\rm{rel}}+ {\cal O}(v^4_{\rm{rel}}) $ and found as follows:

where we have neglected the masses of charged-leptons. The above cross section suggests that it is p-wave dominant. The relic density is then given by

where $ g^*\approx100 $, $ M_P\approx 1.22\times 10^{19}{\rm{GeV}} $, $ x_f\approx20 $. In our numerical analysis below, we use rather relaxed experimental range $ 0.11\le \Omega h^2\le 0.13 $, since we simplify our analysis of the relic density.

Relic density: Our DM is identified as the lightest Majorana fermion $ N_1 $ where we denote $ N_1 $ as X hereafter and its mass is $ m_\chi $. In order to analyze it simpler, we impose the following condition, $ 1.2 m_\chi\lesssim D_{N_2}\le D_{N_3} $, in order to evade an effect of co-annihilation interactions for the relic density of DM. ³ Under the condition, the dominant contribution to the relic density arises from Y. Then, the non-relativistic cross section is expanded by relative velocity $ v_{\rm{rel}}^2 $; $ (\sigma v_{\rm{rel}})\approx a_{\rm{eff}} + b_{\rm{eff}} v^2_{\rm{rel}}+ {\cal O}(v^4_{\rm{rel}}) $ and found as follows:

where we have neglected the masses of charged leptons. The above cross section suggests that it is p-wave dominant. The relic density is then given by

where $ g^*\approx100 $, $ M_P\approx 1.22\times 10^{19}{\rm{GeV}} $, $ x_f\approx20 $. In our numerical analysis, we use a rather relaxed experimental range $ 0.11\le \Omega h^2\le 0.13 $ because we simplify our analysis of the relic density.

In this section, we demonstrate numerical analyses based on all the experimental results that we have discussed. Then, we show the results of the LFVs, lepton $ g-2 $, and DM.

In this section, we demonstrate numerical analyses in light of all the experimental results which we discussed before. Then, we show the results of LFVs, lepton $ g-2 $, and DM.

First, we perform a χ square analysis adopting data from NuFit6.0 [15], where we use five reliable observables (three mixings, two mass square differences) for the analysis. The yellow points represent the interval of $ 2\sigma-3\sigma $, and the red ones $ 3\sigma-5\sigma $, where no solutions are obtained within $ 2\sigma $. Our three input parameters are randomly selected within the following range:

where we work on the fundamental region of τ, and $ \tilde c_\nu, \tilde d_\nu $ are complex.

After the numerical analysis, we find that the IH case is not favored in the model, where the minimal $ \chi^2 $ can be at most $ {\cal{O}}(1500) $. Thus, we summarize our results using only the NH case in the next subsection. Note that the parameters $ \{a_e, a_\mu, a_\tau \} $ are chosen to fit the observed charged-lepton masses, and $ \{a_\nu, b_\nu, M_1\} $ are related to fix the scale of the neutrino mass via κ, defined in Eq. (15). Thus, the relative neutrino mass and three mixing angles are fitted using the remaining parameters $ \{\tau, c_\nu, d_\nu, M_2\} $ corresponding to seven real parameters. Three of these real parameters are related to complex phases; therefore, fitting the neutrino data is nontriavial. In fact, we would not be able to obtain any solutions in the IH case. To improve the fitting further, such as for IH, we need to change the assignment of the modular weight to increase the number of free parameters.

At first, we perform χ square analysis adopting data from NuFit6.0 [15], where we have adopt five reliable observables; three mixings, two mass square differences, for the analysis. The yellow points represents the interval of $ 2\sigma-3\sigma $, and red one $ 3\sigma-5\sigma $, where we would not find any solutions within $ 2\sigma $. Our three input parameters randomly select within the following range:

where we work on the fundamental region of τ and $ \tilde c_\nu, \tilde d_\nu $ are complex.

After the numerical analysis, we find that IH case is not favored in the model where minimal $ \chi^2 $ can be at most $ {\cal{O}}(1500) $. We thus summarize our results regarding only NH case in next subsection. Note that the parameters $ \{a_e, a_\mu, a_\tau \} $ are chosen to fit the observed charged lepton masses and $ \{a_\nu, b_\nu, M_1\} $ are related to fix the scale of neutrino mass via κ defined in Eq. (15). Thus relative neutrino mass and three mixing angles are fitted using remaining parameters $ \{\tau, c_\nu, d_\nu, M_2\} $ corresponding to 7 real parameters. Since three of these real parameters are related to complex phases it is not trivial if we can fit the neutrino data. In fact, we would not be able to obtain any solutions in IH case. To improve the fitting further such as IH, we need to change assignment of modular weight to increase the number of free parameter.

In Fig. 1, we show the allowed region of τ, and find that the allowed region is concentrated at nearby $ |{\rm{Re}}[\tau]|=[0.0-0.2] $ and $ {\rm{Im}}[\tau]=[1.26-1.28] $ where the value is close to the fixed point $ \tau =i $. We also find a few points near the fixed point $ \tau = \omega $. ⁴

Figure 1.  (color online) Allowed region for real τ and imaginary τ in NH.

In Fig. 2, we demonstrate the allowed regions for the absolute values (left) and argument ones (right) of $ \tilde d_\nu $ and $ \tilde c_\nu $ in NH. We show that the allowed region is at approximately $ |\tilde c_\nu|=[10^{-4}-10^{5}] $ and $ |\tilde d_\nu|=[10^{-5}-10^{4}] $, where $ |\tilde{d}_\nu | \ll |\tilde{c}_\nu | $ is preferred, and $ {\rm{Arg}}[\tilde c_\nu] $ and $ {\rm{Arg}}[\tilde d_\nu] $ can be any value with little correlation.

Figure 2.  (color online) Allowed regions for absolute values (left) and argument ones (right) of $ \tilde d_\nu $ and $ \tilde c_\nu $ in NH.

In Fig. 3, we display the allowed region for $\delta_{\rm CP}$ deg (left) and $ \langle m_{ee}\rangle $ meV (right) in terms of $ \sum D_\nu $ meV. We show that most of the points are located at $|\delta_{\rm CP}| = [90-200]$ deg and few points are at approximately $\delta_{\rm CP} = [40-60]$ deg, and $ \langle m_{ee}\rangle \approx [1-4] $ meV. The vertical magenta dotted line is the upper bound of the results of Planck+DESI [19] $ \sum D_\nu\le $72 meV, and the $ \sum D_\nu $ range of our model is $ [58-60] $ meV, which is a trivial consequence of two nonzero mass eigenvalues of active neutrinos.

Figure 3.  (color online) Allowed regions for $\delta_{\rm CP}$ deg (left) and $ \langle m_{ee}\rangle $ meV (right) in terms of $ \sum D_\nu $ meV in NH. The vertical magenta dotted line is the upper bound of the results of Planck+DESI [19], $ \sum D_\nu\le $72 meV. The cyan region in the left panel indicates the allowed region based on the experimental results of Nufit 6.0.

In Fig. 4, we show the allowed region for $ \langle m_{ee}\rangle $ meV (left) and $ \alpha_{21} $ deg (right) in terms of $\delta_{\rm CP}$ deg in NH. We show that the allowed region of $ \alpha_{21} $ is concentrated at approximately $ [80-270] $ deg with few points outside the region.

Figure 4.  (color online) Allowed region for $ \langle m_{ee}\rangle $ meV (left) and $\delta_{\rm CP}$ deg (right) in terms of $\delta_{\rm CP}$ deg in NH.

We show a benchmark point (BP) that has the minimum $ \Delta \chi^2 $ in Table 4 and this BP will be employed to analyze the LFV, $ g-2 $, and DM in the next subsection.

In Fig. 1, we show the allowed region of τ, and find that the allowed region is concentrated at nearby $ |{\rm{Re}}[\tau]|=[0.0-0.2] $ and $ {\rm{Im}}[\tau]=[1.26-1.28] $ where the value is close to the fixed point $ \tau =i $. We also find a few points near the fixed point $ \tau = \omega $. ⁴

Figure 1.  (color online) Allowed region for real τ and imaginary τ in NH.

In Fig. 2, we demonstrate the allowed regions for absolute values (left) and argument ones (right) of $ \tilde d_\nu $ and $ \tilde c_\nu $ in NH. We find that the allowed region is about $ |\tilde c_\nu|=[10^{-4}-10^{5}] $ and $ |\tilde d_\nu|=[10^{-5}-10^{4}] $ where $ |\tilde{d}_\nu | \ll |\tilde{c}_\nu | $ is preferred, and $ {\rm{Arg}}[\tilde c_\nu] $ and $ {\rm{Arg}}[\tilde d_\nu] $ can be any value with little correlation.

Figure 2.  (color online) Allowed regions for absolute values (left) and argument ones (right) of $ \tilde d_\nu $ and $ \tilde c_\nu $ in NH.

In Fig. 3, we display the allowed region for $ \delta_{CP} $ deg (left) and $ \langle m_{ee}\rangle $ meV (right) in terms of $ \sum D_\nu $ meV. We find the most of points are located at $ |\delta_{CP}| = [90-200] $ deg and few points are around $ \delta_{CP} = [40-60] $ deg. $ \langle m_{ee}\rangle \approx [1-4] $ meV. The vertical magenta dotted line is upper bound on results of Planck+DESI [19] $ \sum D_\nu\le $72 meV while $ \sum D_\nu $ of our model is $ [58-60] $ meV which is nothing but a trivial consequence of two nonzero mass eigenvalues of active neutrinos.

Figure 3.  (color online) Allowed regions for $ \delta_{CP} $ deg (left) and $ \langle m_{ee}\rangle $ meV (right) in terms of $ \sum D_\nu $ meV in NH. The vertical magenta dotted line is upper bound on results of Planck+DESI [19] $ \sum D_\nu\le $72 meV. The cyan region in the left panel indicates allowed region by experimental result of Nufit 6.0.

In Fig. 4, we show the allowed region for $ \langle m_{ee}\rangle $ meV (left) and $ \alpha_{21} $ deg (right) in terms of $ \delta_{CP} $ deg in NH. We find the allowed region of $ \alpha_{21} $ is concentrated around $ [80-270] $ deg with few points outside the region.

Figure 4.  (color online) Allowed region for $ \langle m_{ee}\rangle $ meV (left) and $ \delta_{CP} $ deg (right) in terms of $ \delta_{CP} $ deg in NH.

We show a benchmark point (BP) that has the minimum $ \Delta \chi^2 $ in Table 4 and this BP will be employed to analyze the LFV, $ g-2 $, and DM in the next subsection.

	NH
τ	$ 0.137 + 1.26 i $
$ \tilde M_2 $	$ 5.34 \times 10^{-4} $
$ \tilde c_\nu $	$ -2.85 \times 10^{3} - 1.69 \times 10^3 i $
$ \tilde d_\nu $	$ -26.3 + 64.5 i $
$ [a_e,a_\mu,a_\tau] $	$ [7.21 \times10^{-6}, -0.00139, 0.0206] $
$ \Delta m^2_{\rm{atm}} $	$ 2.51 \times10^{-3} {\rm{eV}}^2 $
$ \Delta m^2_{\rm{sol}} $	$ 7.56 \times10^{-5} {\rm{eV}}^2 $
$ \sin\theta_{12} $	0.553
$ \sin\theta_{23} $	0.683
$ \sin\theta_{13} $	0.147
$ [\delta_{\rm{CP}}^\ell,\ \alpha_{21}] $	$ [-170^\circ,\, 257^\circ] $
$ \sum m_i $	58.8 meV
$ \langle m_{ee} \rangle $	2.94 meV
κ	$ 3.49\times 10^{-14} $
$ \sqrt{\Delta\chi^2} $	2.24

Table 4.  Numerical benchmark point (BP) of our input parameters and observables in NH. Here, this BP is taken $ \sqrt{\Delta \chi^2} $ to be minimum.

Before our numerical analysis, we present some definitions. The neutrino-mass matrix does not depend on all the masses inside the loop, but the chi square analysis of the neutrino-oscillation data provides the value of κ. Their masses inside the loop determine the values of the LFVs, muon $ g-2 $, and relic density of DM. Thus, we rewrite Eq. (15) as follows:

When $ a_\nu, b_\nu $, and $ M_1 $ are numerically fixed, $ \lambda_0 $ is numerically determined. Then, we impose the perturbative limit in our numerical analysis to be

In addition, we restrict ourselves to the following conditions to forbid co-annihilation processes and obtain the mass-independent loop function of the neutrino-mass matrix:

where we have defined $ \epsilon_3 $ to be $ \dfrac{D_{N_3}}{m_{S_1}} $.

Our input parameters are randomly selected from the following range:

where $ a_\nu, b_\nu $ are real and the other required parameters are employed by the BP in the previous section.

Our numerical analysis showed that Yukawa coupling $ |b_\nu\times Y| $ exceeds the perturbative limit $ \sim4\pi $ to obtain the observed relic density of DM while satisfying the constraints of LFVs and lepton universalities. The correct relic density requires $ {\cal{O}}(100) \lesssim{\rm{Max}}[|b_\nu\times Y|] $ for the NH case, applying allowed parameters that can fit the neutrino data. This implies that co-annihilations do not help to reduce the Yukawa couplings to the perturbative limit. We may move to one of the next minimum models by changing the modular weight of $ N_R $ to –2 instead of 0 to obtain one more mass parameter. This provides a wider region of allowed parameters, where the other assignments are the same as our model. However, we would still encounter difficulty in realizing the correct relic density while keeping the perturbative limit for the Yukawa couplings. This is because the DM annihilation cross section, Eq. (29), is p-wave dominant and we need a relatively larger coupling constant than that of the s-wave case. In addition, neutrino data and LFV constraints require heavy DM and new scalars that also suppress the DM annihilation cross section. Thus, obtaining the correct relic density in our minimal setting is difficult, and some extension is necessary.

If we do not satisfy the observed relic density and we perform our numerical analysis under the perturbative limit, we obtain the tendencies for electron $ g-2 $, muon $ g-2 $, and LFVs, as shown in Fig. 5. These figures suggests that $ -\Delta a_e $ and $ -\Delta a_\mu $ are at most $ 10^{-20} $ and $ 10^{-15} $, respectively. However, LFVs, especially the $ \mu\to e\gamma $ branching ratio, would be testable in the near future because its maximum value is close to the experimental limit.

Figure 5.  (color online) Allowed region for electron $ g-2 $ (left), muon $ g-2 $ (center), and LFVs (right), where these points do not satisfy the observed relic density.

Before our numerical analysis, we prepare some definitions. The neutrino mass matrix does not depend on all the masses inside the loop, but the chi square analysis of the neutrino oscillation data gives us the value of κ. While their masses inside the loop determine the values of LFVs, muon $ g-2 $, and the relic density of DM. Thus, we rewrite Eq. (15) as follows:

When $ a_\nu, b_\nu $, and $ M_1 $ are numerically fixed, $ \lambda_0 $ is numerically determined. Then we impose the perturbative limit in our numerical analysis to be

In addition, we restrict ourselves to be following conditions in order to forbid co-annihilation processes and obtain the mass-independent loop function of the neutrino mass matrix:

where we have defined $ \epsilon_3 $ to be $ \frac{D_{N_3}}{m_{S_1}} $.

Our input parameters randomly select within the following range:

where $ a_\nu, b_\nu $ are real and the other needed parameters are employed by BP in the previous section.

In our numerical analysis, we found that Yukawa coupling $ |b_\nu\times Y| $ exceeds the perturbative limit $ \sim4\pi $ to obtain the observed relic density of DM while satisfying the constraints of LFVs and lepton universalities. The correct relic density requires $ {\cal{O}}(100) \lesssim{\rm{Max}}[|b_\nu\times Y|] $ for NH case applying allowed parameters that can fit the neutrino data. This implies that co-annihilations do not help reducing the Yukawa couplings to be perturbative limit. We could move to one of the next minimum model by changing the modular weight of $ N_R $ to -2 instead of zero to obtain one more mass parameters giving wider region of allowed parameters, where the other assignments are exactly the same as our model. However we would still have difficulty in realizing correct relic density while keeping perturbative limit for the Yukawa couplings. This is due to the fact that the DM annihilation cross section, Eq. (29), is p-wave dominant and we need relatively larger coupling constant than s-wave case. In addition, neutrino data and LFV constraints require heavy DM and new scalars that also suppress DM annihilation cross section. It is thus difficult to obtain correct relic density in our minimal setting and some extension is necessary.

If we forget satisfying the observed relic density and we work our numerical analysis under the perturbative limit, we obtain the tendency for electron $ g-2 $, muon $ g-2 $, and LFVs as shown in Fig. 5. These figures suggests that $ -\Delta a_e $ and $ -\Delta a_\mu $ are at most $ 10^{-20} $ and $ 10^{-15} $, respectively. On the other hand, LFVs; especially $ \mu\to e\gamma $ branching ratio, would be testable near future due to its maximum value is close to the experimental limit.

Figure 5.  (color online) Allowed region for electron $ g-2 $(left), muon $ g-2 $(center), and LFVs(right), where these points do not satisfy the observed relic density.

We briefly illustrate one of the simplest solutions to explain the observed relic density without breaking our predictions for the neutrino sector, making use of a new interaction. We can introduce a singlet scalar boson $ S_0 $ that leads to new interactions

where its modular weight is assigned to be zero for simplicity, assuming it is singlet under the $ A_4 $ symmetry, and we omit terms with $ S_{1,2}^{\pm} $. We then have Higgs portal to the SM by mixing between $ S^0 $ and h induced by the last term of $ {\cal{L}}_{\rm{new}} $. Note also that the addition of these interactions do not modify neutrino mass and the predictions in our analysis will not be changed.

As a result we have additional DM annihilation processes such as $ \chi \chi \to S^0 \to f_{\rm{SM}} f_{\rm{SM}} $ and $ \chi \chi \to S_0 S_0 $. In particular, s-channel cross section is useful to explain the relic density since annihilation cross section is enhanced at nearby $ m_\chi\approx m_{S_0}/2 $ where $ m_{S_0} $ is the mass of $ S_0 $. The annihilation cross section of $ \chi \chi \to S^0 \to f_{\rm{SM}} f_{\rm{SM}} $ process is approximately given by

where $ y_f $ is the SM Yukawa coupling for fermion f and $ \sin \alpha $ indicates the Higgs-$ S_0 $ mixing. The relic density of DM is roughly estimated as $ \Omega h^2 \sim 0.1 \ {\rm{pb}}/(\sigma v_{\rm{rel}}) $ and we obtain

where we considered top quark as f for simplicity. Thus we can realize $ \Omega h^2 \sim 0.12 $ with $ m_\chi = 1 $ TeV, $ y_S =1 $ and $ \sin \alpha \sim 0.1 $ even if we don't have resonant enhancement. With resonant effect, we can fit the relic density for small Higgs mixing case without conflicting constraints of direct detection searches [26].

We briefly illustrate one of the simplest solutions to explain the observed relic density without breaking our predictions for the neutrino sector, making use of a new interaction. We introduce a singlet scalar boson $ S_0 $ that leads to new interactions

where its modular weight is assigned to zero for simplicity, assuming it is a singlet under the $ A_4 $ symmetry, and we omit terms with $ S_{1,2}^{\pm} $. We then have a Higgs portal to the SM by mixing between $ S^0 $ and h induced by the last term of $ {\cal{L}}_{\rm{new}} $. Note that the addition of these interactions do not modify the neutrino mass, and the predictions in our analysis will not change.

As a result we have additional DM annihilation processes such as $ \chi \chi \to S^0 \to f_{\rm{SM}} f_{\rm{SM}} $ and $ \chi \chi \to S_0 S_0 $. In particular, the s-channel cross section is useful for explaining the relic density because the annihilation cross section is enhanced nearby at $ m_\chi\approx m_{S_0}/2 $, where $ m_{S_0} $ is the mass of $ S_0 $. The annihilation cross section of the $ \chi \chi \to S^0 \to f_{\rm{SM}} f_{\rm{SM}} $ process is approximately given by

where $ y_f $ is the SM Yukawa coupling for fermion f and $ \sin \alpha $ indicates the Higgs-$ S_0 $ mixing. The relic density of DM is estimated as $ \Omega h^2 \sim 0.1 \ {\rm{pb}}/(\sigma v_{\rm{rel}}) $, and we obtain

where we consider the top quark as f for simplicity. Thus, we can realize $ \Omega h^2 \sim 0.12 $ with $ m_\chi = 1 $ TeV, $ y_S =1 $, and $ \sin \alpha \sim 0.1 $, even if we do not have resonant enhancement. With the resonant effect, we can fit the relic density for the small Higgs-mixing case without conflicting constraints of direct detection searches [26].

We have investigated a three-loop induced neutrino mass model in a non-holomorphic modular flavor symmetry, in which we found some prediction in a framework that masses inside the loop does not depend on structure of the neutrino mass matrix. Since our model has a rank two Yukawa matrix in the neutrino sector, the lightest neutrino mass eigenvalue vanishes. Here, we have realized a model with minimum free parameters; three complexes $ \tau,\ \tilde c_\nu,\ \tilde d_\nu $ and five reals $ a_e,\ a_\mu,\ a_\tau,\ \tilde M_2,\ \kappa $, due to an appropriate charge-assignments under the modular symmetry. Then, we have performed the chi square analyses considering the neutrino oscillation data, and in particular, we have found rather narrow arrowed regions for the NH case while we could not fit the data in IH case. By adopting the best fit value for NH, we have further analyzed the lepton flavor violation, muon $ g-2 $, lepton flavor universalities, and dark matter, where we have neglected all the complicated processes such as co-annihilation interactions by controlling the related masses. Through the numerical analyses, we have found it is difficult to explain the observed relic density within the perturbative limit. But, it is easy to resolve it by introducing a singlet boson without changing predictions in neutrino sector.

We investigated a three-loop induced neutrino mass model in a non-holomorphic modular flavor symmetry. We observed that some predictions in a framework that masses inside the loop do not depend on the structure of the neutrino mass matrix. Because our model has a rank two Yukawa matrix in the neutrino sector, the lightest neutrino-mass eigenvalue vanishes. Here, we realized a model with minimum free parameters, three complexes $ \tau,\ \tilde c_\nu,\ \tilde d_\nu $ and five reals $ a_e,\ a_\mu,\ a_\tau,\ \tilde M_2,\ \kappa $, due to the appropriate charge assignments under the modular symmetry. Then, we performed chi-square analyses considering the neutrino-oscillation data. In particular, we observed rather narrow arrowed regions for the NH case, and we could not fit the data in the IH case. By adopting the best-fit value for NH, we further analyzed the lepton-flavor violation, muon $ g-2 $, lepton-flavor universalities, and DM. We have neglected all the complicated processes such as co-annihilation interactions by controlling the related masses. The numerical analyses showed that explaining the observed relic density within the perturbative limit is difficult. However, resolving this is easy by introducing a singlet boson without changing predictions in the neutrino sector.

The loop function at the third level is generally obtained only via numerical methods. However, if some conditions are imposed, one can analytically integrate it out. Here, we show the integration under the case of $ D_{N_i} \ll m_{S_{1,2}} $, to which we apply our model, where $ m_{S_1}^2=m_{S_2}^2 \pm \delta m_S^2 $ with $\epsilon_S\equiv \dfrac{\delta m_S}{m_{S_2}}\ll1$.

One can expand the integration in terms of $ \epsilon_i (\equiv D_{N_i}/ m_{S_1}) $and $ \epsilon_S $ as follows:

where $ a_0\approx 0.062 $, $ a_1\approx -2.92 $, $ b_1\approx -0.0281 $, and $\int [{\rm d}x]_3\equiv \int_0^1{\rm d}x\int_0^{1-x}{\rm d}y|_{z=1-x-y}$.

The loop function at the three level is in general obtained only via numerical way. But if some conditions are imposed, one can analytically integrate it out. Here, we show the integration under the case of $ D_{N_i} \ll m_{S_{1,2}} $ to which we apply our model where $ m_{S_1}^2=m_{S_2}^2 \pm \delta m_S^2 $ with $ \epsilon_S\equiv \frac{\delta m_S}{m_{S_2}}\ll1 $.

One can expand the integration in terms of $ \epsilon_i (\equiv D_{N_i}/m_{S_1}) $ and $ \epsilon_S $ as

Here $ a_0\approx 0.062 $, $ a_1\approx -2.92 $, and $ b_1\approx -0.0281 $ and $ \int [dx]_3\equiv \int_0^1dx\int_0^{1-x}dy|_{z=1-x-y} $.