Lepton flavor violation and scotogenic Majorana neutrino mass in a Stueckelberg $U(1)_X$ model

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. [1] 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 Ω DM h 2 = 0.120 ± 0.001 [2].
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 µ → eγ, µ → 3e, µ − 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 µ → eγ decay in MEG II experiment can reach the sensitivity of 6 × 10 −14 [3]; BR(µ → 3e) 10 −16 is expected in Mu3e experiment [4]; the µ−e conversion rate in COMET at J-PARC [5] and Mu2e at FNAL [6] can reach the level of 10 −17 − 10 −18 , whereas the PRISM/PRIME experiment will push the conversion rate down below 10 −18 [7]; and the probability of muonium-antimuonium conversion in the new generation experiment with the sensitivity of O(10 −14 ) is proposed by the MACE collaboration at CSNC [8].
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 ) and, thus, the DM relic density may not be correctly reproduced through the resonant Z production when we focus on the minimal extension of the SM. To guarantee the coexistence of a massive Z 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 Stueck-elberg mechanism was studied in Ref. [9]. 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 (η 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 † η 1 )(H † η 2 ) in the scalar potential. In addition, since the Yukawa couplings to the SM lepton can have different structures, such as LN Rη1 and L C N Lη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 µ − e conversion rate and the µ → 3e decay will exclude most of the parameter space that lead to BR(µ → eγ) ∼ 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 µ → e LFV processes. In addition, the resulting BR(τ → µγ) and BR(τ → 3µ) can reach the significant levels of O(10 −9 ) and O(10 −10 ), respectively, as expected to probe by Belle II experiment [11].
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 and the gauge field of U (1) Y can be induced via quantum corrections; thus, the Z 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 -mediated s-channel process N N → Z → FF and the N -mediated t-channel one where F (F ) denotes possible SM particles. We will demonstrate that in addition to the Yukawa coupling scenario, the induced Z gauge coupling one can also explain the DM relic density Ω DM h 2 . Moreover, the Z 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 [10].
The structure of this paper is organized as follows. We introduce the model and derive the relevant scalar couplings, Yukawa couplings, and Z 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 η 1,2 , where their representations and charge assignments are given in Table I. Note that here we use the convention that the electromagnetic charge Q = T 3 L + Y /2. To break the lepton number in the scalar sector, η 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 η 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.

A. 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 with m 2 Each of the two 2 × 2 mass-square matrices can be diagonalized by the corresponding or- where the O(θ ξ ) matrix can be parametrized as: 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 θ S = −θ A ≡ θ. The eigenvalues for m 2 S are found to be: For the physical pseudoscalars A 1,2 , we have m 2 A 1(2) = m 2 S 1 (2) . The mixing angle θ is given by: 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 λ 5 term violates the lepton number and eventually leads to the Majorana mass, its value has to be sufficiently small, λ 5 1 [1], as alluded to earlier.
As a result, the off-diagonal mass matrix element |m 2 12 | is suppressed and the mixing angle θ 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 [12] and take λ 5 = 10 −9 .
According to the U (1) X charges listed in Table I, the Yukawa couplings of η 1,2 and N to the SM leptons are given by: where the flavor indices are suppressed, L T = (ν , ) denotes the left-handed lepton doublet in the SM,η j = iτ 2 η * j with τ 2 being a Pauli matrix, N C = Cγ 0 N * is the charge conjugation of N , and m N is the mass of N . Although 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 y 1 is real and y 2 is complex, as parametrized by y 2k = |y 2k |e iφ k (k = 1, 2, 3). Using Eqs. (2) and (6), the Yukawa couplings in Eq. (9) can be decomposed as: where we denote c θ (s θ ) ≡ cos θ (sin θ). We note that if we take the charged lepton L as the mass eigenstate, in general, the η ± 1,2 Yukawa couplings are modified as V L y i , where V L is a unitary matrix introduced for diagonalizing the charged lepton mass matrix. Therefore, the Yukawa couplings of η ± i are generally different from those of S i and A i . Since we know nothing about the information of V L in the SM, we adopt V 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 η ± i .

