$\displaystyle\text{the multiplication $"\diamond"$:}\quad e_{2}\diamond
e_{3}=\lambda_{1}e_{4},\quad e_{3}\diamond e_{2}=\lambda_{2}e_{4},\quad
e_{3}\diamond e_{3}=\lambda_{3}e_{1},\quad\lambda_{i}\in\mathbb{K},$
$\displaystyle\text{the even linear map $\alpha:A\rightarrow A$ given
by:}\quad\begin{array}[t]{ll}\alpha(e_{1})=2e_{1},&\quad\alpha(e_{2})=-e_{2},\\\
\alpha(e_{3})=-e_{3},&\quad\alpha(e_{4})=-2e_{4}.\end{array}$

Therefore, $(A,\cdot,[\cdot,\cdot],\varepsilon,\alpha)$ is a Hom-Gelfand-
Dorfman color Hom-algebra with

$[e_{2},e_{3}]=-[e_{3},e_{2}]=(\lambda_{1}-\lambda_{2})e_{4},\quad[e_{3},e_{3}]=2\lambda_{3}e_{3}.$

######  Lemma 4.18 ([5]).

Let $(A,\cdot,\varepsilon,\alpha)$ be a $\varepsilon$-commutative Hom-
associative color Hom-algebra and $D$ be an even derivation. Define a bilinear
operation $\diamond$ on $A$ by

$x\diamond y=x\cdot D(y),\quad\forall x,y\in\mathcal{H}(A).$ (4.33)

Then $(A,\cdot,\diamond,\varepsilon,\alpha)$ is a Hom-Novikov-Poisson color
Hom-algebras.

Combining Theorem 4.16 and Lemma 4.18 leads to the following corollary.

######  Corollary 4.19.

Let $(A,\cdot,\varepsilon,\alpha)$ be a $\varepsilon$-commutative Hom-
associative color Hom-algebra and $D$ be an even derivation. Then
$(A,\cdot,[\cdot,\cdot],\varepsilon,\alpha)$ is a Hom-Gelfand-Dorfman color
Hom-algebra, where for all $x,y\in\mathcal{H}(A)$,

$[x,y]=x\cdot D(y)-\varepsilon(x,y)y\cdot D(x).$ (4.34)

Next theorem provides a procedure for construction of the Hom-Gelfand-Dorfman
color Hom-algebras from Gelfand-Dorfman color Hom-algebras and morphisms.

######  Theorem 4.20.

Let $\mathcal{A}=(A,\cdot,[\cdot,\cdot],\varepsilon)$ be a Gelfand-Dorfman
color Hom-algebra and $\alpha:\mathcal{A}\rightarrow\mathcal{A}$ be a Gelfand-
Dorfman color Hom-algebras morphism. Define $\cdot_{\alpha},\
[\cdot,\cdot]_{\alpha}:A\times A\rightarrow A$ for all $x,y\in\mathcal{H}(A)$,
by $x\cdot_{\alpha}y=\alpha(x\cdot y)$ and $[x,y]_{\alpha}=\alpha([x,y])$.
Then,
$\mathcal{A}_{\alpha}=(A_{\alpha}=A,\cdot_{\alpha},[x,y]_{\alpha},\varepsilon,\alpha)$
is a Hom-Gelfand-Dorfman color Hom-algebra called the $\alpha$-twist or Yau
twist of $(A,\cdot,[\cdot,\cdot],\varepsilon)$. Moreover, assume that
$\mathcal{A}^{\prime}=(A^{\prime},\cdot^{\prime},[\cdot,\cdot]^{\prime},\varepsilon)$
is another Hom-Gelfand-Dorfman color Hom-algebra and
$\alpha^{\prime}:\mathcal{A}^{\prime}\rightarrow\mathcal{A}^{\prime}$ is a
Hom-Gelfand-Dorfman color Hom-algebras morphism. Let
$f:\mathcal{A}\rightarrow\mathcal{A}^{\prime}$ be a Hom-Gelfand-Dorfman color
Hom-algebras morphism satisfying $f\circ\alpha=\alpha^{\prime}\circ f$. Then,
$f:\mathcal{A}_{\alpha}\rightarrow\mathcal{A}^{\prime}_{\alpha}$ is a Hom-
Gelfand-Dorfman color Hom-algebras morphism.

###### Proof.

Being a Gelfand-Dorfman color Hom-algebras morphism, $\alpha:A\rightarrow A$
is an even linear map which is multiplicative with respect to $\cdot$ and
$[\cdot,\cdot]$, that is,

$\forall\ x,y\in\mathcal{H}(A):\quad\alpha(x\cdot
y)=\alpha(x)\cdot\alpha(y),\quad\alpha([x,y])=[\alpha(x),\alpha(y)].$

The equality (4.2) in $\mathcal{A}_{\alpha}$ is proved as follows:

