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

Thanks to a successful construction of non-holomorphic modular symmetry framework via Qu and Ding [1], we can safely deal with a beyond the standard model (BSM) without super-symmetric theories in 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 some advantage of applying non-supersymmetric framework to reduce number of new fields where extra fields are sometimes required to cancel a gauge anomaly in supersymmetric case.

Radiatively induced neutrino mass models are one of the representative scenarios that do not demand the super-symmetric framework and one can naturally connect new particles to the SM model particles. Sometimes, the model can possess the dark matter candidate(DM) [13] that often requests an additional symmetry to stabilize it. Thus, constructing radiative neutrino mass models (with DM) using the non-holomorphic modular symmetry would be natural manner to make a model more attractive realizing more predictability.

In our paper, 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 since the scale of new particles would be smaller compared to lower loop (or tree) level models due to loop suppression. We then expect rich phenomenology such as collider and lepton flavor physics. Non-holomorphic modular symmetry framework is suitable to construct such a three loop model in minimal way; if we consider holomorphic one we need to add more fields to cancel gauge anomaly. Then we try to find the assignment under non-holomorphic modular $ A_4 $ symmetry as possible as minimal number of free parameters to fit the observables in lepton sector. Through the chi square numerical analysis, we successfully search for the minimum model to predict the lepton masses and mixing angles in addition to reproduce the current neutrino observables in Nufit 6.0 [15]. Then, we perform further numerical analyses in order to satisfy lepton flavor violations (LFVs), muon anomalous magnetic moment, (muon $ g-2 $), lepton universality, and the dark matter (DM). As a result, we find that relic density is too large within the limit of perturbation requiring a new interaction where it can be done by adding one singlet scalar boson without changing predictions in neutrino sector.

This paper is organized as follows. In Sec. II, we explain our minimum three-loop neutrino mass model constructing the renormalizable Lagrangian in the lepton sector, Higgs sector, the charged-lepton sector, heavier Majorana fermion sector, and the 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 show some predictions for normal and inverted hierarchies in the neutrino sector. 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 the relic density of DM. We have conclusions and discussion in Sec. IV. In Appendix A, we show the three-loop function in the neutrino sector.

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 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 $.

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 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].

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 $.

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.

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.

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 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 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 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.

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} $.