B. U (1) X Stueckelberg gauge couplings
The Lagrangian invariant under the SU (2) L × U (1) Y × U (1) X gauge symmetry with the Stueckelberg gauge field included is given by: where W i µν ,B µν , andẐ µν 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, β is the gauge fixing parameter, and m Z is the mass of U (1) X Stueckelberg gauge field. We note that the kinetic mixing termB µνẐ µν can be rotated away by redefining the gauge fieldsB µ andẐ µ . However, we will show later that the mixings in γ-Z and Z-Z can be induced through radiative corrections.
The covariant derivatives on N and η i are expressed as: where q N (η i ) denotes the U (1) X charge of N (η i ), and g and g are the gauge couplings of SU (2) L and U (1) Y , respectively. The N N Z interaction can be immediately obtained as: Although the quartic gauge couplings of η i to gauge bosons W ± , Z, γ, andẐ can be derived from the kinetic terms of η i , these couplings are irrelevant to our study here and are not presented explicitly. The various trilinear gauge couplings of η i from (D µ η i ) † (D µ η i ) can be summarized as: The W ± , Z, and γ gauge bosons and the Weinberg angle θ W are defined through the relations: where c W (s W ) ≡ cos θ W (sin θ W ) and e = gs W = g 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: with 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, i → j γ, i → 3 j , µ − 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. To obtain the Majorana neutrino mass in the model, we need both Yukawa couplings y 1 and y 2 , which can ensure that the structure of Majorana mass term L T CL can be achieved.
In addition, since η 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 to generate the Weinberg's dimension-5 operator LHHL [13]. 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: where we have included A i contributions, m A i = m S i is used, and the symmetric Yukawa couplings y ij in flavor indices are defined as: In the Ma model, the involved Yukawa couplings for m ν ij appear in the combination of y i y j ; thus, the mass matrix elements have strong correlations. It needs at least two singlet righthanded fermions to explain the neutrino data. In our model, due to the introduction of one more inert Higgs doublet η 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 ν can be diagonalized through is the diagonal mass matrix. The unitary matrix U can be parametrized as: B.
Radiative lepton decay i → j γ Among the LFV processes, the most constraining process is the radiative muon decay µ → eγ. 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 τ → (e, µ)γ, in the following, we calculate the branching ratios for the i → j γ decays, where j is the lighter lepton and its mass is neglected unless otherwise stated.
The radiative LFV process in the model arises from the η ± i -mediated diagrams, as shown in Fig. 2. Using the Yukawa couplings in Eq. (10), the loop-induced effective Lagrangian for i → j γ ( * ) can be expressed as: where the loop induced Wilson coefficients are obtained as: and the loop integral functions are defined by The emitted photon can be on-shell (k 2 = 0) or off-shell (k 2 = 0), where the latter can be used to study the i → 3 j process, in which a + j − 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 i → j γ can be obtained as [12]: with α em = e 2 /(4π).
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 j by i . As a result, the lepton (g − 2) can be directly obtained as: C. i → 3 j from penguin and box diagrams We now discuss the three-body lepton flavor-changing decays. For simplicity, we only concentrate on the i → 3 j decay. In the model, the LFV processes i → 3 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 i → j γ * has been given in Eq. (21), the transition amplitude with the assigned momenta for the i (p) → j (p 1 ) + (p 2 ) − (p 3 ) decay can be written as: where the photon coupling to the SM charged lepton shown in Eq. (16) is applied.
The Feynman diagrams for i → j Z * are similar to the photon diagrams shown in Fig. 2, where we use Z * instead of γ * . Using the Z gauge couplings to η ± and the SM lepton, given in Eqs. (14) and (16), the one-loop induced effective interaction for i j Z can be obtained as: where we have dropped the small contributions from k 2 and k terms [12]. It can be clearly seen that the loop-induced Z coupling is proportional to the product of the lepton masses, i.e., m i m 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 η ± i , appears in the model. When a necessary mass insertion occurs in the initial lepton 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 j . Hence, we get the result proportional to m i m j . Due to the fact that m j m i and m i m j /m 2 The box diagrams contributing to i → 3 j are partly shown in Fig. 3. In addition to the diagrams mediated by the same η ± 1(2) , the other possible diagrams are mediated by both η ± 1 and η ± 2 and involve chirality-flipping factor m 2 N . Using the Yukawa couplings in Eq. (10), the four-fermion interaction for i → 3 j is obtained as: where C box 1 denotes the contributions from diagrams involving the same η ± 1(2) and C box 2 are from those involving both η ± 1 and η ± 2 , and the effective Wilson coefficients are given by: with the loop integral functions given by Although C box 2ji seems to have a chiral enhancement, its contribution is in fact somewhat smaller than C box 1ji due to the assumed condition that m N < m η 2 . Combining the contributions from the photon-penguin and box diagrams, the branching ratio for the i → 3 j decay can be written as [12]: where C box ji = C box 1ji + C box 2ji .