$\displaystyle\alpha(y)\cdot_{\alpha}[x,z]_{\alpha}=\alpha(y)\cdot_{\alpha}\alpha([x,z])=\alpha(\alpha(y)\cdot\alpha([x,z]))$
($\alpha$ morphism)
$\displaystyle=\alpha^{2}(y)\cdot\alpha^{2}([x,z])=\alpha^{2}(y)\cdot[\alpha^{2}(x),\alpha^{2}(z)]$
($\mathcal{A}$ is a Hom-G. D. color alg)
$\displaystyle=\varepsilon(y,x)[\alpha^{2}(x),\alpha^{2}(y)\cdot\alpha^{2}(z)]-\varepsilon(x+y,z)[\alpha^{2}(z),\alpha^{2}(y)\cdot\alpha^{2}(x)]$
$\displaystyle\quad+[\alpha^{2}(y),\alpha^{2}(x)]\cdot\alpha^{2}(z)-\varepsilon(x,z)[\alpha^{2}(y),\alpha^{2}(z)]\cdot\alpha^{2}(x)$
($\alpha$ morphism)
$\displaystyle=\varepsilon(y,x)[\alpha^{2}(x),\alpha(\alpha(y)\cdot\alpha(z))]-\varepsilon(x+y,z)[\alpha^{2}(z),\alpha(\alpha(y),\alpha(x))]$
$\displaystyle\quad+\alpha([\alpha(y),\alpha(x)])\cdot\alpha^{2}(z)-\varepsilon(x,z)\alpha([\alpha(y),\alpha(z)]\cdot\alpha^{2}(x)$
$\displaystyle=\varepsilon(y,x)[\alpha^{2}(x),\alpha(y\cdot_{\alpha}z)]-\varepsilon(x+y,z)[\alpha^{2}(z),\alpha(y\cdot_{\alpha}x)]$
$\displaystyle\quad+\alpha([y,x]_{\alpha}\cdot\alpha^{2}(z)-\varepsilon(x,z)\alpha([y,z]_{\alpha})\cdot\alpha^{2}(x)$
$\displaystyle=\varepsilon(y,x)[\alpha(x),y\cdot_{\alpha}z]_{\alpha}-\varepsilon(x+y,z)[\alpha(z),y\cdot_{\alpha}x]_{\alpha}$
$\displaystyle\quad+[y,x]_{\alpha}\cdot_{\alpha}\alpha(z)-\varepsilon(x,z)[y,z]_{\alpha}\cdot_{\alpha}\alpha(x).$

The second assertion follows from

$\displaystyle f(x\cdot_{\alpha}y)$ $\displaystyle=f(\alpha(x\cdot
y))=\alpha^{\prime}(f(x\cdot
y))=\alpha^{\prime}(f(x)\cdot^{\prime}f(y))=f(x)\cdot^{\prime}_{\alpha^{\prime}}f(y),$
$\displaystyle f([x,y]_{\alpha})$
$\displaystyle=f(\alpha([x,y]))=\alpha^{\prime}(f([x,y]))=\alpha^{\prime}([f(x),f(y)]^{\prime})=[f(x),f(y)]^{\prime}_{\alpha}.\qed$

######  Corollary 4.21.

If $\mathcal{A}=(A,\cdot,[\cdot,\cdot],\varepsilon,\alpha)$ is a
multiplicative Hom-Gelfand-Dorfman color algebra, then for any
$n\in\mathbb{N}^{\ast}$,

  1. (i)

The $n{\rm th}$ derived Hom-Gelfand-Dorfman color Hom-algebra of type $1$ of
$\mathcal{A}$ is defined by

$\mathcal{A}_{1}^{n}=(A,\cdot^{(n)}=\alpha^{n}\circ\cdot,\ast^{(n)}=\alpha^{n}\circ[\cdot,\cdot],\varepsilon,\alpha^{n+1}).$

  2. (ii)

The $n{\rm th}$ derived Hom-Gelfand-Dorfman color Hom-algebra of type $2$ of
$A$ is defined by

$\mathcal{A}_{2}^{n}=(A,\cdot^{(2^{n}-1)}=\alpha^{2^{n}-1}\circ\cdot,[\cdot,\cdot]^{(2^{n}-1)}=\alpha^{2^{n}-1}\circ[\cdot,\cdot],\varepsilon,\alpha^{2^{n}}).$

###### Proof.

Apply Theorem 4.20 with $\alpha^{\prime}=\alpha^{n}$ and
$\alpha^{\prime}=\alpha^{2^{n}-1}$ respectively. ∎

######  Example 4.22.

Let $\Gamma=\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ and $A$ be a $4$-dimensional
$\Gamma$-graded linear space with
$A_{(0,0)}=<e_{1}>,A_{(0,1)}=<e_{2}>,A_{(1,0)}=<e_{3}>,A_{(1,1)}=<e_{4}>.$
Then there is a multiplicative admissible Hom-Gelfand-Dorfman color Hom-
algebra $(A,\cdot,[\cdot,\cdot],\varepsilon,\alpha)$ with the bicharacter
$\varepsilon\big{(}(i_{1},i_{2}),(j_{1},j_{2})\big{)}=(-1)^{i_{1}j_{1}+i_{2}j_{2}},$
and the multiplications tables for a basis $\\{e_{1},e_{2},e_{3},e_{4}\\}$:

$\cdot$ | $e_{1}$ | $e_{2}$ | $e_{3}$ | $e_{4}$
---|---|---|---|---
$e_{1}$ | $0$ | $0$ | $0$ | $0$
$e_{2}$ | $0$ | $0$ | $2e_{4}$ | $0$
$e_{3}$ | $0$ | $2e_{4}$ | $e_{1}$ | $0$
$e_{4}$ | $0$ | $0$ | $0$ | $0$
$[\cdot,\cdot]$ | $e_{1}$ | $e_{2}$ | $e_{3}$ | $e_{4}$
---|---|---|---|---
$e_{1}$ | $0$ | $0$ | $0$ | $0$
$e_{2}$ | $0$ | $0$ | $-2e_{4}$ | $0$
$e_{3}$ | $0$ | $2e_{4}$ | $0$ | $0$
$e_{4}$ | $0$ | $0$ | $0$ | $0$
$\begin{array}[]{ll}\alpha(e_{1})=e_{1},&\alpha(e_{2})=-2e_{2},\\\
\alpha(e_{3})=-e_{3},&\alpha(e_{4})=2e_{4},\end{array}$

Then there are Hom-Gelfand-Dorfman color Hom-algebras $\mathcal{A}_{1}^{n}$
and $\mathcal{A}_{2}^{n}$ with multiplications tables respectively:

$\cdot^{(n)}$ | $e_{1}$ | $e_{2}$ | $e_{3}$ | $e_{4}$
---|---|---|---|---
$e_{1}$ | $0$ | $0$ | $0$ | $0$
$e_{2}$ | $0$ | $0$ | $2^{n}e_{4}$ | $0$
$e_{3}$ | $0$ | $2^{n}e_{4}$ | $e_{1}$ | $0$
$e_{4}$ | $0$ | $0$ | $0$ | $0$
$[\cdot,\cdot]^{(n)}$ | $e_{1}$ | $e_{2}$ | $e_{3}$ | $e_{4}$
---|---|---|---|---
$e_{1}$ | $0$ | $0$ | $0$ | $0$
$e_{2}$ | $0$ | $0$ | $(-2)^{n}e_{4}$ | $0$
$e_{3}$ | $0$ | $2^{n}e_{4}$ | $0$ | $0$
$e_{4}$ | $0$ | $0$ | $0$ | $0$
$\begin{array}[]{ll}\alpha^{n+1}(e_{1})=e_{1},&\alpha^{n+1}(e_{2})=(-2)^{n+1}e_{2},\\\ \alpha^{n+1}(e_{3})=(-1)^{n+1}e_{3},&\alpha^{n+1}(e_{4})=2^{n+1}e_{4},\end{array}$ $\cdot^{(2^{n}-1)}$ | $e_{1}$ | $e_{2}$ | $e_{3}$ | $e_{4}$
---|---|---|---|---
$e_{1}$ | $0$ | $0$ | $0$ | $0$
$e_{2}$ | $0$ | $0$ | $2^{2^{n}-1}e_{4}$ | $0$
$e_{3}$ | $0$ | $2^{2^{n}-1}e_{4}$ | $e_{1}$ | $0$
$e_{4}$ | $0$ | $0$ | $0$ | $0$
$[\cdot,\cdot]^{(2^{n}-1)}$ | $e_{1}$ | $e_{2}$ | $e_{3}$ | $e_{4}$
---|---|---|---|---
$e_{1}$ | $0$ | $0$ | $0$ | $0$
$e_{2}$ | $0$ | $0$ | $-2^{2^{n}-1}e_{4}$ | $0$
$e_{3}$ | $0$ | $2^{2^{n}-1}e_{4}$ | $0$ | $0$
$e_{4}$ | $0$ | $0$ | $0$ | $0$
$\alpha^{2^{n}}(e_{1})=e_{1},\quad\alpha^{2^{n}}(e_{2})=2^{2^{n}}e_{2},\quad\alpha^{2^{n}}(e_{3})=e_{3},\quad\alpha^{2^{n}}(e_{4})=2^{2^{n}}e_{4}.$

## References

  * [1] Abdaoui, K., Ammar, F., Makhlouf, A.: Constructions and cohomology of Hom-Lie color Hom-algebras, Comm. Algebra, 43, 4581-4612 (2015)
  * [2] Abdaoui, E., Mabrouk, S., Makhlouf, A.: Rota-Baxter Operators on Pre-Lie Superalgebras, Bulletin of the Malaysian Math. Sci. Soc. (2) 42 (4), 1567-1606 (2019)
  * [3] Ammar, F., Ejbehi, Z., Makhlouf, A.: Cohomology and deformations of Hom-algebras, J. Lie Theory, 21(4), 813-836 (2011)
  * [4] Arnold, V.I.: Mathematical methods of classical mechanics, Grad. Texts in Math, Springer, Berlin, 60, (1978)
  * [5] Attan, S.: Some constructions of color Hom-Novikov-Poisson algebras, Mathematical Sciences and Applications E-Notes, 1, 78-86 (2019)
  * [6] Attan, S., Laraiedh, I.: Construtions and bimodules of BiHom-alternative and BiHom-Jordan algebras, arXiv:2008.07020 [math.RA] (2020)
  * [7] Bahturin, Y. A., Mikhalev, A. A., Petrogradsky, V. M., Zaicev, M. V.: Infinite-dimensional Lie superalgebras, De Gruyter Expositions in Mathematics 7, Walter de Gruyter, Berlin (1992)
  * [8] Bakayoko, I.: Hom-Novikov color Hom-algebras, arXiv:1609.07813 [math.RA], 16pp (2016)
  * [9] Bakayoko, I.: Laplacian of Hom-Lie quasi-bialgebras, International Journal of Algebra, 8(15), 713-727 (2014)
  * [10] Bakayoko, I.: $L$-modules, $L$-comodules and Hom-Lie quasi-bialgebras, African Diaspora Journal of Mathematics, 17, 49-64 (2014)
  * [11] Bakayoko, I. Attan, S.: Some structures of Hom-Poisson color Hom-algebras, arXiv:1912.01358v1[math.RA] (2019).
  * [12] Bakayoko, I., Banagoura, M.: Bimodules and Rota-Baxter Relations. J. Appl. Mech. Eng. 4(5) (2015)
  * [13] Bakayoko, I., Silvestrov, S.: Multiplicative $n$-Hom-Lie color Hom-algebras, In: Silvestrov, S., Malyarenko, A., Ranc̆ić, M. (Eds.), Algebraic Structures and Applications, Springer Proceedings in Mathematics and Statistics 317, Ch. 7, 159-187, Springer (2020). (arXiv:1912.10216[math.QA] (2019))
  * [14] Bakayoko, I., Silvestrov, S.: Hom-left-symmetric color dialgebras, Hom-tridendriform color Hom-algebras and Yau’s twisting generalizations, Afr. Mat. (2021).
(arXiv:1912.01441[math.RA] (2019))

  * [15] Balinskii, A.A., Novikov, S.P.: Poisson brackets of hydrodynamic type, Frobenius algebras and Lie algebras, Soviet Math. Dokl. 32, 228-231 (1985) Translated from Russian: Dokl. Akad. Nauk SSSR, 283:5 (1985), 1036-1039)
  * [16] Benayadi, S., Makhlouf, A.: Hom-Lie algebras with symmetric invariant nondegenerate bilinear forms, J. Geom. Phys. 76, 38-60 (2014)
  * [17] Ben Abdeljelil, A., Elhamdadi, M., Kaygorodov, I., Makhlouf, A.: Generalized Derivations of $n$-BiHom-Lie algebras, In: Silvestrov, S., Malyarenko, A., Ranc̆ić, M. (Eds.), Algebraic Structures and Applications, Springer Proceedings in Mathematics and Statistics 317, Ch. 4, 81-97, Springer (2020). (arXiv:1901.09750[math.RA] (2019))
  * [18] Ben Hassine, A., Chtioui, T., Mabrouk, S., Ncib, O.: Cohomology and linear deformation of BiHom-left-symmetric algebras, arXiv:1907.06979 [math.RA], 19 pp (2019)
  * [19] Calderon, A., Delgado, J. S.: On the structure of split Lie color Hom-algebras, Linear Algebra Appl. 436, 307-315 (2012)
  * [20] Caenepeel, S., Goyvaerts, I.: Monoidal Hom-Hopf Algebras, Comm. Algebra, 39(6), 2216-2240 (2011)
  * [21] Charim, V., Pressley, A.N.: A guide to quantum groups, Cambridge Univ. Press, Cambridge, (1994).
  * [22] Chen, X.-W., Silvestrov, S. D., van Oystaeyen, F.: Representations and Cocycle Twists of Color Lie Algebras, Algebr. Represent. Theory, 9(6), 633-650 (2006)
  * [23] Chen, C. W., Petit, T., Van Oystaeyen, F.: Note on cohomology of color Hopf and Lie algebras, J. Algebra, 299, 419-442 (2006)
  * [24] Cao, Y., Chen, L.: On split regular Hom-Lie color Hom-algebras, Comm. Algebra 40, 575-592 (2012)
  * [25] Chtioui, T., Mabrouk, S., Makhlouf, A.: BiHom-pre-alternative algebras and BiHom-alternative quadri-algebras, Bull. Math. Soc. Sci. Math. Roumanie. 63 (111)(1), 3-21 (2020)
  * [26] Dassoundo, M. L., Silvestrov, S.: Nearly associative and nearly Hom-associative algebras and bialgebras, arXiv:2101.12377 [math.RA], 24pp (2021)
  * [27] Drinfel’d, V. G.: Quantum groups, in: Proc. ICM (Berkeley, 1986), AMS, Providence, RI, 798-820 (1987).
  * [28] Dubrovin, B. A., Novikov, S. P.: Hamiltonian formalism of one-dimensional systems of hydrodynamic type and the Bogolyubov-Whitham averaging method, Soviet Math. Dokl.27, 665-669 (1983) (Translated from Russian: Dokl. Akad. Nauk SSSR 270, 781-785 (1983))
  * [29] Dubrovin, B. A., Novikov, S. P.: On Poisson brackets of hydrodynamic type, Soviet Math. Dokl. 30, 651-654 (1984) (Transalated from Russian: Dokl. Akad. Nauk SSSR 279:2, 294–297 (1984))
  * [30] Ebrahimi-Fard, K., Guo, L.: Rota-Baxter algebras and dendriform algebras, J. Pure Appl. Algebra, 212(2), 320-339 (2008)
  * [31] Frenkel, E., Ben-Zvi, D.: Vertex algebras and algebraic curves, Math Surveys and Monographs 88, 2nd ed., AMS, Providence, RI, (2004)
  * [32] Gel’fand, I. M., Dickey, L. A.: Asymptotic behavior of the resolvent of Sturm-Liouville equations and the Lie algebras of the Korteweg-de Vries equations, Russian Math. Sur. 30(5(185)), 77-113 (1975) (Translated from Russian: Uspekhi Mat. Nauk, 30(5), 67-100 (1975))
  * [33] Gel’fand, I. M., Dickey, L. A.: A Lie algebra structure in a formal variational calculations, Funct. Anal. Its Appl. 10(1), 16-22 (1976) (Translated from Russian: Funktsional’nyi Analiz i Ego Prilozheniya, 10(1), 16-25 (1976))
  * [34] Gel’fand, I. M., Dorfman, I. Y. Hamiltonian operators and algebraic structures related to them, Funct. Anal. Its Appl. 13(4), 248-262 (1979) (Translated from Russian: Funktsional’nyi Analiz i Ego Prilozheniya, 13(4), 13-30 (1979))
  * [35] Gerstenhaber, M.: On the deformation of rings and algebras, Ann. Math. 79, 59-103 (1964)
  * [36] Graziani, G., Makhlouf, A., Menini, C., Panaite, F.: BiHom-Associative Algebras, BiHom-Lie Algebras and BiHom-Bialgebras, Symmetry, Integrability Geom.: Methods Appl. (SIGMA), 11(086), 34 pp (2015)
  * [37] Hartwig, J. T., Larsson, D., Silvestrov, S. D.: Deformations of Lie algebras using $\sigma$-derivations, J. Algebra, 295, 314-361 (2006) (Preprint in Mathematical Sciences 2003:32, LUTFMA-5036-2003, Centre for Mathematical Sciences, Department of Mathematics, Lund Institute of Technology, 52 pp. (2003))
  * [38] Hassanzadeh, M., Shapiro, I., Sütlü, S.: Cyclic homology for Hom-associative algebras, J. Geom. Phys. 98, 40-56 (2015)
  * [39] Hounkonnou, M. N., Dassoundo, M. L.: Center-symmetric Algebras and Bialgebras: Relevant Properties and Consequences. In: Kielanowski P., Ali S., Bieliavsky P., Odzijewicz A., Schlichenmaier M., Voronov T. (eds) Geometric Methods in Physics. Trends in Mathematics. 2016, pp. 281-293. Birkhäuser, Cham (2016)
  * [40] Hounkonnou, M. N., Houndedji, G. D., Silvestrov, S.: Double constructions of biHom-Frobenius algebras, arXiv:2008.06645 [math.QA], 45pp (2020)
  * [41] Hounkonnou, M. N., Dassoundo, M. L.: Hom-center-symmetric algebras and bialgebras. arXiv:1801.06539 [math.RA], 19 pp (2018)
  * [42] Kitouni, A., Makhlouf, A., Silvestrov, S.: On $n$-ary generalization of BiHom-Lie algebras and BiHom-associative algebras, In: Silvestrov, S., Malyarenko, A., Ranc̆ić, M. (Eds.), Algebraic Structures and Applications, Springer Proceedings in Mathematics and Statistics 317, Springer, Ch. 5, 99-126 (2020)
  * [43] Kontsevich, M.: Deformation quantization of Poisson manifolds, Lett. Math. Phys. 66, 157-216 (2003)
  * [44] Mabrouk, S., Ncib, O., Silvestrov, S.: Generalized Derivations and Rota-Baxter Operators of $n$-ary Hom-Nambu Superalgebras, Adv. Appl. Clifford Algebras, 31, 32, (2021) (arXiv:2003.01080[math.QA] (2020))
  * [45] Makhlouf, A., Silvestrov, S. D.: Hom-algebra structures. J. Gen. Lie Theory Appl. 2(2), 51–64 (2008) (Preprints in Mathematical Sciences 2006:10, LUTFMA-5074-2006, Centre for Mathematical Sciences, Department of Mathematics, Lund Institute of Technology, Lund University (2006))
  * [46] Makhlouf, A., Silvestrov, S.: Notes on $1$-parameter formal deformations of Hom-associative and Hom-Lie algebras, Forum Math. 22(4), 715-739 (2010) (Preprints in Mathematical Sciences, Lund University, Centre for Mathematical Sciences, Centrum Scientiarum Mathematicarum, (2007:31) LUTFMA-5095-2007. arXiv:0712.3130v1 [math.RA] (2007))
  * [47] Makhlouf, A., Silvestrov, S. D.: Hom-algebras and Hom-coalgebras, J. Algebra Appl. 9(04), 553-589 (2010) (Preprints in Mathematical Sciences, Lund University, Centre for Mathematical Sciences, Centrum Scientiarum Mathematicarum, (2008:19) LUTFMA-5103-2008. arXiv:0811.0400 [math.RA] (2008))
  * [48] Makhlouf, A., Yau, D.: Rota-Baxter Hom-Lie admissible algebras, Comm. Alg., 23(3), 1231-1257 (2014)
  * [49] Ma, T., Zheng, H.: Some results on Rota-Baxter monoidal Hom-algebras, Results Math. 72 (1-2), 145-170 (2017)
  * [50] Mikhalev, A. A., Zolotykh, A. A.: Combinatorial Aspects of Lie Superalgebras, CRC Press 1995
  * [51] Laraiedh, I.: Bimodules and matched pairs of noncommutative BiHom-(pre)-Poisson algebras, arXiv:2102.11364 [math.RA], 30 pp (2021)
  * [52] Larsson, D., Sigurdsson, G., Silvestrov, S. D.: Quasi-Lie deformations on the algebra $\mathbb{F}[t]/(t^{N})$, J. Gen. Lie Theory Appl. 2(3), 201-205 (2008)
  * [53] Larsson, D., Silvestrov, S. D.: Quasi-deformations of $sl_{2}(\mathbb{F})$ using twisted derivations, Comm. Algebra, 35, 4303-4318 (2007) (Preprint in Mathematical Sciences 2004:26, LUTFMA-5047-2004, Centre for Mathematical Sciences, Lund Institute of Technology, Lund University (2004). arXiv:math/0506172 [math.RA] (2005))
  * [54] Larsson, D., Silvestrov, S. D.: Quasi-Hom-Lie algebras, central extensions and $2$-cocycle-like identities, J. Algebra, 288, 321-344 (2005) (Preprints in Mathematical Sciences 2004:3, LUTFMA-5038-2004, Centre for Mathematical Sciences, Department of Mathematics, Lund Institute of Technology, Lund University (2004))
  * [55] Larsson, D., Silvestrov, S. D.: Quasi-Lie algebras, In: Fuchs, J., Mickelsson, J., Rozenblioum, G., Stolin, A., Westerberg, A. (Eds.), Noncommutative Geometry and Representation Theory in Mathematical Physics, Contemp. Math. 391, Amer. Math. Soc., Providence, RI, 241-248 (2005) (Preprints in Mathematical Sciences 2004:30, LUTFMA-5049-2004, Centre for Mathematical Sciences, Department of Mathematics, Lund Institute of Technology, Lund University (2004))
  * [56] Larsson, D., Silvestrov, S. D.: Graded quasi-Lie agebras, Czechoslovak J. Phys. 55, 1473-1478 (2005)
  * [57] Larsson, D., Silvestrov, S. D.: On generalized $N$-complexes comming from twisted derivations, In: Silvestrov, S., Paal, E., Abramov, V., Stolin, A. (Eds.), Generalized Lie Theory in Mathematics, Physics and Beyond, Springer-Verlag, Ch. 7, 81-88 (2009)
  * [58] Liu, L., Makhlouf, A., Menini, C., Panaite, F.: Rota-Baxter operators on BiHom-associative algebras and related structures, Colloq. Math. 161(2), 263-294 (2020) arXiv:1703.07275 [math.RA], 27pp (2017)
  * [59] Piontkovski, D., Silvestrov, S.: Cohomology of 3-dimensional color Lie algebras, J. Algebra, 316(2), 499-513 (2007)
  * [60] Richard, L., Silvestrov, S. D.: Quasi-Lie structure of $\sigma$-derivations of $\mathbb{C}[t^{\pm 1}]$, J. Algebra, 319(3), 1285-1304 (2008) (arXiv:math/0608196[math.QA] (2006). Preprints in mathematical sciences (2006:12), LUTFMA-5076-2006, Centre for Mathematical Sciences, Lund University (2006))
  * [61] Richard, L., Silvestrov, S. D.: A note on quasi-Lie and Hom-Lie structures of $\sigma$-derivations of ${\mathbb{C}}[z_{1}^{\pm 1},\ldots,z_{n}^{\pm 1}]$, In: Silvestrov, S., Paal, E., Abramov, V., Stolin, A. (Eds.), Generalized Lie Theory in Mathematics, Physics and Beyond, Springer-Verlag, Ch. 22, 257-262, (2009)
  * [62] Saadaou, N., Silvestrov, S.: On $(\lambda,\mu,\gamma)$-derivations of BiHom-Lie algebras, arXiv:2010.09148 [math.RA], (2020)
  * [63] Schaller, P., Strobl, T.: Poisson structure induced (topological) field theories, Mod. Phys. Lett. 9, 3129-3136 (1994)
  * [64] Scheunert, M.: Generalized Lie algebras, J. Math. Phys. 20(4), 712-720 (1979)
  * [65] Scheunert, M.: Introduction to the cohomology of Lie superalgebras and some applications, Res. Exp. Math. 25, 77-107 (2002)
  * [66] Scheunert, M., Zhang, R. B.: Cohomology of Lie superalgebras and their generalizations, J. Math. Phys. 39, 5024-5061 (1998)
  * [67] Sheng, Y.: Representations of Hom-Lie algebras, Algebr. Reprensent. Theory 15, 1081-1098 (2012)
  * [68] Sheng, Y., Bai, C.: A new approach to Hom-Lie bialgebras, J. Algebra, 399, 232-250 (2014)
  * [69] Sigurdsson, G., Silvestrov, S.: Graded quasi-Lie algebras of Witt type, Czechoslovak J. Phys. 56, 1287-1291 (2006)
  * [70] Sigurdsson, G., Silvestrov, S.: Lie color and Hom-Lie algebras of Witt type and their central extensions, In: Silvestrov, S., Paal, E., Abramov, V., Stolin, A. (Eds.), Generalized Lie Theory in Mathematics, Physics and Beyond, Springer-Verlag, Berlin, Heidelberg, Ch. 21, 247-255 (2009)
  * [71] Silvestrov, S. D.: On the classification of $3$-dimensional coloured Lie algebras, In: Quantum groups and quantum spaces (Warsaw, 1995), Banach Center Publ. 40, Polish Acad. Sci., Warsaw, 159-170 (1997)
  * [72] Silvestrov, S.: Paradigm of quasi-Lie and quasi-Hom-Lie algebras and quasi-deformations, In ”New techniques in Hopf algebras and graded ring theory”, K. Vlaam. Acad. Belgie Wet. Kunsten (KVAB), Brussels, 165-177 (2007)
  * [73] Silvestrov, S., Zargeh, C.: HNN-extension of involutive multiplicative Hom-Lie algebras, arXiv:2101.01319 [math.RA], 14pp (2021)
  * [74] Sun, Q.: On Hom-Prealternative Bialgebras, Algebr. Represent. Theor. 19, 657-677 (2016)
  * [75] Sun, Q., Li, H.: On parakähler Hom-Lie algebras and Hom-left-symmetric bialgebras, Comm. Algebra, 45(1), 105-120 (2017)
  * [76] Vaisman, I.: Lectures on the geometry of Poisson manifolds, Birkhäuser, Basel, 1994
  * [77] Xu, X.: On simple Novikov algebras and their irreducible modules, J. Algebra, 185, 905-934 (1996)
  * [78] Xu, X.: Novikov-Poisson algebras, J. Algebra, 190, 253-279 (1997)
  * [79] Xu, X.: Quadratic Conformal Superalgebras, J. Algebra 231, 1-38 (2000)
  * [80] Yuan, L.: Hom-Lie color Hom-algebras, Comm. Alg., 40(2), 575-592 (2012)
  * [81] Yau, D.: Module Hom-algebras, arXiv:0812.4695[math.RA], 10pp (2008)
  * [82] Yau D.: The Hom-Yang-Baxter equation, Hom-Lie algebras and quasi-triangular bialgebras, J. Phys. A.: Math. Theor. 42(16), 165-202 (2009)
  * [83] Yau, D.: Hom-algebras and homology, J. Lie Theory 19(2), 409-421 (2009)
  * [84] Yau, D.: Hom-Novikov algebras, arXiv:0909.0726 (2009)
  * [85] Yau, D.: A twisted generalization of Novikov-Poisson algebras, arXiv:1010.3410 [math.RA] (2010)
  * [86] Yau, D.: Hom-bialgebras and comodule Hom-algebras, Int. Electron. J. Algebra, 8, 45-64 (2010) (arXiv:0810.4866[math.RA] (2008))
  * [87] Yau, D.: Hom-Malcev, Hom-alternative, and Hom-Jordan algebras, Int. Electron. J. Algebra, 11, 177-217 (2012)
  * [88] Yuan, L.: Hom Gel’fand-Dorfman bialgebras and Hom-Lie conformal algebras, J. Math. Phys. 55(4), 043507 (2014)


# Lepton flavor violation and scotogenic Majorana neutrino mass in a
Stueckelberg $U(1)_{X}$ model

Chuan-Hung Chen<EMAIL_ADDRESS>Department of Physics, National
Cheng-Kung University, Tainan 70101, Taiwan  Physics Division, National Center
for Theoretical Sciences, Taipei 10617, Taiwan  Cheng-Wei Chiang
<EMAIL_ADDRESS>Department of Physics and Center for Theoretical
Physics, National Taiwan University, Taipei 10617, Taiwan  Physics Division,
National Center for Theoretical Sciences, Taipei 10617, Taiwan  Takaaki Nomura
<EMAIL_ADDRESS>College of Physics, Sichuan University, Chengdu 610065,
China  Chun-Wei Su<EMAIL_ADDRESS>Department of Physics and Center for
Theoretical Physics, National Taiwan University, Taipei 10617, Taiwan

###### Abstract

We construct a scotogenic Majorana neutrino mass model in a gauged $U(1)_{X}$
extension of the standard model, where the mass of the gauge boson and the
unbroken gauge symmetry, which leads to a stable dark matter (DM), can be
achieved through the Stueckelberg mechanism. It is found that the simplest
version of the extended model consists of the two inert-Higgs doublets and one
vector-like singlet fermion. In addition to the Majorana neutrino mass, we
study the lepton flavor violation (LFV) processes, such as
$\ell_{i}\to\ell_{j}\gamma$, $\ell_{i}\to 3\ell_{j}$, $\mu-e$ conversion rate
in nucleus, and muonium-antimuonium oscillation. We show that the
sensitivities of $\mu\to 3e$ and $\mu-e$ conversion rate designed in Mu3e and
COMET/Mu2e experiments make both decays the most severe constraints on the
$\mu\to e$ LFV processes. It is found that $\tau\to\mu\gamma$ and $\tau\to
3\mu$ can reach the designed significance level of Belle II. In addition to
explaining the DM relic density, we also show that the DM-nucleon scattering
cross section can satisfy the currently experimental limit of DM direct
detection.

##  I Introduction

A scotogenic mechanism, which generate the neutrino mass radiatively, was
proposed in Ref. Ma:2006km  and is referred to as the Ma model hereinafter.
The model not only resolves the origin of neutrino mass but also provides a
way to explain the dark matter (DM) relic density, where the current value
observed by Planck collaboration is $\Omega_{\rm DM}h^{2}=0.120\pm 0.001$
Planck:2018vyg .

Since the active neutrinos and the charged leptons form $SU(2)_{L}$ doublets
in the standard model (SM), the mechanism proposed to explain the neutrino
mass usually also induces interesting phenomena with lepton flavor violation
(LFV), such as $\mu\to e\gamma$, $\mu\to 3e$, $\mu-e$ conversion in nucleus,
and muonium-antimuonium oscillation, that are highly suppressed in the SM.

Due to their high sensitivities to new physics effects, many ongoing or
planning experiments are designed to achieve an unprecedented precision to
study the LFV processes. For instance, the branching ratio (BR) for the
$\mu\to e\gamma$ decay in MEG II experiment can reach the sensitivity of
$6\times 10^{-14}$ MEGII:2018kmf ; $BR(\mu\to 3e)\lesssim 10^{-16}$ is
expected in Mu3e experiment Perez:2021gnr ; the $\mu-e$ conversion rate in
COMET at J-PARC COMET:2018auw  and Mu2e at FNAL Diociaiuti:2020yvo  can reach
the level of $10^{-17}-10^{-18}$, whereas the PRISM/PRIME experiment will push
the conversion rate down below $10^{-18}$ Barlow:2011zza ; and the probability
of muonium-antimuonium conversion in the new generation experiment with the
sensitivity of ${\cal O}(10^{-14})$ is proposed by the MACE collaboration at
CSNC MACE:CSNC .

In order to explain the observed neutrino mass and DM relic density and to
generate testable LFV processes in a scotogenic model, we study a gauged
$U(1)_{X}$ extension of the SM, where only the newly introduced particles
carry the $U(1)_{X}$ charges. The neutrino mass then arises from the quantum
corrections where the particles carrying the $U(1)_{X}$ charges are running in
the loops. The DM candidate is stable if the $U(1)_{X}$ gauge symmetry is
unbroken. However, the unbroken $U(1)_{X}$ would normally leads to a massless
dark photon ($Z^{\prime}$) and, thus, the DM relic density may not be
correctly reproduced through the resonant $Z^{\prime}$ production when we
focus on the minimal extension of the SM. To guarantee the coexistence of a
massive $Z^{\prime}$ gauge boson and gauge invariance, we consider the
Stueckelberg mechanism by introducing a Stueckelberg scalar field instead of
the mechanism of spontaneous symmetry breaking.

Using the scotogenic model to radiatively generate the Dirac neutrino mass in
the Stueckelberg mechanism was studied in Ref. Leite:2020wjl . In this work,
using a different framework, we construct a model that can induce the Majorana
neutrino mass. To retain the $U(1)_{X}$ gauge anomaly-free, we use a vector-
like singlet Dirac fermion ($N_{L,R}$) instead of the right-handed singlet
Majorana fermion in the Ma model. Utilizing the unbroken $U(1)_{X}$ gauge
symmetry to stabilize DM is done at the cost of introducing two inert-Higgs
doublets. Because the two inert-Higgs doublets ($\eta_{1,2}$) carry the lepton
number and different $U(1)_{X}$ charges, lepton number violation (LNV) for
generating the Majorana neutrino mass arises from the quartic coupling
$(H^{\dagger}\eta_{1})(H^{\dagger}\eta_{2})$ in the scalar potential. In
addition, since the Yukawa couplings to the SM lepton can have different
structures, such as $\overline{L}N_{R}\tilde{\eta}_{1}$ and
$\overline{L^{C}}N_{L}\tilde{\eta}_{2}$, it is found that the neutrino data
can be fitted with only one generation $N_{L,R}$.

When the neutrino oscillation data and the upper limits on LFV processes are
taken into account, it is found that the $\mu-e$ conversion rate and the
$\mu\to 3e$ decay will exclude most of the parameter space that lead to
$BR(\mu\to e\gamma)\sim O(10^{-15})$, where the former is mediated by the
photon-penguin diagrams and the letter is dominated by the box diagrams in the
model, respectively. In other words, both experiments will give the strongest
constraints on new physics for the $\mu\to e$ LFV processes. In addition, the
resulting $BR(\tau\to\mu\gamma)$ and $BR(\tau\to 3\mu)$ can reach the
significant levels of $O(10^{-9})$ and $O(10^{-10})$, respectively, as
expected to probe by Belle II experiment Belle-II:2018jsg .

Two sources can contribute to the DM relic density in the model. One is
through the Yukawa couplings and the other arises from the gauge coupling of
the kinetic mixing. The kinetic mixing between $Z^{\prime}$ and the gauge
field of $U(1)_{Y}$ can be induced via quantum corrections; thus, the
$Z^{\prime}$ gauge boson can couple to the SM particles. The singlet fermion
$N$, which is the DM candidate in the model, can then annihilate into the SM
particles through the $Z^{\prime}$-mediated s-channel process $NN\to
Z^{\prime}\to F\bar{F}$ and the $N$-mediated t-channel one $NN\to
Z^{\prime}(\to F\bar{F})Z^{\prime}(\to F^{\prime}\bar{F}^{\prime})$, where
$F(F^{\prime})$ denotes possible SM particles. We will demonstrate that in
addition to the Yukawa coupling scenario, the induced $Z^{\prime}$ gauge
coupling one can also explain the DM relic density $\Omega_{\rm DM}h^{2}$.
Moreover, the $Z^{\prime}$ couplings to quarks will contribute to the DM-
nucleon scattering. It is found that the resulting DM-nucleon scattering cross
section is under the current upper bound of direct search of DM from the
XENON1T experiment XENON:2018voc .

The structure of this paper is organized as follows. We introduce the model
and derive the relevant scalar couplings, Yukawa couplings, and $Z^{\prime}$
gauge couplings in Sec. II. In Sec. III, we formulate the phenomena for
neutrino mass and LFV processes. Sec. IV contains the detailed numerical
analysis. Finally, we summarize our findings of the study in Sec. V.

##  II The Model

To obtain the Majorana mass for neutrinos through radiative corrections in a
scotogenic $U(1)_{X}$ Stueckelberg gauge model, we find that a minimal
extension of the SM is to include one singlet vector-like fermion $N$ and two
new inert scalar doublets $\eta_{1,2}$, where their representations and charge
assignments are given in Table 1. Note that here we use the convention that
the electromagnetic charge $Q=T_{L}^{3}+Y/2$. To break the lepton number in
the scalar sector, $\eta_{1,2}$ need to carry one unit of lepton number while
the singlet vector-like fermion has no lepton number. As a result, the new
particles are $R$-parity odd and the SM particles are $R$-parity even, where
the $R$-parity quantum number is defined as $R_{p}=(-1)^{3B+L+2S}$, with $B$,
$L$, and $S$ bring the baryon number, lepton number, and spin of the particle,
respectively. Due to the odd $R$-parity property, $N$ and the neutral
components in $\eta_{1,2}$ can be DM candidates. Based upon the charge
assignments, in the following we discuss the relevant Yukawa and gauge
couplings for the phenomenological analysis.

Table 1:  Representations and charged assignments of new particles.  |  $SU(2)_{L}$ |  $U(1)_{Y}$ |  $U(1)_{X}$ |  Lepton # | $R_{p}$
---|---|---|---|---|---
$N$ | 1 | 0 | $Q_{X}$ | 0 | $-1$
$\eta_{1}$ | 2 | 1 | $Q_{X}$ | 1 | $-1$
$\eta_{2}$ | 2 | 1 | $-Q_{X}$ | 1 | $-1$

###  II.1 Scalar masses and mixings, and Yukawa couplings

Apart from the Yukawa couplings, the most important effect to generate the
Majorana mass from the scotogenic mechanism is the appearance of LNV term in
the scalar potential, which dictates the scalar masses, couplings, and
mixings. To examine these effects in the model, we write the scalar potential
for the SM Higgs $H$ and $\eta_{1,2}$ as:

$\displaystyle\begin{split}V=&V_{\rm SM}+V(H,\eta_{1},\eta_{2})~{},\\\ V_{\rm
SM}=&-\mu^{2}_{H}H^{\dagger}H+\lambda_{H}(H^{\dagger}H)^{2}~{},\\\
V(H,\eta_{1},\eta_{2})=&\mu^{2}_{1}\eta^{\dagger}_{1}\eta_{1}+\mu^{2}_{2}\eta^{\dagger}_{2}\eta_{2}+\lambda_{1}(\eta^{\dagger}_{1}\eta_{1})^{2}+\lambda_{2}(\eta^{\dagger}_{2}\eta_{2})^{2}+\lambda_{3}(\eta^{\dagger}_{1}\eta_{1})(\eta^{\dagger}_{2}\eta_{2})\\\
&+\lambda_{4}(\eta^{\dagger}_{1}\eta_{2})(\eta^{\dagger}_{2}\eta_{1})+\left(\lambda_{5}(H^{\dagger}\eta_{1})(H^{\dagger}\eta_{2})+\mbox{H.c.}\right)+\lambda_{6}(H^{\dagger}\eta_{1})(\eta^{\dagger}_{1}H)\\\
&+\lambda_{7}(H^{\dagger}\eta_{2})(\eta^{\dagger}_{2}H)+\lambda_{8}(H^{\dagger}H)(\eta^{\dagger}_{1}\eta_{1})+\lambda_{9}(H^{\dagger}H)(\eta^{\dagger}_{2}\eta_{2})~{}.\end{split}$
(1)

It can be seen that the only non-self-Hermitian term comes from the
$\lambda_{5}$ term, which violates the lepton number by two units and plays an
important role on the radiative generation of the Majorana neutrino mass. The
tiny neutrino mass can be achieved when $\lambda_{5}\ll 1$, which is the same
as that shown in the Ma model Ma:2006km . For spontaneously breaking the
electroweak gauge symmetry, we take $\mu^{2}_{H},\lambda_{H}>0$ as in the SM.
The masses of $\eta_{1,2}$ can be irrelevant to the electroweak symmetry
breaking, and we thus require $\mu^{2}_{1,2}(\lambda_{1,2})>0$. To preserve
the $U(1)_{X}$ and $R_{p}$ symmetries, the vacuum expectation values (VEVs) of
$\eta_{1,2}$ have to vanish. We therefore parametrize the components of the
three doublet scalars as:

$H=\left(\begin{array}[]{c}G^{+}\\\ \frac{1}{\sqrt{2}}(v+h+iG^{0})\\\
\end{array}\right)~{},~{}~{}\eta_{j}=\left(\begin{array}[]{c}\eta^{+}_{j}\\\
\frac{1}{\sqrt{2}}(s_{j}+ia_{j})\\\ \end{array}\right)~{},$ (2)

where $G^{\pm,0}$ are the Goldstone bosons, $v$ is the VEV of $H$, and $h$ is
the SM Higgs boson.

Since $\eta_{1}$ and $\eta_{2}$ carry different $U(1)_{X}$ charges, the
charged scalars $\eta^{\pm}_{1}$ and $\eta^{\pm}_{2}$ do not mix and their
masses are respectively obtained as:

$\displaystyle\begin{split}m^{2}_{\eta^{\pm}_{1}}&=\mu^{2}_{1}+\frac{\lambda_{8}v^{2}}{2}~{},\\\
m^{2}_{\eta^{\pm}_{2}}&=\mu^{2}_{2}+\frac{\lambda_{9}v^{2}}{2}~{}.\end{split}$
(3)

Unlike the charged scalars, the neutral components of $\eta_{1,2}$ can mix via
the $\lambda_{5}$ term. According to the scalar potential with Eq. (2), the
mass-square matrices for $(s_{1},s_{2})$ and $(a_{1},a_{2})$ are respectively:

$m^{2}_{S}=\left(\begin{array}[]{cc}m^{2}_{s_{1}}&m^{2}_{12}\\\
m^{2}_{12}&m^{2}_{s_{2}}\\\
\end{array}\right)~{},~{}~{}m^{2}_{A}=\left(\begin{array}[]{cc}m^{2}_{s_{1}}&-m^{2}_{12}\\\
-m^{2}_{12}&m^{2}_{s_{2}}\\\ \end{array}\right)~{},$ (4)

with

$\displaystyle m^{2}_{s_{1}}$
$\displaystyle=\mu^{2}_{1}+\frac{v^{2}}{2}(\lambda_{6}+\lambda_{8})~{},$
$\displaystyle m^{2}_{s_{2}}$
$\displaystyle=\mu^{2}_{2}+\frac{v^{2}}{2}\left(\lambda_{7}+\lambda_{9}\right)~{},$
$\displaystyle m^{2}_{12}$ $\displaystyle=\frac{v^{2}}{2}\lambda_{5}~{}.$ (5)

Each of the two $2\times 2$ mass-square matrices can be diagonalized by the
corresponding orthogonal matrix $O(\theta_{\xi})$ ($\xi=S,A$) through
$O(\theta_{\xi})m^{2}_{\xi}O^{T}(\theta_{\xi})$, where the $O(\theta_{\xi})$
matrix can be parametrized as:

$O(\theta_{\xi})=\left(\begin{array}[]{cc}\cos\theta_{\xi}&\sin\theta_{\xi}\\\
-\sin\theta_{\xi}&\cos\theta_{\xi}\\\ \end{array}\right)~{}.$ (6)

Since the matrix elements in $m^{2}_{A}$ are the same as those in $m^{2}_{S}$
except the sign change in the off-diagonal elements, we therefore take
$\theta_{S}=-\theta_{A}\equiv\theta$. The eigenvalues for $m^{2}_{S}$ are
found to be:

$m^{2}_{S_{1,2}}=\frac{1}{2}\left[m^{2}_{s_{1}}+m^{2}_{s_{2}}\pm\sqrt{(m^{2}_{s_{2}}-m^{2}_{s_{1}})^{2}-4(m^{4}_{12})}\right]~{}.$
(7)

For the physical pseudoscalars $A_{1,2}$, we have
$m^{2}_{A_{1(2)}}=m^{2}_{S_{1(2)}}$. The mixing angle $\theta$ is given by:

$\sin(2\theta)=-\frac{\lambda_{5}v^{2}}{m^{2}_{S_{2}}-m^{2}_{S_{1}}}~{}.$ (8)

Because $S_{i}$ and $A_{i}$ are degenerate in mass, if one of them is the DM
candidate, the large gauge interaction $A_{i}-S_{i}-Z$ will render too large a
DM-nucleon scattering cross section. Hence, the possibility of using a scalar
as the DM candidate in this model is excluded by the direct detection
experiments. Instead, the singlet vector-like fermion $N$ becomes a promising
DM candidate in the model. Since $\lambda_{5}$ term violates the lepton number
and eventually leads to the Majorana mass, its value has to be sufficiently
small, $\lambda_{5}\ll 1$ Ma:2006km , as alluded to earlier. As a result, the
off-diagonal mass matrix element $|m_{12}^{2}|$ is suppressed and the mixing
angle $\theta\ll 1$. In order to make the Yukawa couplings sufficiently large
so that the LFV processes can be possibly detectable in the ongoing and
planning experiments, we follow the approach in Toma:2013zsa  and take
$\lambda_{5}=10^{-9}$.

According to the $U(1)_{X}$ charges listed in Table 1, the Yukawa couplings of
$\eta_{1,2}$ and $N$ to the SM leptons are given by:

$\displaystyle-{\cal L}_{Y}$ $\displaystyle=\overline{L}\,{\bf
y}_{1}\,\tilde{\eta}_{1}N_{R}+\overline{L}\,{\bf
y}_{2}\,\tilde{\eta}_{2}N^{C}_{L}+m_{N}\overline{N}_{L}N_{R}+\mbox{H.c.}~{},$
(9)

where the flavor indices are suppressed, $L^{T}=(\nu_{\ell},\ell)$ denotes the
left-handed lepton doublet in the SM,
${\tilde{\eta}_{j}}=i\tau_{2}\eta^{*}_{j}$ with $\tau_{2}$ being a Pauli
matrix, $N^{C}=C\gamma_{0}N^{*}$ is the charge conjugation of $N$, and $m_{N}$
is the mass of $N$. Although ${\bf y}_{1,2}$ can be generally complex, we can
rotate away the phase of one of them by redefining the phases of the complex
lepton doublet $L$. In our following analysis, we take the convention that
${\bf y}_{1}$ is real and ${\bf y}_{2}$ is complex, as parametrized by
$y_{2k}=|y_{2k}|e^{i\phi_{k}}$ ($k=1,2,3$). Using Eqs. (2) and
$(\ref{eq:omatrix})$, the Yukawa couplings in Eq. (9) can be decomposed as:

$\displaystyle\begin{split}-{\cal L}_{Y}=&-\overline{\ell}_{L}{\bf
y}_{1}N_{R}\eta^{-}_{1}-\overline{\ell}_{L}{\bf
y}_{2}N^{C}_{L}\eta^{-}_{2}+m_{N}\overline{N}_{L}N_{R}\\\
&+\frac{1}{\sqrt{2}}\overline{\nu}_{\ell L}{\bf
y}_{1}N_{R}\left[c_{\theta}(S_{1}-iA_{1})-s_{\theta}(S_{2}+iA_{2})\right]\\\
&+\frac{1}{\sqrt{2}}\overline{\nu}_{\ell L}{\bf
y}_{2}N^{C}_{L}\left[s_{\theta}(S_{1}+iA_{1})+c_{\theta}(S_{2}-iA_{2})\right]+\mbox{H.c.}~{},\end{split}$
(10)

where we denote $c_{\theta}~{}(s_{\theta})\equiv\cos\theta~{}(\sin\theta)$. We
note that if we take the charged lepton $\ell_{L}$ as the mass eigenstate, in
general, the $\eta^{\pm}_{1,2}$ Yukawa couplings are modified as
$V^{\ell}_{L}{\bf y}_{i}$, where $V^{\ell}_{L}$ is a unitary matrix introduced
for diagonalizing the charged lepton mass matrix. Therefore, the Yukawa
couplings of $\eta^{\pm}_{i}$ are generally different from those of $S_{i}$
and $A_{i}$. Since we know nothing about the information of $V^{\ell}_{L}$ in
the SM, we adopt $V^{\ell}_{L}=1$ to reduce the number of free parameters in
the study. As a result, the effects contributing to the neutrino mass through
the Yukawa couplings of $S_{i}$ and $A_{i}$ now have a strong correlation to
the LFV processes that are mediated by $\eta^{\pm}_{i}$.

###  II.2 $U(1)_{X}$ Stueckelberg gauge couplings

The Lagrangian invariant under the $SU(2)_{L}\times U(1)_{Y}\times U(1)_{X}$
gauge symmetry with the Stueckelberg gauge field included is given by:

$\displaystyle\begin{split}{\cal L}_{\rm
kin}=&-\frac{1}{4}\vec{W}_{\mu\nu}\cdot\vec{W}^{\mu\nu}-\frac{1}{4}\hat{B}_{\mu\nu}\hat{B}^{\mu\nu}-\frac{1}{4}\hat{Z^{\prime}}_{\mu\nu}\hat{Z^{\prime}}^{\mu\nu}+\frac{m^{2}_{Z^{\prime}}}{2}\left(\hat{Z^{\prime}}^{\mu}-\frac{1}{m_{Z^{\prime}}}\partial^{\mu}B\right)^{2}\\\
&-\frac{1}{2\beta}\left(\partial_{\mu}\hat{Z^{\prime}}^{\mu}+m_{Z^{\prime}}\beta
B\right)^{2}+\overline{N}i\not{D}N+(D_{\mu}\eta_{1})^{\dagger}(D^{\mu}\eta_{1})+(D_{\mu}\eta_{2})^{\dagger}(D^{\mu}\eta_{2})~{},\end{split}$
(11)

where $W^{i}_{\mu\nu}$, $\hat{B}_{\mu\nu}$, and $\hat{Z^{\prime}}_{\mu\nu}$
are the gauge field stress tensors associated with the $SU(2)_{L}$,
$U(1)_{Y}$, and $U(1)_{X}$ gauge groups, respectively, $B$ field is the
Stueckelberg scalar field, $\beta$ is the gauge fixing parameter, and
$m_{Z^{\prime}}$ is the mass of $U(1)_{X}$ Stueckelberg gauge field. We note
that the kinetic mixing term $\hat{B}_{\mu\nu}\hat{Z^{\prime}}^{\mu\nu}$ can
be rotated away by redefining the gauge fields $\hat{B}_{\mu}$ and
$\hat{Z^{\prime}}_{\mu}$. However, we will show later that the mixings in
$\gamma$-$Z^{\prime}$ and $Z$-$Z^{\prime}$ can be induced through radiative
corrections.

The covariant derivatives on $N$ and $\eta_{i}$ are expressed as:

$\displaystyle\begin{split}D_{\mu}N&=(\partial_{\mu}+ig_{X}q_{N}\hat{Z^{\prime}}_{\mu})N~{},\\\
D_{\mu}\eta_{i}&=\left(\partial_{\mu}+i\frac{g}{2}\vec{\tau}\cdot\vec{W}_{\mu}+i\frac{g^{\prime}}{2}\hat{B}_{\mu}+ig_{X}q_{\eta_{i}}\hat{Z^{\prime}}_{\mu}\right)\eta_{i}~{},\end{split}$
(12)

where $q_{N(\eta_{i})}$ denotes the $U(1)_{X}$ charge of $N(\eta_{i})$, and
$g$ and $g^{\prime}$ are the gauge couplings of $SU(2)_{L}$ and $U(1)_{Y}$,
respectively. The $NNZ^{\prime}$ interaction can be immediately obtained as:

${\cal
L}_{NNZ^{\prime}}=-g_{X}q_{N}\overline{N}\gamma^{\mu}N\hat{Z^{\prime}}_{\mu}~{}.$
(13)

Although the quartic gauge couplings of $\eta_{i}$ to gauge bosons $W^{\pm}$,
$Z$, $\gamma$, and $\hat{Z^{\prime}}$ can be derived from the kinetic terms of
$\eta_{i}$, these couplings are irrelevant to our study here and are not
presented explicitly. The various trilinear gauge couplings of $\eta_{i}$ from
$(D_{\mu}\eta_{i})^{\dagger}(D^{\mu}\eta_{i})$ can be summarized as:

$\displaystyle\begin{split}{\cal L}_{\rm
kin}\supset&~{}i\frac{g}{\sqrt{2}}\left\\{\left[\partial_{\mu}\eta^{-}_{i}(S_{i}+iA_{i})-\eta^{-}_{i}(\partial_{\mu}S_{i}+i\partial_{\mu}A_{i})\right]W^{+\mu}+\mbox{H.c.}\right\\}\\\
&+\left(\frac{g}{2c_{W}}Z^{\mu}-g_{X}q_{X}\hat{Z^{\prime}}^{\mu}\right)\left[(\partial_{\mu}S_{i})A_{i}-S_{i}\partial_{\mu}A_{i}\right]\\\
&+i\left(eA^{\mu}+\frac{1-2s^{2}_{W}}{2c_{W}}gZ^{\mu}+g_{X}q_{X}\hat{Z^{\prime}}^{\mu}\right)\left[(\partial_{\mu}\eta^{-}_{i})\eta^{+}_{i}-\eta^{-}_{i}\partial_{\mu}\eta^{+}_{i}\right]~{}.\end{split}$
(14)

The $W^{\pm}$, $Z$, and $\gamma$ gauge bosons and the Weinberg angle
$\theta_{W}$ are defined through the relations:

$\displaystyle\begin{split}W^{\pm}&=\frac{1}{\sqrt{2}}(W^{1}_{\mu}\mp
W^{2}_{\mu})~{},\\\ A_{\mu}&=c_{W}\hat{B}_{\mu}+s_{W}W^{3}_{\mu}~{},\\\
Z_{\mu}&=-s_{W}\hat{B}_{\mu}+c_{W}W^{3}_{\mu}~{},\end{split}$ (15)

where $c_{W}~{}(s_{W})\equiv\cos\theta_{W}~{}(\sin\theta_{W})$ and
$e=gs_{W}=g^{\prime}c_{W}$ are used.

In the study of LFV processes, we need to know the photon and $Z$-boson gauge
couplings to the SM particles. Therefore, the relevant photon and $Z$ gauge
couplings are given by:

${\cal
L}_{Z}=-eQ_{f}\overline{f}\gamma_{\mu}fA^{\mu}-\frac{g}{2c_{W}}\overline{f}\left(C^{f}_{R}P_{R}+C^{f}_{L}P_{L}\right)\gamma_{\mu}fZ^{\mu}~{},$
(16)

with

$C^{f}_{R}=-2Q_{f}s^{2}_{W}~{},~{}~{}C^{f}_{L}=2T^{3}_{f}-2Q_{f}s^{2}_{W}~{},$
(17)

where $Q_{f}$ is the electric charge of the fermion $f$, and $T^{3}_{f}$ is
its third component of the weak isospin.

##  III Neutrino mass, lepton-flavor violation, kinetic mixing of gauge boson

Based on the introduced interactions, in this section we investigate the new
physics effects on the rare lepton flavor related processes, such as neutrino
mass, $\ell_{i}\to\ell_{j}\gamma$, $\ell_{i}\to 3\ell_{j}$, $\mu-e$ conversion
in nuclei, and muonium-antimuonium oscillation. Note that throughout this
section, we assume that the new physics effects dominate and ignore the SM
contributions.

###  III.1 Scotogenic neutrino mass matrix

Figure 1: Sketched Feynman diagram mediated by $Z_{2}$-odd particles for the
radiative neutrino mass.

To obtain the Majorana neutrino mass in the model, we need both Yukawa
couplings ${\bf y}_{1}$ and ${\bf y}_{2}$, which can ensure that the structure
of Majorana mass term $L^{T}CL$ can be achieved. In addition, since
$\eta_{1,2}$ have no VEVs and the lepton number is preserved in the ground
state, we are left with the choice of using the lepton number-violating
interaction $\lambda_{5}(H^{\dagger}\eta_{1})(H^{\dagger}\eta_{2})$ to
generate the Weinberg’s dimension-5 operator $LHHL$ Weinberg:1979sa . For the
purpose of clarity, we show a representative one-loop Feynman diagram in Fig.
1. According to the Yukawa couplings in Eq. (10), the neutrino mass matrix
elements can be obtained as:

$m^{\nu}_{ij}=\frac{\sin(2\theta)}{16\pi^{2}}\overline{y}_{ij}m_{N}\left[\frac{m^{2}_{S_{1}}}{m^{2}_{S_{1}}-m^{2}_{N}}\ln\left(\frac{m^{2}_{S_{1}}}{m^{2}_{N}}\right)-\frac{m^{2}_{S_{2}}}{m^{2}_{S_{2}}-m^{2}_{N}}\ln\left(\frac{m^{2}_{S_{2}}}{m^{2}_{N}}\right)\right]~{},$
(18)

where we have included $A_{i}$ contributions, $m_{A_{i}}=m_{S_{i}}$ is used,
and the symmetric Yukawa couplings $\overline{y}_{ij}$ in flavor indices are
defined as:

$\overline{y}_{ij}=y^{*}_{1i}y^{*}_{2j}+y^{*}_{2i}y^{*}_{1j}~{}.$ (19)

In the Ma model, the involved Yukawa couplings for $m^{\nu}_{ij}$ appear in
the combination of $y_{i}y_{j}$; thus, the mass matrix elements have strong
correlations. It needs at least two singlet right-handed fermions to explain
the neutrino data. In our model, due to the introduction of one more inert
Higgs doublet $\eta_{2}$, the induced neutrino mass matrix elements appear in
the combination of $y_{1i}y_{2j}+y_{2i}y_{1j}$ and have less correlations
among the matrix elements. It is for this reason that the neutrino data can be
accommodated using only one singlet fermion.

The symmetric mass matrix $m^{\nu}$ can be diagonalized through $m^{\nu}_{\rm
dia}=U^{\dagger}m^{\nu}U^{*}$, where $m^{\nu}_{\rm dia}={\rm
dia}(m_{1},m_{2},m_{3})$ is the diagonal mass matrix. The unitary matrix $U$
can be parametrized as:

$U=\left(\begin{array}[]{ccc}c_{12}c_{13}&s_{12}c_{13}&s_{13}e^{-i\delta_{\rm
CP}}\\\ -s_{12}c_{23}-c_{12}s_{13}s_{23}e^{i\delta_{\rm
CP}}&c_{12}c_{23}-s_{12}s_{13}s_{23}e^{i\delta_{\rm CP}}&c_{13}s_{23}\\\
s_{12}s_{23}-c_{12}s_{13}c_{23}e^{i\delta_{\rm
CP}}&-c_{12}s_{23}-s_{12}s_{13}c_{23}e^{i\delta_{\rm
CP}}&c_{13}c_{23}\end{array}\right)\left(\begin{array}[]{ccc}e^{i\alpha_{1}}&&\\\
&e^{i\alpha_{2}}&\\\ &&1\end{array}\right)~{},$

where $c_{ij}~{}(s_{ij})\equiv\cos\theta_{ij}~{}(\sin\theta_{ij})$,
$\delta_{\rm CP}$ is the Dirac CP-violating phase, and $\alpha_{1,2}$ are the
Majorana phases. Consequently, the theoretical mass matrix obtained in Eq.
(18) can be determined by the observables from the neutrino oscillations, and
their relations are given by:

$\displaystyle\begin{split}(m^{\nu})_{11}=&~{}m_{1}c^{2}_{12}c^{2}_{13}e^{i2\alpha_{1}}+m_{2}s^{2}_{12}c^{2}_{13}e^{i2\alpha_{2}}+m_{3}s^{2}_{13}e^{-2i\delta_{\rm
CP}}~{},\\\
(m^{\nu})_{12}=&~{}m_{1}e^{2i\alpha_{1}}c_{12}c_{13}\left(-s_{12}c_{23}-e^{i\delta_{\rm
CP}}c_{12}s_{13}s_{23}\right)\\\
&+m_{2}e^{2i\alpha_{2}}s_{12}c_{13}\left(c_{12}c_{23}-e^{i\delta_{\rm
CP}}s_{12}s_{13}s_{23}\right)+m_{3}e^{-i\delta_{\rm
CP}}s_{13}c_{13}s_{23}~{},\\\
(m^{\nu})_{13}=&~{}m_{1}e^{2i\alpha_{1}}c_{12}c_{13}\left(s_{12}s_{23}-e^{i\delta_{\rm
CP}}c_{12}s_{13}s_{23}\right)\\\
&+m_{2}e^{2i\alpha_{2}}s_{12}c_{13}\left(-c_{12}s_{23}-e^{i\delta_{\rm
CP}}s_{12}s_{13}c_{23}\right)+m_{3}e^{-i\delta_{\rm
CP}}s_{13}c_{13}c_{23}~{},\\\
(m^{\nu})_{22}=&~{}m_{1}e^{2i\alpha_{1}}\left(-s_{12}c_{23}-e^{i\delta_{\rm
CP}}c_{12}s_{13}s_{23}\right)^{2}\\\
&+m_{2}e^{2i\alpha_{2}}\left(c_{12}c_{23}-e^{i\delta_{\rm
CP}}s_{12}s_{13}s_{23}\right)^{2}+m_{3}c^{2}_{13}s^{2}_{23}~{},\\\
(m^{\nu})_{23}=&~{}m_{1}e^{2i\alpha_{1}}\left(s_{12}s_{23}-e^{i\delta_{\rm
CP}}c_{12}s_{13}c_{23}\right)\left(-s_{12}c_{23}-e^{i\delta_{\rm
CP}}c_{12}s_{13}s_{23}\right)\\\
&+m_{2}e^{2i\alpha_{2}}\left(-c_{12}s_{23}-e^{i\delta_{\rm
CP}}s_{12}s_{13}c_{23}\right)\left(c_{12}c_{23}-e^{i\delta_{\rm
CP}}s_{12}s_{13}s_{23}\right)\\\ &+m_{3}e^{-i\delta_{\rm
CP}}c^{2}_{13}s_{23}c_{23}~{},\\\
(m^{\nu})_{33}=&~{}m_{1}e^{2i\alpha_{1}}\left(s_{12}s_{23}-e^{i\delta_{\rm
CP}}c_{12}s_{13}c_{23}\right)^{2}\\\
&+m_{2}e^{2i\alpha_{2}}\left(-c_{12}s_{23}-e^{i\delta_{\rm
CP}}s_{12}s_{13}c_{23}\right)^{2}+m_{3}c^{2}_{13}c^{2}_{23}~{}.\end{split}$
(20)

###  III.2  Radiative lepton decay $\ell_{i}\to\ell_{j}\gamma$

Among the LFV processes, the most constraining process is the radiative muon
decay $\mu\to e\gamma$. It can be the most important process to discover LFV
due to the high sensitivity to the new physics effects. In order to study
other radiative lepton decays, such as $\tau\to(e,\mu)\gamma$, in the
following, we calculate the branching ratios for the
$\ell_{i}\to\ell_{j}\gamma$ decays, where $\ell_{j}$ is the lighter lepton and
its mass is neglected unless otherwise stated.

Figure 2: One-loop Feynman diagrams for $\ell_{i}\to\ell_{j}\gamma^{(*)}$

The radiative LFV process in the model arises from the
$\eta^{\pm}_{i}$-mediated diagrams, as shown in Fig. 2. Using the Yukawa
couplings in Eq. (10), the loop-induced effective Lagrangian for
$\ell_{i}\to\ell_{j}\gamma^{(*)}$ can be expressed as:

${\cal
L}^{\gamma}_{\ell_{i}\ell_{j}\gamma}=ek^{2}C^{\gamma}_{1ji}\,\overline{\ell}_{j}\gamma_{\mu}P_{L}\ell_{i}A^{\mu}-\frac{em_{\ell_{i}}}{2}C^{\gamma}_{2ji}\,\overline{\ell}_{j}\sigma_{\mu\nu}P_{R}\ell_{i}F^{\mu\nu}~{},$
(21)

where the loop induced Wilson coefficients are obtained as:

$\displaystyle\begin{split}C^{\gamma}_{1ji}&=\sum^{2}_{\alpha=1}\frac{y^{*}_{\alpha
i}y_{\alpha
j}}{(4\pi)^{2}m^{2}_{\eta^{\pm}_{\alpha}}}J^{\gamma}_{1}\left(\frac{m^{2}_{N}}{m^{2}_{\eta^{\pm}_{\alpha}}}\right)~{},\\\
C^{\gamma}_{2ji}&=\sum^{2}_{\alpha=1}\frac{y^{*}_{\alpha i}y_{\alpha
j}}{(4\pi)^{2}m^{2}_{\eta^{\pm}_{\alpha}}}J^{\gamma}_{2}\left(\frac{m^{2}_{N}}{m^{2}_{\eta^{\pm}_{\alpha}}}\right)~{},\end{split}$
(22)

and the loop integral functions are defined by

$\displaystyle\begin{split}J^{\gamma}_{1}(x)&=\frac{1}{36(1-x)^{4}}\left(2-9x+18x^{2}-11x^{3}+6x^{3}\ln
x\right)~{},\\\
J^{\gamma}_{2}(x)&=\frac{1}{12(1-x)^{4}}\left(1-6x+3x^{2}+2x^{3}-6x^{2}\ln
x\right)~{}.\end{split}$ (23)

The emitted photon can be on-shell ($k^{2}=0$) or off-shell ($k^{2}\neq 0$),
where the latter can be used to study the $\ell_{i}\to 3\ell_{j}$ process, in
which a $\ell_{j}^{+}\ell_{j}^{-}$ pair is produced by the off-shell photon.
(The contribution of the $Z$-mediated diagrams is found to be small and will
be neglected, as discussed in the next section.) Using the effective
Lagrangian in Eq. (21), the branching ratio for $\ell_{i}\to\ell_{j}\gamma$
can be obtained as Toma:2013zsa :

$BR(\ell_{i}\to\ell_{j}\gamma)=\frac{3(4\pi)^{3}\alpha_{\rm
em}}{4G^{2}_{F}}\left|C^{\gamma}_{2ji}\right|^{2}BR(\ell_{i}\to\ell_{j}\overline{\nu}_{j}\nu_{i})~{},$
(24)

with $\alpha_{\rm em}=e^{2}/(4\pi)$.

In addition to the radiative LFV decays, the dipole operator in Eq. (21) can
contribute to the lepton anomalous magnetic dipole moment ($g-2$) when we
replace $\ell_{j}$ by $\ell_{i}$. As a result, the lepton $(g-2)$ can be
directly obtained as:

$a_{\ell}\equiv\left(\frac{g-2}{2}\right)_{\ell}=-m^{2}_{\ell}C^{\gamma}_{2\ell\ell}~{}.$
(25)

###  III.3 $\ell_{i}\to 3\ell_{j}$ from penguin and box diagrams

We now discuss the three-body lepton flavor-changing decays. For simplicity,
we only concentrate on the $\ell_{i}\to 3\ell_{j}$ decay. In the model, the
LFV processes $\ell_{i}\to 3\ell_{j}$ can arise from photon-penguin,
$Z$-penguin, and box diagrams. In the following, we discuss their
contributions in sequence.

Since the effective interaction for $\ell_{i}\to\ell_{j}\gamma^{*}$ has been
given in Eq. (21), the transition amplitude with the assigned momenta for the
$\ell_{i}(p)\to\ell_{j}(p_{1})\ell^{+}(p_{2})\ell^{-}(p_{3})$ decay can be
written as:

$\displaystyle\begin{split}{\cal
M}_{\gamma}=&~{}e^{2}C^{\gamma}_{1ji}\overline{u}(p_{1})\gamma_{\mu}u(p)\,\overline{u}(p_{3})\gamma^{\mu}v(p_{2})\\\
&+\frac{e^{2}m_{\ell_{i}}}{k^{2}}C^{\gamma}_{2ji}\overline{u}(p_{1})i\sigma_{\mu\nu}k^{\nu}u(p)\,\overline{u}(p_{3})\gamma^{\mu}v(p_{2})-(p_{1}\leftrightarrow
p_{3})~{},\end{split}$ (26)

where the photon coupling to the SM charged lepton shown in Eq. (16) is
applied.

The Feynman diagrams for $\ell_{i}\to\ell_{j}Z^{*}$ are similar to the photon
diagrams shown in Fig. 2, where we use $Z^{*}$ instead of $\gamma^{*}$. Using
the $Z$ gauge couplings to $\eta^{\pm}$ and the SM lepton, given in Eqs. (14)
and (16), the one-loop induced effective interaction for $\ell_{i}\ell_{j}Z$
can be obtained as:

$\displaystyle{\cal L}^{Z}_{\ell_{i}\ell_{j}Z}$
$\displaystyle=\frac{gc_{2W}}{c_{W}}\sum^{2}_{\alpha=1}\frac{m_{\ell_{i}}m_{\ell_{j}}}{(4\pi)^{2}m^{2}_{\eta_{\alpha}}}J^{\gamma}_{2}\left(\frac{m^{2}_{N}}{m^{2}_{\eta_{\alpha}}}\right)\overline{\ell}_{j}\gamma_{\mu}P_{R}\ell_{i}Z^{\mu}~{},$
(27)

where we have dropped the small contributions from $k^{2}$ and $k$ terms
Toma:2013zsa . It can be clearly seen that the loop-induced $Z$ coupling is
proportional to the product of the lepton masses, i.e.,
$m_{\ell_{i}}m_{\ell_{j}}$. We note that because $N$ does not couple to $Z$,
the $Z$ boson cannot be emitted from the fermion line in the loop; therefore,
no chiral enhancement factor, e.g., $m^{2}_{N}/m^{2}_{\eta^{\pm}_{i}}$,
appears in the model. When a necessary mass insertion occurs in the initial
lepton $\ell_{i}$, in order to balance the chirality of the lepton vector
current that couples to the vector gauge boson $Z$, a mass factor is
necessarily inserted in the final lepton $\ell_{j}$. Hence, we get the result
proportional to $m_{\ell_{i}}m_{\ell_{j}}$. Due to the fact that
$m_{\ell_{j}}\ll m_{\ell_{i}}$ and $m_{\ell_{i}}m_{\ell_{j}}/m^{2}_{Z}\ll 1$,
we thus neglect the $Z$-penguin contributions to $\ell_{i}\to 3\ell_{j}$.

Figure 3: Selected box diagrams for $\ell_{i}\to 3\ell_{j}$

The box diagrams contributing to $\ell_{i}\to 3\ell_{j}$ are partly shown in
Fig. 3. In addition to the diagrams mediated by the same $\eta^{\pm}_{1(2)}$,
the other possible diagrams are mediated by both $\eta^{\pm}_{1}$ and
$\eta^{\pm}_{2}$ and involve chirality-flipping factor $m^{2}_{N}$. Using the
Yukawa couplings in Eq. (10), the four-fermion interaction for $\ell_{i}\to
3\ell_{j}$ is obtained as:

${\cal M}_{\rm box}=(C^{\rm box}_{1ji}+C^{\rm
box}_{2ji})\,\overline{\ell}_{j}\gamma_{\mu}P_{L}\ell_{j}\,\overline{\ell}_{j}\gamma^{\mu}P_{L}\ell_{i}~{},$
(28)

where $C^{\rm box}_{1}$ denotes the contributions from diagrams involving the
same $\eta^{\pm}_{1(2)}$ and $C^{\rm box}_{2}$ are from those involving both
$\eta^{\pm}_{1}$ and $\eta^{\pm}_{2}$, and the effective Wilson coefficients
are given by:

$\displaystyle\begin{split}C^{\rm
box}_{1ji}&=-\sum^{2}_{\alpha=1}\frac{y^{*}_{\alpha i}y_{\alpha j}|y_{\alpha
j}|^{2}}{(4\pi)^{2}m^{2}_{\eta_{\alpha}}}J^{\rm
box}_{1}\left(\frac{m^{2}_{N}}{m^{2}_{\eta_{\alpha}}}\right)~{},\\\ C^{\rm
box}_{2ji}&=\frac{y^{*}_{1i}y_{1j}|y_{2j}|^{2}+y^{*}_{2i}y_{2j}|y_{1j}|^{2}}{(4\pi)^{2}m^{2}_{\eta_{2}}}\frac{m^{2}_{N}}{m^{2}_{\eta_{2}}}J^{\rm
box}_{2}\left(\frac{m^{2}_{N}}{m^{2}_{\eta_{2}}},\frac{m^{2}_{\eta_{1}}}{m^{2}_{\eta_{2}}}\right)~{},\end{split}$
(29)

with the loop integral functions given by

$\displaystyle\begin{split}J^{\rm
box}_{1}(x)&=\frac{1+x}{2(1-x)^{2}}+\frac{\ln x}{(1-x)^{3}}~{},\\\ J^{\rm
box}_{2}(x,y)&=-\frac{1}{(1-x)(y-x)}-\frac{(y-x^{2})\ln
x}{(1-x)^{2}(y-x)^{2}}+\frac{y\ln y}{(1-y)(y-x)^{2}}~{}.\end{split}$ (30)

Although $C^{\rm box}_{2ji}$ seems to have a chiral enhancement, its
contribution is in fact somewhat smaller than $C^{\rm box}_{1ji}$ due to the
assumed condition that $m_{N}<m_{\eta_{2}}$.

Combining the contributions from the photon-penguin and box diagrams, the
branching ratio for the $\ell_{i}\to 3\ell_{j}$ decay can be written as
Toma:2013zsa :

$\displaystyle\begin{split}BR(\ell_{i}\to
3\ell_{j})=&~{}\frac{3(4\pi\alpha_{\rm
em})^{2}}{8G^{2}_{F}}\Bigg{[}|C^{\gamma}_{1ji}|^{2}+|C^{\gamma}_{2ji}|^{2}\left(\frac{16}{3}\ln\left(\frac{m_{i}}{m_{j}}\right)-\frac{22}{3}\right)+\frac{|C^{\rm
box}_{ji}|^{2}}{6(4\pi\alpha_{\rm em})^{2}}\\\
&+\frac{2}{3}\left(-6\operatorname{Re}(C^{\gamma}_{1ji}C^{\gamma*}_{2ji})+\frac{1}{(4\pi\alpha_{\rm
em})^{2}}\operatorname{Re}\left((C^{\gamma}_{1ji}-2C^{\gamma}_{2ji})C^{\rm
box*}_{ji}\right)\right)\Bigg{]}\\\ &\times
BR(\ell_{i}\to\ell_{j}\overline{\nu}_{j}\nu_{i})~{},\end{split}$ (31)

where $C^{\rm box}_{ji}=C^{\rm box}_{1ji}+C^{\rm box}_{2ji}$.

###  III.4 $\mu-e$ conversion in nuclei and muonium-antimuonium oscillation

The $\mu-e$ conversion process describes the muon capture by a nucleus through
the process $\mu(A,Z)\to e(A,Z)$. At the quark level, the process can be
represented as $\mu q\to eq$, as induced by the photon- and $Z$-penguin
diagrams in the model. Following the results in Kuno:1999jp ; Arganda:2007jw ,
the conversion rate, which is relative to the muon capture rate, is given by:

$\displaystyle\begin{split}CR(\mu-e,{\rm
Nucleus})=&~{}\frac{p_{e}E_{e}m^{3}_{\mu}G^{2}_{F}\alpha^{3}_{\rm
em}Z^{4}_{\rm eff}F^{2}_{p}}{8\pi^{2}Z\Gamma_{\rm cap}}\\\
&\times\left[\left|A_{+}\left(g^{(0)}_{LV}+g^{(0)}_{LS}\right)+A_{-}\left(g^{(1)}_{LV}+g^{(1)}_{LS}\right)\right|^{2}+\left(R\leftrightarrow
L\right)\right]~{},\end{split}$ (32)

where $p_{e}~{}(E_{e})$ is the momentum (energy) of electron and is taken to
be $m_{\mu}$ in the numerical estimates, $Z_{\rm eff}$ denotes the effective
atomic charge of the nucleus, $F_{p}$ is the nuclear matrix element,
$\Gamma_{\rm cap}$ is the total muon capture rate, $A_{\pm}=Z\pm\tilde{N}$
with $Z~{}(\tilde{N})$ being the proton (neutron) number in the nucleus, and
$g^{(j)}_{\chi K}$ with $\chi=L,R$ and $K=V,S$ are defined by

$\displaystyle\begin{split}g^{(0)}_{\chi
K}&=\frac{1}{2}\sum_{q=u,d,s}\left(g^{q}_{\chi K}G^{(q,p)}_{K}+g^{q}_{\chi
K}G^{(q,n)}_{K}\right)~{},\\\ g^{(1)}_{\chi
K}&=\frac{1}{2}\sum_{q=u,d,s}\left(g^{q}_{\chi K}G^{(q,p)}_{K}-g^{q}_{\chi
K}G^{(q,n)}_{K}\right)~{}.\end{split}$ (33)

Since the $Z$-penguin contribution is negligible, the dominant effect is from
the photon-penguin diagrams. Thus, the only nonzero $g^{q}_{\chi K}$ is:

$g^{q}_{LV}=\frac{\sqrt{2}e^{2}Q_{q}}{G_{F}}(C^{\gamma}_{1e\mu}-C^{\gamma}_{2e\mu})~{}.$
(34)

The nucleon matrix element $G^{(q,N^{\prime})}_{K}$ is defined by $\langle
N^{\prime}|\overline{q}\Gamma_{K}q|N^{\prime}\rangle=G^{(q,N^{\prime})}_{K}\overline{N}^{\prime}\Gamma_{K}N^{\prime}$,
where $\Gamma_{K}=(1,\gamma_{\mu})$ when $K=(S,V)$, respectively. Their values
can be found in Refs. Kuno:1999jp ; Kosmas:2001mv  as:

$\displaystyle
G^{(u,p)}_{V}=G^{(d,n)}_{V}=2~{},~{}~{}G^{(d,p)}_{V}=G^{(u,n)}_{V}=1~{},~{}~{}G^{(s,N^{\prime})}=0~{}.$
(35)

For heavy nuclei, because $A_{-}\ll A_{+}$ and $g^{(1)}_{\chi K}<g^{(0)}_{\chi
K}$, the $\mu-e$ conversion rate in the model can be simplified as:

$CR(\mu-e,{\rm
Nucleus})\approx\frac{p_{e}E_{e}m^{3}_{\mu}G^{2}_{F}\alpha^{3}_{\rm
em}Z^{4}_{\rm eff}F^{2}_{p}}{8\pi^{2}Z\Gamma_{\rm
cap}}A^{2}_{+}\left|g^{(0)}_{LV}\right|^{2}~{}.$ (36)

In addition to the $\mu-e$ conversion and LFV processes, which are $\Delta
L_{\mu}=1$ processes, another interesting process involving $\Delta L_{\mu}=2$
($\Delta L=2$ will be used hereinafter) is the muonium-antimuonium
oscillation, where the muonium is a bound state of $\mu^{+}$ and $e^{-}$,
i.e., $|M_{\mu}\rangle\equiv|\mu^{+}e^{-}\rangle$. Similar to the $\mu\to 3e$
decay, the $M_{\mu}$-$\overline{M}_{\mu}$ mixing matrix element can arise from
the box diagrams, where the Feynman diagrams are similar to Fig. 3 with
$(\ell_{i},\ell_{j})=(\mu,e)$. As in the case of meson oscillations,
$M_{\mu}$-$\overline{M}_{\mu}$ mixing effect can be taken as a perturbative
effect in quantum mechanics and is formulated by Conlin:2020veq

$\left(m-\frac{i}{2}\Gamma\right)_{12}=\frac{1}{2m_{M_{\mu}}}\langle\overline{M}_{\mu}|{\cal
H}_{\rm
eff}|M_{\mu}\rangle+\frac{1}{m_{M_{\mu}}}\sum_{n}\frac{\langle\overline{M}_{\mu}|{\cal
H}_{\rm eff}|n\rangle\langle n|{\cal H}_{\rm
eff}|M_{\mu}\rangle}{m_{M_{\mu}}-E_{n}+i\epsilon}~{},$ (37)

where $m_{12}$ and $\Gamma_{12}$ lead to the mass and lifetime differences
between the two physical states of muonium, and $n$ in the second term denotes
the possible intermediate state. In order to produce the lifetime difference
$\Delta\Gamma$, we need a resonant intermediate state. In the model, the new
particle masses are heavier than the muonium mass $m_{M_{\mu}}$. As such,
$\Delta\Gamma$ by the new effects is negligible, and we concentrate on the
mass difference $\Delta m=m_{1}-m_{2}\approx 2\operatorname{Re}(m_{12})$. From
Eq. (37), the parameter $x\equiv\Delta m/\Gamma_{\mu}$ used to show the
probability of $M_{\mu}\to\overline{M}_{\mu}$ can be written as:

$x\approx\frac{1}{m_{M_{\mu}}\Gamma_{\mu}}\operatorname{Re}(\langle\overline{M}_{\mu}|{\cal
H}_{\rm eff}|M_{\mu}\rangle)~{}.$ (38)

Following the formulation obtained in Ref. Conlin:2020veq , the oscillation
probability is:

$P(M_{\mu}\to\overline{M}_{\mu})=\frac{1}{2}(x^{2}+y^{2})\approx\frac{x^{2}}{2}~{},$
(39)

while the experimental upper limit is:

$P(M_{\mu}\to\overline{M}_{\mu})^{\rm exp}\leq 8.3\times
10^{-11}/S_{B}(B_{0})~{},$ (40)

with $S_{B}(B_{0})=0.75$ Willmann:1998gd ; Conlin:2020veq  taken in this work.
Since the spin-0 para-muonium and spin-1 ortho-muonium are produced in the
experiment Willmann:1998gd , in order to compare the theoretical estimate with
the experimental data, we follow the prescription given in Ref. Conlin:2020veq
and take the spin average by combing both spin-0 and spin-1 muonia as:

$P(M_{\mu}\to\overline{M}_{\mu})^{\rm
exp}=\sum_{i=P,V}\frac{1}{2s_{i}+1}P(M^{i}_{\mu}\to\overline{M}^{i}_{\mu})~{},$
(41)

where $s_{i}$ denotes the spin of the muonium $M^{i}_{\mu}$.

According to the Yukawa couplings shown in Eq. (10) and the Feynman diagrams
in Fig. 3, the effective interaction for the $\Delta L=2$ process can be
written as:

${\cal H}_{\rm eff}=C^{\Delta
L=2}_{1}\left(\overline{\mu}\gamma_{\mu}P_{L}e\right)\left(\overline{\mu}\gamma^{\mu}P_{L}e\right)~{},$
(42)

where the Wilson coefficient $C^{\Delta L=2}_{1}$ is given by:

$C^{\Delta L=2}_{1}=\sum^{2}_{\alpha=1}\frac{(y_{\alpha\mu}y^{*}_{\alpha
e})^{2}}{(4\pi)^{2}m^{2}_{\eta_{\alpha}}}J^{\rm
box}_{1}\left(\frac{m^{2}_{N}}{m^{2}_{\eta_{\alpha}}}\right)-2\frac{y_{1\mu}y^{*}_{1e}y_{2\mu}y^{*}_{2e}}{(4\pi)^{2}m^{2}_{\eta_{2}}}\frac{m^{2}_{N}}{m^{2}_{\eta_{2}}}J^{\rm
box}_{2}\left(\frac{m^{2}_{N}}{m^{2}_{\eta_{2}}},\frac{m^{2}_{\eta_{1}}}{m^{2}_{\eta_{2}}}\right)~{}.$
(43)

Using the transition matrix element
$\langle\overline{M}_{\mu}|O_{1}|M_{\mu}\rangle=f^{2}_{M_{\mu}}m^{2}_{M_{\mu}}$,
the $x$-parameter for para-muonium and ortho-muonium are:

$\displaystyle x_{P}$ $\displaystyle=\frac{4(m_{\rm red}\alpha_{\rm
em})^{3}}{\pi\Gamma_{\mu}}C^{\Delta L=2}_{1}~{},$ $\displaystyle x_{V}$
$\displaystyle=-\frac{12(m_{\rm red}\alpha_{\rm
em})^{3}}{\pi\Gamma_{\mu}}C^{\Delta L=2}_{1}~{},$ (44)

where the reduced mass $m_{\rm red}=m_{\mu}m_{e}/(m_{\mu}+m_{e})$, and
$f^{2}_{M_{\mu}}=4(m_{\rm red}\alpha_{\rm em})^{3}/m_{M_{\mu}}$ Conlin:2020veq
. Hence, Eq. (41) can be expressed as:

$P(M_{\mu}\to\overline{M}_{\mu})^{\rm
exp}=\frac{x^{2}_{P}}{2}+\frac{x^{2}_{V}}{6}~{}.$ (45)

###  III.5 Loop-induced kinetic mixing between $U(1)_{Y}$ and $U(1)_{X}$

Figure 4: One-loop diagrams inducing kinetic mixing between $U(1)_{Y}$ and
$U(1)_{X}$ gauge fields.

In our model the SM particles do not interact with $Z^{\prime}$ boson at tree
level since they are not charged under $U(1)_{X}$. However they interact with
$Z^{\prime}$ through kinetic mixing between $U(1)_{Y}$ and $U(1)_{X}$ gauge
fields induced by the one-loop diagrams shown in Fig. 4. We calculate these
diagrams and obtain kinetic mixing term

$\displaystyle\mathcal{L}_{KM}$
$\displaystyle=-\frac{\sin\epsilon}{2}\hat{B}_{\mu\nu}\hat{Z^{\prime}}^{\mu\nu}~{},$
(46) $\displaystyle\sin\epsilon$
$\displaystyle=\frac{Q_{X}g_{X}g^{\prime}}{2(4\pi^{2})}\int_{0}^{1}dx(1-x)(2x-1)^{2}$
$\displaystyle\times\left(\ln\left[\frac{x(1-x)q^{2}-m_{\eta^{\pm}_{2}}^{2}}{x(1-x)q^{2}-m_{\eta^{\pm}_{1}}^{2}}\right]+\ln\left[\frac{x(1-x)q^{2}-m_{S_{2}}^{2}}{x(1-x)q^{2}-m_{S_{1}}^{2}}\right]\right)~{},$
(47)

where $q^{2}$ corresponds to the momentum carried by the gauge bosons. Note
that the diagrams sum up to give a finite result for the kinetic mixing with
divergences canceled between the contributions from the two inert doublet
scalars carrying opposite charges under $U(1)_{X}$. For $q^{2}\ll
m^{2}_{\eta_{1,2}^{\pm},S_{1,2}}$, we can approximate the formula of
$\epsilon$ by

$\sin\epsilon\sim\frac{Q_{X}g_{X}g^{\prime}}{3(4\pi)^{2}}\left(\ln\frac{m_{\eta^{\pm}_{2}}}{m_{\eta^{\pm}_{1}}}+\ln\frac{m_{S_{1}}}{m_{S_{2}}}\right).$
(48)

We can diagonalize the kinetic terms for $\hat{Z^{\prime}}_{\mu}$ and
$B_{\mu}$ via the following transformations Babu:1997st :

$\displaystyle\begin{split}&\hat{B}_{\mu}=B_{\mu}-\tan{\epsilon}Z^{\prime}_{\mu}~{},\\\
&\hat{Z^{\prime}_{\mu}}=\frac{1}{\cos{\epsilon}}Z^{\prime}_{\mu}~{}.\end{split}$
(49)

where we approximate $\sin\epsilon\simeq\epsilon$ and $\cos\epsilon\simeq 1$
hereafter. After electroweak symmetry breaking, we obtain the mass terms for
neutral gauge bosons as

$\mathcal{L}_{M}=\frac{1}{2}m_{Z}^{2}\tilde{Z}_{\mu}\tilde{Z}^{\mu}+m_{Z}^{2}\epsilon
s_{W}Z^{\prime}_{\mu}\tilde{Z}^{\mu}+\frac{1}{2}m_{Z^{\prime}}^{2}Z^{\prime}_{\mu}Z^{\prime\mu}~{},$
(50)

where $m_{Z}=v\sqrt{g^{2}+g^{\prime 2}}/2$ and
$\tilde{Z}_{\mu}=c_{W}W^{3}_{\mu}-s_{W}B_{\mu}$. We then obtain mass
eigenstates and eigenvalues of the neutral gauge bosons as

$\displaystyle\begin{split}&\\{m_{Z_{1}}^{2},m_{Z_{2}}^{2}\\}=\frac{1}{2}(m_{Z^{\prime}}^{2}+m_{Z}^{2})\mp\frac{1}{2}\sqrt{(m_{Z^{\prime}}^{2}-m_{Z}^{2})^{2}+4\epsilon^{2}s_{W}^{2}m_{Z}^{4}}~{},\\\
&\begin{pmatrix}Z_{1\mu}\\\
Z_{2\mu}\end{pmatrix}=\begin{pmatrix}\cos\theta_{ZZ^{\prime}}&-\sin\theta_{ZZ^{\prime}}\\\
\sin\theta_{ZZ^{\prime}}&\cos\theta_{ZZ^{\prime}}\end{pmatrix}\begin{pmatrix}\tilde{Z}_{\mu}\\\
Z^{\prime}_{\mu}\end{pmatrix}~{},\quad\tan
2\theta_{ZZ^{\prime}}=\frac{2s_{W}\epsilon
m_{Z}^{2}}{m_{Z}^{2}-m_{Z^{\prime}}^{2}}~{},\end{split}$ (51)

where we can approximate $m_{Z_{1}}\simeq m_{Z}$ and $m_{Z_{2}}\simeq
m_{Z^{\prime}}$ for tiny $\epsilon$. The $Z_{1}$ and $Z_{2}$ bosons are to be
identified as the physical massive gauge bosons and, to avoid the pedantry
while being generally not confusing, will be referred to as the $Z$ and
$Z^{\prime}$ bosons. The $Z^{\prime}$ interaction with SM fermions $f$ is now
given by

$\displaystyle\mathcal{L}_{Z^{\prime}\overline{f}f}$
$\displaystyle=\sum_{\chi=L,R}\frac{g}{\cos\theta_{W}}Z^{\prime}_{\mu}\overline{f}_{\chi}\gamma^{\mu}\left[S_{ZZ^{\prime}}(T_{3}-Q\sin^{2}\theta_{W})+C_{ZZ^{\prime}}\epsilon
Y\sin\theta_{W}\right]f_{\chi}$ $\displaystyle\equiv
Z^{\prime}_{\mu}\overline{f}\gamma^{\mu}(C_{L}^{f}P_{L}+C_{R}^{f}P_{R})f~{},$
(52)

where $T_{3}$ is diagonal generator of $SU(2)_{L}$, $Q$ is the electric charge
operator, $S_{ZZ^{\prime}}\equiv\sin\theta_{ZZ^{\prime}}$, and
$C_{ZZ^{\prime}}\equiv\cos\theta_{ZZ^{\prime}}$.

##  IV Numerical Analysis

In this section, we perform numerical scans for the allowed parameter space
and the corresponding ranges of various observables, which are then compared
with current experimental bounds or the sensitivities that ongoing/future
experiments can probe. We also calculate the DM relic density and check
against the constraint of the direct search limit.

###  IV.1 Inputs and constraints

From Eq. (20), it can be seen that when the neutrino mixing angles and masses
are determined, the parameters of ${\bf y}_{1,2}$, $m_{S_{i},A_{i}}$, and
$m_{N}$ can be bounded. In order to get the allowed ranges for the free
parameters, we take the values of the neutrino oscillation parameters obtained
from the global fit in Refs. Esteban:2020cvm ; Gonzalez-Garcia:2021dve  as the
inputs, where the global fit results with $3\sigma$ errors are given in Table
2. Based upon the fit results, the ranges of the Majorana mass matrix elements
in units of eV for the normal ordering (NO) and inverted ordering (IO) can be
respectively estimated as:

$\displaystyle\begin{split}\begin{pmatrix}|m^{\nu}_{11}|&|m^{\nu}_{12}|&|m^{\nu}_{13}|\\\
|m^{\nu}_{21}|&|m^{\nu}_{22}|&|m^{\nu}_{23}|\\\
|m^{\nu}_{31}|&|m^{\nu}_{32}|&|m^{\nu}_{33}|\end{pmatrix}_{\rm
NO}&\simeq\begin{pmatrix}0.09-0.42&0.095-0.909&0.087-0.906\\\
0.095-0.909&1.51-3.31&2.03-2.81\\\
0.087-0.906&2.03-2.81&1.46-3.27\end{pmatrix}\times 10^{-2}~{},\\\
\begin{pmatrix}|m^{\nu}_{11}|&|m^{\nu}_{12}|&|m^{\nu}_{13}|\\\
|m^{\nu}_{21}|&|m^{\nu}_{22}|&|m^{\nu}_{23}|\\\
|m^{\nu}_{31}|&|m^{\nu}_{32}|&|m^{\nu}_{33}|\end{pmatrix}_{\rm
IO}&\simeq\begin{pmatrix}1.50-4.97&0.-3.81&0.-3.88\\\
0.-3.81&0.-3.05&0.65-2.60\\\ 0.-3.88&0.65-2.60&0.-3.16\end{pmatrix}\times
10^{-2}~{},\end{split}$ (53)

where $m_{1(3)}=0$ for NO (IO) is applied. Since we do not have any
information on the Majorana phases at the moment, we take
$\alpha_{1,2}\in[-\pi,\pi]$. We then use Eq. (53) to bound the free
parameters.

Table 2: Global fit with $3\sigma$ ranges based upon the observations of neutrino oscillations Esteban:2020cvm ; Gonzalez-Garcia:2021dve . | $\theta_{12}/^{\circ}$ | $\theta_{13}/^{\circ}$ | $\theta_{23}/^{\circ}$ | $\delta_{\rm CP}/^{\circ}$ |  $\Delta m^{2}_{21}\times 10^{5}$ [eV2] |  $\Delta m^{2}_{3\ell}\times 10^{3}$ [eV2]
---|---|---|---|---|---|---
NO |  $(31.27,35.87)$ |  $(8.25,8.98)$ |  $(39.7,50.9)$ |  $(144,350)$ |  $(6.82,8.04)$ |  $(2.43,2.59)$
IO | $(31.27,35.87)$ | $(8.24,9.02)$ | $(39.8,51.6)$ | $(194,345)$ | $(6.82,8.04)$ | $(-2.574,-2.410)$

The parameter space can be further constrained by various LFV processes, whose
current upper bounds are given in Table 3. With $S_{B}(B_{0})=0.75$, the upper
limit on the probability of the muonium oscillation is taken as
$P(M_{\mu}\to\overline{M}_{\mu})^{\rm exp}<11.1\times 10^{-11}$
Willmann:1998gd ; Conlin:2020veq .

Table 3: Current upper bounds on various LFV processes. Except the $\mu-e$ conversion rate in titanium quoted from Ref. SINDRUMII:1993gxf , all other values are obtained from the Particle Data Group PDG . | $\mu\to e\gamma$ | $\mu\to 3e$ | $\tau\to e\gamma$ | $\tau\to\mu\gamma$ | $\tau\to 3e$ | $\tau\to 3\mu$ | $CR(\mu-e,{\rm Ti})$
---|---|---|---|---|---|---|---
U.L. | $4.2\times 10^{-13}$ | $1.0\times 10^{-12}$ | $3.3\times 10^{-8}$ | $4.4\times 10^{-8}$ | $2.7\times 10^{-8}$ | $2.1\times 10^{-8}$ |  $4.3\times 10^{-12}$ SINDRUMII:1993gxf

There are $12$ free parameters involved in the model, and we only have $6$
observables from the neutrino oscillation experiments. In order to make the
parameter scans more efficient, following the strict constraint from the
$\mu\to e\gamma$ decay shown in Eqs. (22) and (24), we assume

$y_{2\mu}=-\frac{m^{2}_{S_{2}}}{m^{2}_{S_{1}}}\frac{y^{*}_{1e}}{y^{*}_{2e}}\frac{J^{\gamma}_{2}\left(\frac{m^{2}_{N}}{m^{2}_{S_{1}}}\right)}{J^{\gamma}_{2}\left(\frac{m^{2}_{N}}{m^{2}_{S_{2}}}\right)}+\xi
e^{i\phi}~{},$ (54)

where $\xi$ is a real parameter and $\phi$ is its phase. The parameter ranges
used to scan the parameter space are taken as follows:

$\displaystyle\begin{split}&m_{N}\in[100,1000]~{}{\rm
GeV}~{},~{}m_{S_{1,2}}\in[800,2500]~{}{\rm GeV}~{},\\\
&y_{1i,2i\neq\mu}\in[-1,1]~{},~{}\phi_{1,3}\in[-\pi,\pi]~{},~{}\xi\in[-0.1,0.1]~{},~{}\phi\in[-\pi,\pi]~{}.\end{split}$
(55)

In addition, we also assume $m_{N}<m_{S_{1}}<m_{S_{2}}$ and
$m_{\eta^{\pm}_{i}}=m_{S_{i}}$ in the numerical estimates.

Under the constraints of Table 2 and Table 3, we have sampled $10^{9}$ points
with the ranges of involved parameters shown in Fig. 5, where the black (red)
points denote the NO (IO) neutrino mass. Figs. 5(a)–(c) show the allowed
ranges of various products of Yukawa couplings in absolute values that will
appear in the calculations of LFV processes, while Fig. 5(d) shows those of
$m_{N,S_{1}}$. The plots indicate that the constrained Yukawa couplings in the
model are not sensitive to the neutrino mass ordering.

Figure 5:  (a)-(c) Constraints on various products of Yukawa couplings, and
(d) allowed ranges of $m_{N,S_{1}}$. The black and red points denote the
results of NO and IO, respectively.

###  IV.2 Correlations among the LFV processes

In this subsection, we numerically discuss the influence of parameters on the
rare LFV processes and the correlations among the considered LFV processes in
detail. Since the $\mu\to e\gamma$ decay is the most constraining LFV process
at present and in foreseeable future, we first show its branching ratio
dependence on the parameters and then investigate various correlations among
the LFV processes.

According to Eq. (24), the branching ratio for $\mu\to e\gamma$ can be taken
as a function of $|y_{ie}y^{*}_{i\mu}|/m^{2}_{\eta^{\pm}_{i}}$. Based upon the
scanning results exhibited in Fig. 5, we show the scatter plots for $BR(\mu\to
e\gamma)$ as a function of $|y_{1e}y^{*}_{1\mu}|/m^{2}_{\eta^{\pm}_{i}}$ and
$|y_{2e}y^{*}_{2\mu}|/m^{2}_{\eta^{\pm}_{i}}$ in Fig. 6(a) and (b),
respectively, where the horizontal dashed line denotes the sensitivity of MEG
II experiment MEGII:2018kmf . Since the MEG II experiment can only probe
$\mu\to e\gamma$ down to the level of $6\times 10^{-14}$, we have further
restricted the model parameters to satisfy $BR(\mu\to
e\gamma)\in(0.05,5)\times 10^{-13}$ in the plots. As shown in the plots, the
difference between the dependence on
$|y_{1e}y^{*}_{1\mu}|/m^{2}_{\eta^{\pm}_{1}}$ and that on
$|y_{2e}y^{*}_{2\mu}|/m^{2}_{\eta^{\pm}_{2}}$ is insignificant. As given in
Eqs. (32) and (34), the $\mu-e$ conversion rate arises from the photon-penguin
diagram in the model. Thus, we show $CR(\mu-e,{\rm Ti})$ as a function of
$|y_{1e}y^{*}_{1\mu}|/m^{2}_{\eta^{\pm}_{1}}$ in Fig. 6(c), where the dashed
line is the sensitivity of COMET COMET:2018auw  and Mu2e Diociaiuti:2020yvo
experiments. As such, these experiments have the capability to probe most of
the considered parameter space through the $\mu-e$ conversion rate. The
correlation between $BR(\mu\to e\gamma)$ and $CR(\mu\to e,{\rm Ti})$ is
plotted in Fig. 6(d). Again, when the measurement of $CR(\mu-e,{\rm Ti})\sim
10^{-18}$ achieves the expected sensitivity, the parameter space with
$BR(\mu\to e\gamma)\gtrsim O(10^{-13})$ is mostly covered. Hence, the $\mu-e$
conversion process has the potential to become the most stringent constraint
among the LFV processes.

Figure 6:  (a) and (b) $BR(\mu\to e\gamma)$ (in units of $10^{-13}$) as a
function of $|y_{ie}y^{*}_{i\mu}|/m^{2}_{\eta^{\pm}_{i}}$; (c) $\mu-e$
conversion rate as a function of of
$|y_{1e}y^{*}_{1\mu}|/m^{2}_{\eta^{\pm}_{1}}$; and (d) correlation between
$BR(\mu\to e\gamma)$ and $CR(\mu-e,{\rm Ti})$, where the dashed lines are the
sensitivities of currently ongoing experiments. The black and red points
denote the results of NO and IO, respectively.

As stated earlier, in addition to the photon-penguin diagrams, the $\mu\to 3e$
decay can be induced by the box diagrams in the model. To exhibit the role of
the box diagrams, the ratio of $BR(\mu\to 3e)$ purely from the photon-penguin
contribution to $BR(\mu\to e\gamma)$ can be expressed as:

$R_{3e/e\gamma}=\frac{BR(\mu\to 3e)_{\gamma}}{BR(\mu\to
e\gamma)}\approx\frac{\alpha_{\rm
em}}{8\pi}\left[\frac{|C^{\gamma}_{1e\mu}|^{2}}{|C^{\gamma}_{2e\mu}|^{2}}+\frac{16}{3}\ln\frac{m_{\mu}}{m_{e}}-\frac{22}{3}-4\operatorname{Re}\left(\frac{C^{\gamma}_{1e\mu}}{C^{\gamma}_{2e\mu}}\right)\right]~{}.$
(56)

With $C^{\gamma}_{1e\mu}/C^{\gamma}_{2e\mu}\sim{\cal O}(1)$, the ratio in Eq.
(56) can be estimated to be $R_{3e/e\gamma}\sim 5\times 10^{-3}$. Using the
constrained parameter values shown in Fig. 5, indeed, we approximately obtain
$R_{3e/e\gamma}\in(0.5,5)\times 10^{-2}$. We can conclude that if $BR(\mu\to
3e)>BR(\mu\to e\gamma)$ is observed in experiments, the enhancement should
arise from other effects than from the photon-penguin diagrams.

If we assume that the photon-penguin contribution to $\mu\to 3e$ is
subleading, it is interesting to compare the $\mu\to 3e$ decay with the
munoium oscillation, where both processes are dominated by the box diagrams.
Implementing the constrained parameter values shown in Fig. 5 to the formulas
in Eqs. (31) and (45), we plot the correlation between
$P(M_{\mu}-\overline{M}_{\mu})$ and $BR(\mu\to 3e)$ in Fig. 7(a), where the
dashed lines label the sensitivities of Mu3e and MACE experiments. It is seen
that the predicted values of $P(M_{\mu}-\overline{M}_{\mu})$ in the model are
mostly below than the current experimental limit of $4.3\times 10^{-12}$ when
the $\mu\to 3e$ decay is constrained to satisfy the current upper limit of
$1.0\times 10^{-12}$. However, when $BR(\mu\to 3e)$ reaches the sensitivity of
$10^{-16}$, the decay can cover most of the considered parameter space, i.e.,
under the presumption that $BR(\mu\to e\gamma)\in(0.05,5)10^{-13}$. We thus
see that the two strictest constraints on the $\mu\to e$ transitions come from
the $\mu-e$ conversion in nucleus and the $\mu\to 3e$ decay. Their correlation
in the model can be found in Fig. 7(b). As such, if we do not see any evidence
of the model in the LFV experiments, the highly sensitive measurements of
$BR(\mu\to 3e)$ and $\mu-e$ conversion rate severely constrain the $\mu\to
e\gamma$ decay and munoium oscillation processes.

Figure 7: Correlation between $BR(\mu\to 3e)$ and (a)
$P(M_{\mu}-\overline{M}_{\mu})$ and (b) $CR(\mu-e,{\rm Ti})$, where the dashed
lines label the experimental sensitivities. The black and red points denote
the results of NO and IO, respectively.

After discussing the rare $\mu\to e$ transition, we discuss in the following
analysis the lepton flavor-changing effects on the heavy $\tau$-lepton decays,
where the sensitivities assuming 50 ab-1 of data at Belle II can reach Belle-
II:2018jsg :

$\displaystyle\begin{split}BR(\tau\to\ell\gamma)&\sim 10^{-8}-10^{-9}~{},\\\
BR(\tau\to 3\ell)&\sim 10^{-9}-10^{-10}~{},\end{split}$ (57)

with $\ell=e,\mu$. Analogous to $\mu\to e\gamma$, the radiative $\tau$ decay
also arises from the same types of diagrams. Therefore, using Eq. (24) and the
constrained parameter values, we plot the correlation between $BR(\tau\to
e\gamma)$ and $BR(\mu\to e\gamma)$ in Fig. 8(a), where the dashed lines are
the sensitivities of MEG II and Belle II experiments. Since the resulting
branching ratio for $\tau\to e\gamma$ is lower than the sensitivity of Belle
II, it is difficult to observe the $\tau\to e\gamma$ decay in the model. In
Fig. 8(b), we show the correlation between $BR(\tau\to\mu\gamma)$ and
$BR(\tau\to e\gamma)$. It can be found that unlike the $\tau\to e\gamma$
decay, the branching ratio for $\tau\to\mu\gamma$ can be as large as
$O(10^{-8})$. As such, $\tau\to\mu\gamma$ serves as a good candidate to probe
the new physics effects of the model in Belle II experiment.

Figure 8:  Correlation of BR in LFV processes: (a) $\tau\to e\gamma$ and
$\mu\to e\gamma$, and (b) $\tau\to\mu\gamma$ and $\tau\to e\gamma$. The black
and red points denote the results of NO and IO, respectively.

As discussed earlier, although the $\ell_{i}\to 3\ell_{j}$ decay is induced by
the photon-penguin and box diagrams, the effects of the latter are more
dominant. In order to reveal their contributions to $\tau\to 3e$ and $\tau\to
3\mu$, we show their branching ratio correlations with $\mu\to 3e$ in Fig. 9.
Similar to $\tau\to e\gamma$, it is difficult for Belle II to reach the
predicted $BR(\tau\to 3e)$ in the model. Nevertheless, the $\tau\to 3\mu$
decay is more promising to detect at Belle II because the value can reach
$O(10^{-9}-10^{-10})$.

Figure 9:  Correlation of BR between $\mu\to 3e$ and (a) $\tau\to 3e$ and (b)
$\tau\to 3\mu$. The black and red points denote the results of NO and IO,
respectively.

We finally make some remarks on the lepton $(g-2)$’s shown in Eq. (25). It is
known that the lepton $(g-2)$ mediated by the charged scalar boson is usually
negative. Although the sign of $a_{\mu}$ in the model contradicts with that
given by the recent E989 experiment at Fermilab Abi:2021gix , its value
$|a_{\mu}|<10^{-11}$ due to the $m^{2}_{\mu}$ dependence and
$m_{\eta^{\pm}_{i}}\sim O(1)$ TeV. Therefore, if the muon $g-2$ anomaly can be
resolved within the SM Borsanyi:2020mff , the negative $a_{\mu}$ of the model
will not cause a serious problem. On the other hand, the electron $(g-2)$ in
the model can be estimated as $|a_{e}|<10^{-17}$ and is negligible.

###  IV.3 DM relic density and DM-nucleon scattering cross section

In this subsection we discuss the DM phenomenology in our model. Our DM
candidate is the vector-like neutral fermion $N$ which interacts with the SM
particles via the Yukawa couplings in Eq. (9) and the $Z^{\prime}$ exchange
with the kinetic mixing effect. In the estimate of relic density, the dominant
DM annihilation processes are summarized as follows:

  * •

$N\overline{N}\to Z^{\prime}\to f_{SM}\overline{f}_{SM},W^{+}W^{-}$.

  * •

$N\overline{N}\to Z^{\prime}Z^{\prime}$.

  * •

$\\{N\overline{N},NN,\overline{N}\overline{N}\\}\to\\{\nu\overline{\nu},\nu\nu,\overline{\nu}\overline{\nu},\ell^{+}\ell^{-}\\}$,

where $\ell^{\pm}$ and $f_{SM}$ denote the SM charged leptons and fermions,
respectively. The first two processes are mediated by the new gauge
interactions, while the last channels rely mostly on the new Yukawa couplings.
We estimate the relic density of DM using micrOMEGAs 5.2.4 Belanger:2014vza
implemented with the new interaction vertices in the model.

We first discuss the relic density by considering only $Z^{\prime}$
interactions and choosing vanishing Yukawa couplings for the DM mass in range
of 100 GeV to 1 TeV. For illustration purposes, we consider two specific
cases; (1) $m_{N}=2m_{Z^{\prime}}$ where the $N\overline{N}\to
Z^{\prime}Z^{\prime}$ process is dominant, (2) $m_{Z^{\prime}}=2.025m_{N}$
where the $N\overline{N}\to
Z^{\prime}\to\\{f_{SM}\overline{f}_{SM},W^{+}W^{-}\\}$ processes are dominant.
The relic density is then estimated by scanning the values of
$\\{m_{N},g_{X}\\}$. The left and right plots in Fig. 10 show the parameter
regions in the $m_{N}$-$g_{X}$ plane that give the relic density $0.11<\Omega
h^{2}<0.13$ PDG  for cases (1) and (2), respectively. For case (1), the gauge
coupling around $0.2\lesssim g_{X}\lesssim 0.7$ can accommodate the observed
relic density in the assumed DM mass range. For case (2), on the other hand, a
larger gauge coupling is required and the observed relic density can be
explained only when $m_{N}\lesssim 400$ GeV imposing perturbative condition
$g_{X}<\sqrt{4\pi}$. This behavior is due to that fact that the small kinetic
mixing for the $s$-channel annihilation via $Z^{\prime}$ exchange is
suppressed.

Figure 10: Scatter plots of the parameter region in the $\\{m_{N}$-$g_{X}\\}$
plane that satisfies $0.11<\Omega h^{2}<0.13$, considering only the
$Z^{\prime}$ interactions. The Left and the right plots correspond to cases of
$m_{N}=2m_{Z^{\prime}}$ and $m_{Z^{\prime}}=2.025m_{N}$, respectively.

Secondly, we discuss the relic density by use of the Yukawa couplings that
satisfy the neutrino oscillation observations and flavor-changing constraints
given in the previous subsections. For illustration purposes, we consider a
vanishing gauge coupling $g_{X}$ and focus on the effects of Yukawa couplings.
Fig. 11 shows the DM relic density as a function of $m_{N}$ where the black
and red points correspond to allowed parameter sets in the NO and IO
scenarios, respectively. We find that the observed relic density can be
explained by the Yukawa interactions for $m_{N}\lesssim 650$ GeV and that it
gets smaller than the observed value in the heavier DM region. This is because
the Yukawa couplings get larger in the heavier mass region, as required to fit
the neutrino oscillation data. We note in passing that it is possible to
utilize the $Z^{\prime}$ interactions to shift the scatter points downwards by
turning on the gauge coupling $g_{X}$ and tuning its value and the
$Z^{\prime}$ mass.

Figure 11: Scatter plot of the relic density of $N$ for the allowed parameter
sets with $g_{X}\to 0$. The black and red points denote the results of NO and
IO, respectively.

Finally we discuss the DM-nucleon scattering cross section via the
$Z^{\prime}$ exchange. Here we estimate the cross section in the non-
relativistic limit and obtain

$\displaystyle\begin{split}\sigma_{n(p)}\simeq&~{}\frac{g_{X}^{2}}{2\pi}\frac{\mu_{n(p)}^{2}}{m^{4}_{Z^{\prime}}}\left[\left(C_{L}^{n(p)}\right)^{2}+\left(C_{R}^{n(p)}\right)^{2}+C_{L}^{n(p)}C_{R}^{n(p)}\right]~{},\\\
\mbox{with }&C_{L(R)}^{n}=C_{L(R)}^{u}+2C_{L(R)}^{d}~{},\quad
C_{L(R)}^{p}=2C_{L(R)}^{u}+C_{L(R)}^{d}~{},\end{split}$ (58)

where $n~{}(p)$ stands for the neutron (proton), $\mu_{n(p)}\equiv
m_{N}m_{n(p)}/(m_{N}+m_{n(p)})$ and $C_{L(R)}^{f}$ is given in Eq. (III.5). We
show in Fig. 12 the DM-neutron scattering cross section as a function of
$g_{X}$ for $m_{Z^{\prime}}=\\{50,100,250\\}$ GeV and $m_{N}=500$ GeV. The
gray region has $\Omega h^{2}\lesssim 0.11$ and the cyan region is excluded by
XENON1T XENON:2018voc . In making this plot, we assume the dominance
$Z^{\prime}$ interactions and ignore the Yukawa couplings. The DM-proton
scattering cross section is very similar and thus omitted here.

Figure 12: DM-neutron scattering cross section as a function of $g_{X}$ for
$m_{Z^{\prime}}=\\{50,100,250\\}$ GeV where ray region indicates $g_{X}$
region providing $\Omega h^{2}\lesssim 0.11$ and light blue region indicate
excluded region by XENON1T, for $m_{N}=500$ GeV.

##  V Summary

Based upon the concept of radiative neutrino mass in scotogenic model proposed
in Ma:2006km , we have studied the possibility for the Majorana neutrino mass
to arise from an unbroken Stueckelberg $U(1)_{X}$ gauge model. It is found
that although we need two inert Higgs doublets to generate the neutrino mass,
just one extra vector-like singlet fermion is sufficient to fit the neutrino
data while in contrast at least two right-handed singlet fermions are required
in the Ma model. Using the proposed model, we have studied its implications on
the low-energy lepton flavor violating (LFV) processes.

We have found that in the model, the resulting $BR(\mu\to e\gamma)$ can fit
the currently experimental upper limit. However, the planned sensitivities of
the experiments on the $\mu-e$ conversion and the $\mu\to 3e$ decay in
COMET/Mu2e and Mu3e can cover most of the parameter space that can lead to
$BR(\mu\to e\gamma)\sim 6\times 10^{-14}$, the designed sensitivity of MEG II
experiment. Hence, the $\mu-e$ conversion in nucleus and $\mu\to 3e$ will be
the most promising processes to detect the new physics in the $\mu\to e$
transitions.

In heavy lepton decays, although the predicted $BR(\tau\to e\gamma)$ and
$BR(\tau\to 3e)$ are lower than $O(10^{-9})$ and $O(10^{-10})$, the branching
ratios for $\tau\to\mu\gamma$ and $\tau\to 3\mu$ in the model can cover the
range from $O(10^{-8})$ to $O(10^{-9})$ and $O(10^{-10})$, respectively, and
can be probed by Belle II experiment.

We have discussed the dark matter (DM) phenomenology, including estimates of
the relic density and the DM-nucleon scattering cross section for the DM mass
in the range of 100 GeV to 1 TeV. The relic density has been estimated for two
cases: (i) the $Z^{\prime}$ interaction is dominant (ii) the Yukawa
interaction is dominant. In case (i), the observed relic density can be
explained with $0.2\lesssim g_{X}\lesssim 0.7$ when DM annihilates into
$Z^{\prime}Z^{\prime}$ while a larger gauge coupling is required when DM
annihilate into the SM particles via the $s$-channel $Z^{\prime}$ exchange. In
case (ii), we make use of the parameter sets obtained from neutrino
oscillation observations and flavor-changing constraints. We find that the
relic density can be explained for a DM mass less than around 650 GeV. The DM-
nucleon cross section is obtained by considering the $Z^{\prime}$ exchange
process and evaluated by fixing some parameters. We have found that it is
possible to avoid constraints from direct detection while satisfying the relic
density by properly choosing the gauge coupling and $Z^{\prime}$ mass.

## Acknowledgments

This work was supported in part by the Ministry of Science and Technology,
Taiwan under the Grant Nos. MOST-110-2112-M-006-010-MY2 and
MOST-108-2112-M-002-005-MY3. The work is also supported by the Fundamental
Research Funds for the Central Universities (T. N.).

## References

  * (1) E. Ma, Phys. Rev. D 73, 077301 (2006) [hep-ph/0601225].
  * (2) N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)] [arXiv:1807.06209 [astro-ph.CO]].
  * (3) A. M. Baldini et al. [MEG II], Eur. Phys. J. C 78, no.5, 380 (2018) [arXiv:1801.04688 [physics.ins-det]].
  * (4) C. M. Perez and L. Vigani, Universe 7, no.11, 420 (2021)
  * (5) R. Abramishvili et al. [COMET], PTEP 2020, no.3, 033C01 (2020) [arXiv:1812.09018 [physics.ins-det]].
  * (6) E. Diociaiuti, PoS EPS-HEP2019, 232 (2020)
  * (7) R. J. Barlow, Nucl. Phys. B Proc. Suppl. 218, 44-49 (2011)
  * (8) J. Tang et al., Letter of interest contribution to snowmass21, https://www.snowmass21.org/docs/files/summaries/RF/SNOWMASS21-RF5_RF0_Jian_ang-126.pdf.
  * (9) J. Leite, A. Morales, J. W. F. Valle and C. A. Vaquera-Araujo, Phys. Lett. B 807, 135537 (2020) [arXiv:2003.02950 [hep-ph]].
  * (10) E. Aprile et al. [XENON], Phys. Rev. Lett. 121 (2018) no.11, 111302 [arXiv:1805.12562 [astro-ph.CO]].
  * (11) E. Kou et al. [Belle-II], PTEP 2019, no.12, 123C01 (2019) [erratum: PTEP 2020, no.2, 029201 (2020)] [arXiv:1808.10567 [hep-ex]].
  * (12) T. Toma and A. Vicente, JHEP 01, 160 (2014) [arXiv:1312.2840 [hep-ph]].
  * (13) S. Weinberg, Phys. Rev. Lett. 43, 1566-1570 (1979).
  * (14) J. Hisano, T. Moroi, K. Tobe and M. Yamaguchi, Phys. Rev. D 53, 2442-2459 (1996) [arXiv:hep-ph/9510309 [hep-ph]].
  * (15) Y. Kuno and Y. Okada, Rev. Mod. Phys. 73, 151-202 (2001) [arXiv:hep-ph/9909265 [hep-ph]].
  * (16) E. Arganda, M. J. Herrero and A. M. Teixeira, JHEP 10, 104 (2007) [arXiv:0707.2955 [hep-ph]].
  * (17) T. S. Kosmas, S. Kovalenko and I. Schmidt, Phys. Lett. B 511, 203 (2001) [arXiv:hep-ph/0102101 [hep-ph]].
  * (18) R. Conlin and A. A. Petrov, Phys. Rev. D 102, no.9, 095001 (2020) [arXiv:2005.10276 [hep-ph]].
  * (19) L. Willmann, P. V. Schmidt, H. P. Wirtz, R. Abela, V. Baranov, J. Bagaturia, W. H. Bertl, R. Engfer, A. Grossmann and V. W. Hughes, et al. Phys. Rev. Lett. 82, 49-52 (1999) [arXiv:hep-ex/9807011 [hep-ex]].
  * (20) K. S. Babu, C. F. Kolda and J. March-Russell, Phys. Rev. D 57, 6788-6792 (1998) [arXiv:hep-ph/9710441 [hep-ph]].
  * (21) I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz and A. Zhou, JHEP 09, 178 (2020) [arXiv:2007.14792 [hep-ph]].
  * (22) M. C. Gonzalez-Garcia, M. Maltoni and T. Schwetz, Universe 7, no.12, 459 (2021) doi:10.3390/universe7120459 [arXiv:2111.03086 [hep-ph]].
  * (23) C. Dohmen et al. [SINDRUM II], Phys. Lett. B 317, 631-636 (1993)
  * (24) P. A. Zyla et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2020, 083C01 (2020).
  * (25) B. Abi et al. [Muon g-2], Phys. Rev. Lett. 126, 141801 (2021) [arXiv:2104.03281 [hep-ex]].
  * (26) S. Borsanyi et al., Nature (2021), 2002.12347. [arXiv:2002.12347 [hep-lat]].
  * (27) G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, Comput. Phys. Commun. 192, 322 (2015) [arXiv:1407.6129 [hep-ph]].