D. µ − e conversion in nuclei and muonium-antimuonium oscillation
The µ−e conversion process describes the muon capture by a nucleus through the process µ(A, Z) → e(A, Z). At the quark level, the process can be represented as µq → eq, as induced by the photon-and Z-penguin diagrams in the model. Following the results in [15,16], the conversion rate, which is relative to the muon capture rate, is given by: where p e (E e ) is the momentum (energy) of electron and is taken to be m µ in the numerical estimates, Z eff denotes the effective atomic charge of the nucleus, F p is the nuclear matrix element, Γ cap is the total muon capture rate, A ± = Z ±Ñ with Z (Ñ ) being the proton (neutron) number in the nucleus, and g Since the Z-penguin contribution is negligible, the dominant effect is from the photonpenguin diagrams. Thus, the only nonzero g q χK is: The nucleon matrix element G , respectively. Their values can be found in Refs. [15,17] as: For heavy nuclei, because A − A + and g (1) χK , the µ − e conversion rate in the model can be simplified as: In addition to the µ − e conversion and LFV processes, which are ∆L µ = 1 processes, another interesting process involving ∆L µ = 2 (∆L = 2 will be used hereinafter) is the muonium-antimuonium oscillation, where the muonium is a bound state of µ + and e − , i.e., |M µ ≡ |µ + e − . Similar to the µ → 3e decay, the M µ -M µ mixing matrix element can arise from the box diagrams, where the Feynman diagrams are similar to Fig. 3 with ( i , j ) = (µ, e). As in the case of meson oscillations, M µ -M µ mixing effect can be taken as a perturbative effect in quantum mechanics and is formulated by [18] m where m 12 and Γ 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 ∆Γ, we need a resonant intermediate state. In the model, the new particle masses are heavier than the muonium mass m Mµ . As such, ∆Γ by the new effects is negligible, and we concentrate on the mass difference ∆m = m 1 − m 2 ≈ 2 Re(m 12 ).
From Eq. (37), the parameter x ≡ ∆m/Γ µ used to show the probability of M µ → M µ can be written as: Following the formulation obtained in Ref. [18], the oscillation probability is: while the experimental upper limit is: with S B (B 0 ) = 0.75 [18,19] taken in this work. Since the spin-0 para-muonium and spin-1 ortho-muonium are produced in the experiment [19], in order to compare the theoretical estimate with the experimental data, we follow the prescription given in Ref. [18] and take the spin average by combing both spin-0 and spin-1 muonia as: where s i denotes the spin of the muonium M i µ . According to the Yukawa couplings shown in Eq. (10) and the Feynman diagrams in Fig. 3, the effective interaction for the ∆L = 2 process can be written as: where the Wilson coefficient C ∆L=2 1 is given by: Using the transition matrix element M µ |O 1 |M µ = f 2 Mµ m 2 Mµ , the x-parameter for paramuonium and ortho-muonium are: where the reduced mass m red = m µ m e /(m µ +m e ), and f 2 Mµ = 4(m red α em ) 3 /m Mµ [18]. Hence, Eq. (41) can be expressed as: 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 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 , we can approximate the formula of by We can diagonalize the kinetic terms forẐ µ and B µ via the following transformations [20]: where we approximate sin and cos 1 hereafter. After electroweak symmetry breaking, we obtain the mass terms for neutral gauge bosons as where We then obtain mass eigenstates and eigenvalues of the neutral gauge bosons as where we can approximate m Z 1 m Z and m Z 2 m Z for tiny . 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 bosons. The Z interaction with SM fermions f is now given by where T 3 is diagonal generator of SU (2) L , Q is the electric charge operator, S ZZ ≡ sin θ ZZ , and C ZZ ≡ cos θ ZZ .

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.

A. Inputs and constraints
From Eq. (20), it can be seen that when the neutrino mixing angles and masses are determined, the parameters of 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. [21,22] as the inputs, where the global fit results with 3σ errors are given in Table II. 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: 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 α 1,2 ∈ [−π, π]. We then use Eq. (53) to bound the free parameters. The parameter space can be further constrained by various LFV processes, whose current upper bounds are given in Table III. With S B (B 0 ) = 0.75, the upper limit on the probability of the muonium oscillation is taken as P (M µ → M µ ) exp < 11.1 × 10 −11 [18,19].
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 µ → eγ decay shown in Eqs. (22) and (24), we  [24].  [23] assume where ξ is a real parameter and φ is its phase. The parameter ranges used to scan the parameter space are taken as follows: In addition, we also assume m N < m S 1 < m S 2 and m η ± i = m S i in the numerical estimates. Under the constraints of Table II and Table III,

B. 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 µ → eγ 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.   Fig. 6(c), where the dashed line is the sensitivity of COMET [5] and Mu2e [6] experiments. As such, these experiments have the capability to probe most of the considered parameter space through the µ − e conversion rate. The correlation be-tween BR(µ → eγ) and CR(µ → e, Ti) is plotted in Fig. 6(d). Again, when the measurement of CR(µ − e, Ti) ∼ 10 −18 achieves the expected sensitivity, the parameter space with BR(µ → eγ) O(10 −13 ) is mostly covered. Hence, the µ − e conversion process has the potential to become the most stringent constraint among the LFV processes. As stated earlier, in addition to the photon-penguin diagrams, the µ → 3e decay can be induced by the box diagrams in the model. To exhibit the role of the box diagrams, the ratio of BR(µ → 3e) purely from the photon-penguin contribution to BR(µ → eγ) can be expressed as: With C γ 1eµ /C γ 2eµ ∼ O(1), the ratio in Eq. (56) can be estimated to be R 3e/eγ ∼ 5 × 10 −3 . Using the constrained parameter values shown in Fig. 5, indeed, we approximately obtain R 3e/eγ ∈ (0.5, 5) × 10 −2 . We can conclude that if BR(µ → 3e) > BR(µ → eγ) 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 µ → 3e is subleading, it is interesting to compare the µ → 3e decay with the munoium oscillation, where both processes are dominated by the box diagrams. Implementing the constrained parameter values shown in i.e., under the presumption that BR(µ → eγ) ∈ (0.05, 5)10 −13 . We thus see that the two strictest constraints on the µ → e transitions come from the µ − e conversion in nucleus and the µ → 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(µ → 3e) and µ − e conversion rate severely constrain the µ → eγ decay and munoium oscillation processes.
After discussing the rare µ → e transition, we discuss in the following analysis the lepton flavor-changing effects on the heavy τ -lepton decays, where the sensitivities assuming 50 ab −1 of data at Belle II can reach [11]: with = e, µ. Analogous to µ → eγ, the radiative τ decay also arises from the same types of diagrams. Therefore, using Eq. (24) and the constrained parameter values, we plot the correlation between BR(τ → eγ) and BR(µ → eγ) in Fig. 8(a), where the dashed lines are the sensitivities of MEG II and Belle II experiments. Since the resulting branching ratio for τ → eγ is lower than the sensitivity of Belle II, it is difficult to observe the τ → eγ decay in the model. In Fig. 8(b), we show the correlation between BR(τ → µγ) and BR(τ → eγ). It can be found that unlike the τ → eγ decay, the branching ratio for τ → µγ can be as large as O(10 −8 ). As such, τ → µγ serves as a good candidate to probe the new physics effects of the model in Belle II experiment.
As discussed earlier, although the i → 3 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 τ → 3e and τ → 3µ, we show their branching ratio correlations with µ → 3e in Fig. 9. Similar to τ → eγ, it is difficult for Belle II to reach the predicted BR(τ → 3e) in the model. Nevertheless, the τ → 3µ decay is more promising to detect at Belle II because the value can reach O(10 −9 − 10 −10 ).
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 µ in the model contradicts with that given by the recent E989 experiment at Fermilab [25], its value |a µ | < 10 −11 due to the m 2 µ dependence and m η ± Therefore, if the muon g − 2 anomaly can be resolved within the SM [26], the negative a µ  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.

C. 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 exchange with the kinetic mixing effect. In the estimate of relic density, the dominant DM annihilation processes are summarized as follows: where ± 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 [27] implemented with the new interaction vertices in the model.  Fig. 10 show the parameter regions in the m N -g X plane that give the relic density 0.11 < Ωh 2 < 0.13 [24] for cases (1) and (2), respectively. For case (1), the gauge coupling around 0.2 g X 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 400 GeV imposing perturbative condition g X < √ 4π. This behavior is due to that fact that the small kinetic mixing for the s-channel annihilation via Z exchange is suppressed.
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 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 interactions to shift the scatter points downwards by turning on the gauge coupling g X and tuning its value and the Z mass.
Finally we discuss the DM-nucleon scattering cross section via the Z exchange. Here we estimate the cross section in the non-relativistic limit and obtain where n (p) stands for the neutron (proton), µ n(p) ≡ m N m n(p) /(m N + m n(p) ) and C f L(R) is given in Eq. (52). We show in Fig. 12 the DM-neutron scattering cross section as a function of g X for m Z = {50, 100, 250} GeV and m N = 500 GeV. The gray region has Ωh 2 0.11 and the cyan region is excluded by XENON1T [10]. In making this plot, we assume the dominance Z interactions and ignore the Yukawa couplings. The DM-proton scattering cross section is very similar and thus omitted here.

V. SUMMARY
Based upon the concept of radiative neutrino mass in scotogenic model proposed in [1], 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(µ → eγ) can fit the currently experimental upper limit. However, the planned sensitivities of the experiments on the µ − e conversion and the µ → 3e decay in COMET/Mu2e and Mu3e can cover most of the parameter space that can lead to BR(µ → eγ) ∼ 6 × 10 −14 , the designed sensitivity of MEG II experiment. Hence, the µ − e conversion in nucleus and µ → 3e will be the most promising processes to detect the new physics in the µ → e transitions.
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 interaction is dominant (ii) the Yukawa interaction is dominant. In case (i), the observed relic density can be explained with 0.2 g X 0.7 when DM annihilates into Z Z while a larger gauge coupling is required when DM annihilate into the SM particles via the s-channel Z 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 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 mass.