Impact of sterile neutrinos on nuclear-assisted cLFV processes

We discuss charged lepton flavour violating processes occurring in the presence of muonic atoms, such as muon-electron conversion in nuclei $\text{CR}(\mu -e, \text{ N})$, the (Coulomb enhanced) decay of muonic atoms into a pair of electrons BR($\mu^- e^- \to e^- e^-$, N), as well as Muonium conversion and decay, $\text{Mu}-\bar{\text{Mu}}$ and $\text{Mu}\to e^+ e^-$. Any experimental signal of these observables calls for scenarios of physics beyond the Standard Model. In this work, we consider minimal extensions of the Standard Model via the addition of sterile fermions, providing the corresponding complete analytical expressions for all the considered observables. We first consider an"ad hoc"extension with a single sterile fermion state, and investigate its impact on the above observables. Two well motivated mechanisms of neutrino mass generation are then considered: the Inverse Seesaw embedded into the Standard Model, and the $\nu$MSM. Our study reveals that, depending on their mass range and on the active-sterile mixing angles, sterile neutrinos can give significant contributions to the above mentioned observables, some of them even lying within present and future sensitivity of dedicated cLFV experiments. We complete the analysis by confronting our results to other (direct and indirect) searches for sterile fermions.


Introduction
In addition to the several theoretical caveats of the Standard Model (hierarchy and flavour puzzles for instance), three observations clearly signal the need to extend it. New Physics scenarios must be necessarily considered, since the Standard Model (SM) cannot explain the baryon asymmetry of the Universe (BAU), offers no viable dark matter (DM) candidate, and cannot account for neutrino oscillation phenomena (neutrino masses and mixings).
At present, the search for New Physics (NP) is being actively carried in many fronts: direct searches for new states are being pursued at high-energy colliders, and NP is also being indirectly searched for at high-intensity facilities, looking for rare processes or deviations from the SM predictions.
By construction, lepton flavour violation (LFV) is forbidden in the SM. By itself neutrino oscillation phenomena signal the violation of lepton flavour in the neutral sector, and -neutrinos being part of the SM SU (2) L doublets which include the charged leptons -it is only natural to expect that LFV also occurs in the charged lepton sector (cLFV). There are numerous possible manifestations of cLFV, both at high-and low-energies; many observables are currently being searched for in different high-intensity dedicated facilities, as is the case of J-PARC, PSI, SuperB factories, and many others.
In addition to cLFV radiative (ℓ i → ℓ j γ) and three-body (ℓ i → ℓ j ℓ k ℓ k ) decays, rare muon transitions can take place in the presence of nuclei: when a (negative) muon is stopped in matter, it can be trapped, thus forming a "muonic atom". In this work we consider several rare LFV processes occurring in the presence of nuclei (the most common ones being Aluminium, Gold and Titanium), of which the most widely investigated is the coherent µ − e conversion (see e.g. [1]). Numerous experiments have already set strong bounds on the corresponding rate [2][3][4], and the experimental sensitivity to CR(µ − e, N) is expected to significantly improve in the near future [5][6][7][8][9]. Recently, another interesting observable has been proposed [10], which is the Coulomb-enhanced decay of a muonic atom into a pair of electrons, µ − e − → e − e − . The BR(µ − e − → e − e − , N) can be strongly enhanced in the presence of large nuclei [11], and this process will hopefully be part of the physics programme of the COMET experiment [7]. The Muonium [12], bound state composed by a µ + and an e − , constitues a hydrogen-like atom (despite the absence of hadronic nuclear matter). Formed in matter, muonium phenomena are studied in vacuum, and include several cLFV transitions; the study of the latter has the advantage of being free from nuclear uncertainties. Here we will consider the Muonium conversion Mu − Mu [13] and its cLFV decay Mu → e + e − [14].
Many NP models can provide contributions to these observables either by the introduction of new sources of flavour violation and/or the presence of additional degrees of freedom (extensions of the gauge sector and/or of the particle field content) [15]. In this analysis we will focus on the study of cLFV signals arising in minimal extensions of the SM by sterile fermion states.
The sterile neutrino hypothesis is strongly motivated by a number of observations, from an interpretation of ν oscillation anomalies in reactor/accelerator experiments, to an explanation of several cosmology and astrophysics issues (warm DM candidates, pulsar velocities, BAU, etc.) [16]. In addition to the experimental motivation, the sterile fermion hypothesis is very appealing from a theoretical point of view, as these states are an integral part of numerous SM extensions aiming at accounting for the observed neutrino masses and mixings [17]. Indeed, the most minimal extensions of the SM allowing to accommodate neutrino data call upon the introduction of at least two right-handed (RH) neutrinos (SM gauge singlets) providing a Dirac mass term for the neutral leptons, and allowing flavour violation in leptonic charged currents. Should neutrinos be Majorana fermions, then a seesaw-like mechanism emerges as a simple explanation for the smallness of their masses which are thus linked to a NP high-scale via natural Yukawa couplings. Nevertheless, and despite their simplicity, high-scale seesaws are hard to probe, since their impact on low-energy phenomenology is negligible, and the very massive sterile states cannot be produced at colliders. Low-scale seesaw realisations offer far richer experimental prospects: provided that the couplings of the new states to the SM particles are not excessively tiny, they can be probed both at colliders and in high-intensity facilities.
While one can indeed consider theoretically well-motivated models, such as the "standard" type I seesaw mechanism [18], and its many variations (for example the Inverse Seesaw), a first approach to evaluate the phenomenological impact of the sterile states on a given observable is to consider an "ad-hoc" framework, independent of any model of neutrino mass generation. This approach simply relies in adding a single massive Majorana sterile state to the SM, which would encode the effects of a number n S of possibly existing sterile fermion states. No hypothesis is made regarding the mechanism of neutrino mass generation: a leptonic mixing matrix, which must necessarily accommodate oscillation data, provides the only connection between the interaction and the physical states. This "3+1" toy model will be used in our study for a first discussion of the phenomenological impact of the sterile fermion states regarding the nuclear-assisted muonic LFV processes. We will also study the prospects for two well-motivated low-scale seesaw realisations, the Inverse Seesaw [19,20] and the ν-MSM [21][22][23].
For each of the above mentioned frameworks, and having taken into account all available experimental constraints, our analysis reveals that sterile fermions can give significant contributions to several cLFV nuclear-assisted processes; depending on their mass range and on the active-sterile mixing angles, the new contributions can be well within present and future sensitivities of the corresponding cLFV experiments. Concerning the decay of a muonic atom (µ − e − → e − e − ), our phenomenological analysis has shown that it has a strong experimental potential (especially in the case of heavier targets -such as Lead, or Uranium); should its study be pursued at J-PARC [7,24], it would offer a powerful probe into this class of SM extensions, especially in what concerns the Inverse Seesaw realisation here investigated.
Our work is organised as follows: in Section 2 we motivate and describe the three theoretical frameworks used for our study, also summarising the different experimental and observational constraints imposed on the sterile fermions. In Section 3 we discuss the cLFV observables considered, presenting the analytical formulae and the corresponding experimental status. The results arising from the numerical studies (for each of the considered frameworks) are collected and discussed in Section 4. We complete our analysis with a thorough comparison of several observables, confronting our results with the prospects for direct and (other) indirect searches. Our concluding remarks are given in Section 5. The two appendices collect the relevant analytical expressions for quantities (form factors, loop functions, etc.) entering the computation of the cLFV processes under study (Appendix A), as well as a brief summary of the related µ → eee decay (Appendix B).

SM extensions via sterile states
Sterile fermions, such as RH neutrinos (or other fermions which are singlets under the SM gauge group) appear as building blocks of many SM extensions aiming at accounting for neutrino masses and their mixings. The existence of sterile states is further motivated [16] as they might provide an explanation to neutrino oscillation anomalies, or play a rôle in understanding several cosmological observations (being viable warm DM candidates, explaining pulsar velocities, etc.).
In addition to possibly generating Dirac and/or Majorana masses for the light neutrinos, the sterile states can have a non-negligible impact for a number of processes: due to the mixing with the light (mostly active) neutrinos, the sterile fermions induce modifications to the SM charged and neutral currents. If the new sterile states are not excessively heavy, and have sizeable mixings to the light neutrinos, their phenomenological imprint can be important -many observables will thus be sensitive to the active-sterile mixing couplings, and their current experimental values (or bounds) will thus constrain such SM extensions.

Modified leptonic interactions
Let us consider an extension of the SM via n S additional sterile neutral (Majorana) fermions, which have non-negligible mixings with the active neutrinos. In this framework, both charged and neutral current interactions lead to the violation of lepton flavour; the SM boson and scalar interactions with leptons are modified as follows 1 : , where the different terms have been cast in the physical (mass eigenstates) basis; P L,R = (1∓γ 5 )/2, g w denotes the weak coupling constant, cos 2 θ w = 1 − sin 2 θ w = M 2 W /M 2 Z , and m i are the physical neutrino masses (light and heavy). The coefficients C V and C A parametrise the SM vector and axial-vector currents for the interaction of neutrinos with charged leptons, C V = 1 2 + 2 sin 2 θ w and C A = 1 2 . In Eq. (1), the indices α = 1, . . . , 3 denote the flavour of the charged leptons, while i, j = 1, . . . , 3 + n S correspond to the physical (massive) neutrino states.
The mixing in charged current interactions is parametrised by a rectangular 3×(3+n S ) mixing matrix, U αj (which corresponds to the (unitary) PMNS matrix, U PMNS , in the case of only three neutrino generations). The mixing between the left-handed leptons, which we denote byŨ PMNS , corresponds to a 3 × 3 block of U. It proves convenient to parametrise theŨ PMNS in terms of a matrix η [26], which encodes the deviation ofŨ PMNS from unitarity [19,20], due to the mixing between the active neutrinos and the extra fermion states: For the purpose of phenomenological analyses, the invariant quantityη, which is defined as is often used to illustrate the effect of the new active-sterile mixings (corresponding to a deviation from unitarity of theŨ PMNS ). Notice that the extra neutral Majorana fermions can also lead to a violation of lepton flavour in neutral currents. As seen from Eq. (1), the latter effect is encoded in the squared mixing matrix

Some motivated theoretical frameworks
A first, and simple approach, to evaluate the phenomenological impact of the extra sterile fermions is to consider a minimal "toy" model, in which one adds one massive Majorana state to the three light neutrinos of the SM. This model-independent "ad-hoc" construction makes no hypothesis on the underlying mechanism of mass generation, only assuming that the physical and the interaction neutrino bases are related via a 4 × 4 unitary mixing matrix U ij , whose 3 × 4 sub-matrix U lj appears in Eq. (1). Thus, and in addition to the three (light) active masses and corresponding mixing angles, it is only assumed that the leptonic sector contains the following degrees of freedom: the mass of the new sterile state, m 4 , three active-sterile mixing angles θ i4 , two new (Dirac) CP phases and one extra Majorana phase. Although not directly affected by the experimental and observational constraints which will be mentioned in Section 2.3, should the sterile state be sufficiently heavy to decay into a W ± boson and a charged lepton, or into a light (active) neutrino and either a Z or a Higgs boson, its decay width is indirectly bounded from the fact that the decays should comply with a perturbative unitarity condition [27][28][29][30][31][32], . Since the dominant contribution arises from the charged current term, the following bound on the heavy sterile masses and their couplings to the active states [27][28][29][30][31][32] is obtained: where α w = g 2 w /4π, and C ii is given in Eq. (4). In the numerical analysis carried in Section 4, we will consider a normal ordering (NH) for the light neutrino spectra 2 ; concerning the extra sterile state, we scan over the following mass range while randomly varying the active-sterile mixing angles in the interval [0, 2π] (ensuring that the condition of Eq. (5) is respected). All the CP-violating phases are also taken into account, and likewise randomly varied between 0 and 2π.
As mentioned before, sterile fermions are an integral part of several mechanisms of neutrino mass generation, and appear in the matter field content of many SM extensions (type I seesaw and its variants, Left-Right symmetric models, GUTs). In our study we will also illustrate the impact of sterile fermions on cLFV rare processes in two well motivated theoretical models: the Inverse Seesaw (ISS) mechanism [33] and the Neutrino Minimal SM (νMSM) [21], both characterised by a low seesaw scale (well below the GUT scale). Both models have a rich phenomenological and cosmological impact, and have been subject to extensive studies in recent years.
The ISS mechanism [33] extends the SM via the addition of both RH and sterile neutrinos, and allows to account for neutrino data with (almost) natural values of the Yukawa couplings and for a comparatively low seesaw scale. Depending on its actual realisation, the ISS does allow to accommodate the observed DM relic abundance and (potential) indirect DM detection hints [34,35], also providing a framework where the observed BAU can be generated via leptogenesis [36].
We consider here a specific ISS realisation in which n R = 3 generations of RH neutrinos and n X = 3 generations of extra singlet fermions X are added to the SM content, both carrying lepton number L = +1 [33]. The modified leptonic SM Lagrangian is given by where i, j = 1, 2, 3 are generation indices andH = iσ 2 H * . After electroweak (EW) symmetry breaking, the neutral lepton spectrum is composed of three (mostly active) light states, while the (mostly) sterile states form three nearly degenerate pseudo-Dirac pairs. The interaction and the physical bases are related by a 9 × 9 unitary mixing matrix U, which diagonalises the 9 × 9 full neutrino mass matrix as U T MU = diag(m i ). In the basis where the charged lepton mass matrix is diagonal, the leptonic mixing matrix is given by the rectangular 3 × 9 sub-matrix corresponding to the first three columns of U, with the 3 × 3 block corresponding to the (non-unitary)Ũ PMNS . In order to derive the contributions to the studied observables, a scan is carried over the 9 × 9 neutrino mass matrix (for a detailed discussion of the numerical studies, see [37]): the modulus of the entries of the M R and µ X matrices are randomly taken to lie on the intervals 0.5 GeV (M R ) i 10 6 GeV and 0.01 eV (µ X ) ij 1 MeV, with complex entries for the lepton number violating matrix µ X ; a modified Casas-Ibarra parametrisation [38] for Y ν allows to accommodate neutrino oscillation data (we take complex angles for the R matrix accounting for the additional degrees of freedom, randomly varying their values 3 in the interval [0, 2π]; we further verify that the Yukawa couplings are perturbative, i.e. Y ν < 4π). We again consider a NH for the light neutrino spectrum.
The νMSM consists in a truly minimal extension of the SM via the inclusion of three RH neutrinos, aiming at simultaneously addressing the problems of neutrino mass generation, the BAU and providing a viable DM candidate [21-23, 39, 40].
The addition of three generations of RH Majorana states ν R to the SM particle content allows to add the following terms to the leptonic Lagrangian: where i, j = 1, 2, 3 are generation indices, L is the SU (2) L lepton doublet andH = iσ 2 H * ; Y ν denotes the Yukawa couplings, while M M is a Majorana mass matrix (leading to the violation of total lepton number, ∆L = 2). In addition to the three light (mostly active) neutrinos -whose masses are given by a type I seesaw relation -the neutral lepton spectrum further contains three heavy states (with masses m ν 4−6 ); should the latter three play a rôle regarding DM, or the BAU, their spectrum and couplings to the active states are subject to several (strong) conditions [39][40][41]. As previously done in [37], we parametrise the physical νMSM degrees of freedom in terms of its six mass eigenvalues, encoding all physical mixing angles and CP violating phases (Dirac and Majorana) in an effective 6 × 6 unitary mixing matrix U (which allows to readily implement the already well-established bounds on the νMSM parameter space) 4 . The angles θ lj , l = 1, 2, 3, j = 4, 5, 6 encode the active-sterile mixings, while the mixings between the sterile states are given by three additional angles θ 45,46,56 . The matrix U is further parametrised by 3 additional Majorana and 9 Dirac phases. As before, we will only consider a NH for the light neutrino spectrum. 3 The choice of a larger interval for the complex angles would increase the range of the entries of Y ν , and thus would lead to augmented contributions to the considered cLFV observables; however, this would also increase the other (muonic) cLFV observables which, given the present bounds, would in turn lead to the exclusion of these regimes. 4 In our numerical analysis, we rely on the results of [39], where the most relevant constraints are translated into bounds on the (U 2 , M ) planes, as well as on the splitting δM between the two heaviest states (in the range from ∼ 10 −4 eV to 1 keV). The quantity U 2 encodes the experimentally relevant combination of couplings:

Constraints on sterile fermions
The potentially sizable contributions to a number of processes, induced by the mixings of the sterile states with the active neutrinos, implies that several observables (low-energy, collider and cosmological) can severely constrain these models. Firstly, one ensures that the SM extension complies with ν-oscillation data [42][43][44][45][46][47]; in all the scenarios considered in this work, and for both cases of NH and inverted mass hierarchy (IH), we require compatibility with the corresponding best-fit intervals [45] (no constraints being imposed on the yet undetermined value of the CP violating Dirac phase δ).
In our analysis we also apply -when applicableunitarity bounds on the (non-unitarity) matrix η (cf. Eq. (2)); these arise from non-standard neutrino interactions with matter, and have been derived in [48][49][50] by means of an effective theory approach (valid for sterile masses above the GeV, but below the electroweak scale, Λ EW ).
The addition of sterile states to the SM with a sizeable active-sterile mixing may have an impact on electroweak precision observables either at tree-level (charged currents) or at higher order. We take into account the impact of sterile neutrinos on the invisible Z-decay width (which has been addressed in [51][52][53][54]), requiring compatibility with LEP results on Γ(Z → νν) [55]; in addition, we further require that potential new contributions to the cLFV Z decay width do not exceed the present uncertainty on the total Z width [55]: Γ(Z → ℓ ∓ 1 ℓ ± 2 ) < δΓ tot . LHC data on invisible Higgs decays already allows to constrain regimes where the sterile states are below the Higgs mass; in our study we apply the constraints derived in [56][57][58]. Negative laboratory searches for monochromatic lines in the spectrum of muons from π ± → µ ± ν decays [59,60] also impose robust bounds on sterile neutrino masses in the MeV-GeV range.
The introduction of singlet neutrinos with Majorana masses allows for new processes like lepton number violating interactions, among which neutrinoless double beta decay remains the most important one [61]. In our analysis, we evaluate the contributions of the sterile states to the effective mass m ee according to [62,63]; strong bounds on the effective mass have been set by several experiments, among them GERDA [64], EXO-200 [65,66], KamLAND-ZEN [67]. In our analysis we use the most recent constraint from [66].
Further constraints arise from leptonic and semileptonic decays of pseudoscalar mesons K, D, D s , B (see [68,69] for kaon decays, [70,71] for D and D S decay rates, and [72,73] for B-meson observations). Recent studies suggest that in the framework of the SM extended by sterile neutrinos the most severe bounds arise from the violation of lepton universality in leptonic kaon decays (parametrised by the observable ∆r K ) [54,74].
Other than the rare decays occurring in the presence of nuclei, the new states can contribute to several charged lepton flavour violating processes such as ℓ i → ℓ j γ and ℓ i → ℓ j ℓ k ℓ k . In our analysis we compute the contribution of the sterile states to all these observables [20,25,[75][76][77][78][79][80], imposing compatibility with the bounds summarised in Table 1, also considering the impact of the future experimental sensitivities.
Finally, a number of cosmological observations [59,88] put severe constraints on sterile neutrinos with a mass below the GeV. While very light sterile neutrinos (with a mass eV) appear to be disfavoured by CMB analysis with the Planck satellite [89], a sterile state with mass ∼keV could be a viable DM candidate, also offering a possible explanation for the observed X-ray line in galaxy clusters spectra at an energy ∼ 3.5 keV [90,91], for the origin of pulsar kicks, or even to the BAU (for a review see [17]). Independently of the model under consideration, and for all regimes of sterile masses investigated, in our phenomenological analysis we ensure that all the above referred constraints - Table 1: Current experimental bounds and future sensitivities for radiative and 3-body cLFV decays considered in our study.
theoretical (such as perturbativity of the active-sterile couplings) and experimental -are verified.
Since the cosmological bounds are in general derived by assuming the minimal possible abundance (in agreement with neutrino oscillations) of sterile neutrinos in halos consistent with standard cosmology, a non-standard cosmology with a very low reheating temperature or a scenario where the sterile neutrinos couple to a dark sector [92], could allow to evade some bounds [93]. In our numerical analysis, we adopt a conservative approach, allowing for the violation of these cosmological bounds in some scenarios (which will be identified).

Muonic atoms and cLFV rare processes
In what follows, we present the different cLFV observables which can be studied in relation to rare processes involving muonic atoms.

Muon-electron conversion
One of the most sensitive probes of cLFV is the µ − e conversion occurring in a muonic atom, which is formed when a muon is "stopped", falling into the 1s state of a target nucleus. The observable is defined as The coherent conversion (in which the nuclear final state is in its ground state) increases 5 with the atomic number (Z) for light nuclei, Z 30, and is maximal for 30 Z 60 [94]. For heavier nuclei, Coulomb distortion effects of the wave function lead to a reduction of the corresponding conversion rate. Past and present experiments have mostly explored Titanium, Lead and Gold nuclei; the present bounds for different targets (obtained by the SINDRUM II experiment) are summarised in Table 2.
At present, two projects are aiming at improving the current bounds, Mu2e (at Fermilab) [5] and COMET (J-PARC) [6,7], sharing several common features (nominal muon beam energy ∼ 8 GeV, Aluminium targets, ...). The COMET experiment will carry a two-phase search for µ − e conversion, with the second phase expected to bring the sensitivity down by two orders of magnitude with respect to Phase I. (It is also possible that different targets will be used at this stage [24].) The COMET experiment will be subsequently modified: PRISM/PRIME will further Present CR(µ − e, N) bound material year 4.3 × 10 −12 [2] Ti 1993 4.6 × 10 −11 [3] Pb 1996 7 × 10 −13 [4] Au 2006 improve the sensitivity to these cLFV processes [8]. Another possibility is that of DeeMe [9], which is a parallel project at J-PARC, albeit on a smaller experimental scale: working with a siliconcarbide target, they aim at lowering the present sensitivity by almost two orders of magnitude.
In the long term, Project-X (at Fermilab) [95,96] is expected to benefit from very intense muon beams, with at least ten times more muons than Mu2e (however for a lower energy beam, around 1-3 GeV). We summarise the expected future sensitivities in Table 3. Several models of NP can give rise to the rare cLFV nuclear conversion, and the corresponding rates have been extensively discussed in the literature [1,15,97]. For the present class of SM extensions, contributions to the nuclear µ − e conversion arise from the Z-and photon-penguin diagrams, as well as boxes, which are schematically depicted in Fig. 1. At lowest order, the flavour violating µ−e transition originates from one-loop diagrams involving neutrinos (active and sterile); their non-zero masses and mixings prevent a GIM cancellation (which would otherwise occur).
The interactions depicted in Fig. 1 can be described by the following effective Lagrangian [79] (neglecting the electron mass), where G µe γ refers to the photon-lepton dipole coupling corresponding to an on-shell photon ("non local", long-range contribution), F λρ = ∂ λ A ρ − ∂ ρ A λ denotes the electromagnetic field strength, andF µe q contains the "local" contribution of the monopole F µe γ , as well as those arising from the weak gauge-boson exchange diagrams. The dependence of the process on the new degrees of freedom associated with the sterile neutrinos are contained in the form factorsF µe q and in the dipole term G µe γ . In the SM extended by sterile neutrinos, the ratio for the nuclear assisted µ − e conversion, see Eq. (9), can be cast in the following compact form, (11) In the above, Γ capt (Z) denotes the capture rate of the nucleus characterised by the atomic number Z [94]; G F is the Fermi constant, m µ the muon mass, α = e 2 /(4π), with s w corresponding to the sine of the weak mixing angle. The form factorsF µe q (q = u, d) are given bỹ where Q q denotes the quark electric charge (Q u = 2/3, Q d = −1/3) and I 3 q is the weak isospin (I 3 u = 1/2 , I 3 d = −1/2). The quantities F µe γ , F µe Z and F µeqq Box correspond to the different form factors of the diagrams depicted in Fig. 1, respectively photon-penguins, Z-penguins and box diagrams; their expressions are collected in Appendix A, as is the one concerning the dipole term, G µe γ . The relevant nuclear information (nuclear form factors and averages over the atomic electric field) are encoded in the D, V (p) and V (n) form factors. In our analysis we use the numerical values presented in [94].
We conclude this discussion commenting about another interesting nuclear-assisted LFV process, which is a charge-changing reaction of the type [1,98]: denoted µ − → e + conversion in nuclei. The final nucleus, which is different from the initial onethus preventing a coherent enhancement -can be either in its ground state or in an excited one. This rare process violates the conservation of the total lepton number by two units and is related to neutrinoless double decay process. Although the sterile neutrinos may give a non negligible contribution to the µ − → e + conversion rate, we do not address this observable here.

Decay of muonic atoms to
A different cLFV process was recently proposed in [10]. It consists in the flavour violating decay of a bound µ − in a muonic atom into a pair of electrons, and has been identified as potentially complementary to other cLFV muon decays: In the above transition, the initial states are a µ − and a 1s atomic e − , bound in the Coulomb field of a nucleus [10].
Although the underlying sources of flavour violation giving rise to this observable are the same as those responsible for other non-radiative µ−e transitions (such as µ → eee, or Muonium decay), the µ − e − → e − e − decay in a muonic atom offers significant advantages. From an experimental point of view, and in addition to having a larger phase space, it leads to a cleaner experimental signature than the 3-body µ + → e + e − e − decay; the LFV muonic atom decay has a two-body final state in which the electrons are emitted nearly back-to-back, each with a well-defined energy 6 , E e − ∼ m µ /2. Furthermore, the rate of this process can be enhanced due to the Coulomb attraction from the nucleus, which increases the overlap of the 1s electron and muon wavefunctions, ψ (e),(µ) 1s . Following [10], one can write the transition rate of the µ − e − → e − e − process as where σ µe→ee v rel denotes the electroweak cross-section. Thus, and when compared to the (apparently) similar Muonium decay µ + e − → e + e − , the effects of the Coulomb interaction lead to an enhancement of the muonic atom decay rate by a factor ∼ (Z − 1) 3 , which can become important for nuclei with large atomic numbers.
The most recent results of [11] emphasised the importance of taking into account the distortion effect of the emitted electrons due to the nuclear Coulomb potential, and of performing a relativistic treatment of the wave function of the bound leptons. Whilst for small atoms the enhancement is indeed well described by f Coul. (Z) ≈ (Z − 1) 3 , for large Z atoms the rate can be enhanced by as much as an additional order of magnitude (see Fig. 1 of [11]). In our study we take into account these results (in particular we use the data obtained for a uniform distribution of the nuclear charge).
As was the case for the coherent conversion in nuclei, in extensions of the SM via sterile fermions, the effective Lagrangian describing the muonic atom LFV decay contains γ-dipole (longrange) interactions and "local" (contact) terms [1,10]. The contribution to the cross-section of the (long-range) photonic interactions is subdominant for the present class of NP models 7 ; moreover, in the absence of a complete estimation of the associated nuclear effects [24,99], we will not take into account the "long-range" γ-dipole contributions. The contact interactions include contributions from box diagrams, as well as photon and Z penguin diagrams; in a muonic atom, with an atomic number Z, they contribute to the branching ratio of the process in Eq. (14) as where we have neglected the interference between contact terms (which can be sensitive to CP violating phases that might emerge in relation to flavour violating effects). In the above equatioñ τ µ corresponds to the lifetime of a muonic atom, which depends on the specific element, and is always smaller than the lifetime of free muons, τ µ (τ µ = 2.197 × 10 −6 s [55]). In Table 4 we summarise the values ofτ µ for different elements [100], which will be relevant in our discussion.  In Eq. (16), F µe γ,Z denote the contributions from photon-and Z-penguins (corresponding to those already introduced in Eq. (11) of the previous subsection); F µeee Box is defined in Appendix A. As can be inferred from the nature of the contributing diagrams, the muonic atom µ − e − → e − e − decay arises from the same elementary processes as the rare 3-body decay µ + → e + e − e − . Together with the bounds arising from the coherent conversion in Nuclei (CR(µ−e, N)), the present experimental bounds on BR(µ → eee) may indirectly constrain the possible values of BR(µ − e − → e − e − , N). Thus, and for completeness, we collect the relevant expressions for the 3-body muon decay in Appendix B.
As mentioned in the Introduction, the µ − e − → e − e − process could be investigated by the COMET collaboration (possibly being part of its Phase II programme). In our numerical analysis, we will thus work under the hypothesis of a similar future sensitivity for BR(µ − e − → e − e − , Al) and CR(µ − e, Al).

Muonium: Mu-Mu conversion and decay
The Muonium (Mu) atom is a Coulomb bound state of an electron and an anti-muon (e − µ + ) [12], formed when a µ + slows down inside matter and captures an e − . This hydrogen-like state is actually free of hadronic interactions, and its electromagnetic binding is well described by the SM electroweak interactions. In turn, this renders the Muonium system an interesting laboratory to accurately determine fundamental constants, or test for deviations from the SM induced by the presence of possible NP states and interactions.
The spontaneous conversion of a Muonium atom to its anti-atom (Mu = e + µ − ) has been identified as an interesting class of muon cLFV process [13]. Similar to the other cLFV observables, It is worth noticing that the new interactions will cause a splitting of the energy levels of Muonium and antimuonium (which, in the absence of an external magnetic field, would otherwise have degenerate ground state energy levels). This splitting can be also cast in terms of the effective coupling as [1] δ where n and a 0 denote the principal quantum number and the Bohr radius of the muonic atom. Searches for Mu-Mu conversion were started more than fourty years ago, and have employed different techniques. At present, no positive signal has been found; the best limit has been set at PSI [101], where Muonium atoms are formed by electron capture when positively charged muons (from a very intense muon beam) are stopped in a SiO 2 powder target. Experimental searches put a bound on Re G MM , or on its ratio to G F , and the limit(s) are usually determined assuming that the undelying interaction is of the type (V ± A) × (V ± A); for these cases, the 90%C.L. limit on the effective coupling constant has been reported to be Re G MM ≤ 3.0 × 10 −3 G F [101].
On Table 5 we summarise the different bounds so far obtained for the conversion probability P(Mu-Mu), translated into upper bounds on the effective coupling G MM . With the advent of new, very intense muon sources, it can be expected that the current bounds will be improved in the near future.
In extensions of the SM with sterile neutrinos, the e − µ + → e + µ − transition responsible for Mu-Mu conversion occurs at the loop-level via box diagrams. There are four different diagrams which can be separated in two contributions, those which are common to both Dirac and Majorana neutrinos (Fig. 2, upper boxes) and two other which appear only if neutrinos are Majorana particles (Fig. 2, lower boxes). The computation of the box diagrams of Fig. 2 (in a unitary gauge) allows to write the effective coupling G MM as [105,106]: where x i = This observable has been addressed in the context of several extensions of the SM via additional sterile states, mostly in relation with type I seesaw realisations with heavy RH neutrinos [105][106][107].
Although a Muonium state mostly decays via standard channels (the most frequent one being Mu → e + e −ν µ ν e , or equivalently, a SM muon decay with the electron as a spectator), the presence of NP can also induce the cLFV decay In SM extensions with RH neutrinos, these decays have been addressed at low-energies, and also in the framework of a future muon collider [14]. The cLFV Muonium decay rate can be written as with Γ µ the muon decay width, and where |M tot | denotes the full amplitude, summed (averaged) over final (initial) spins [14]. The full expression for |M tot | is given in Appendix A; we notice that the intervening quantities are those also present for 3-body lepton decays (although the distinct form factors -penguins, long-range dipoles, boxes -contribute differently to the matrix element). The experimental roadmap concerning this last cLFV observable is not clear; we will thus make no reference to experimental limits or sensitivities in the numerical analysis.

cLFV in muonic atoms: numerical results
We now investigate the cLFV observables presented in the previous section, for different classes of SM extensions: we first carry a detailed analysis for a simple toy model ("3+1"), and then highlight the most important points emerging from the study of two well motivated NP models, the ISS and the νMSM.
4.1 "3+1 model" toy model: nuclear-assisted cLFV processes As described in Section 2.2, this simple model allows to illustrate the potential effects of the addition of n S sterile fermion states, "encoded" in the contribution of ν 4 .

µ → e µ → e µ → e conversion in nuclei
We begin our study by considering the impact of the additional sterile state concerning the coherent µ → e conversion in a muonic atom. We choose to illustrate this observable for the specific case of an Aluminium nucleus; similar results would be obtained for other nuclei (as Gold or Titanium) 9 . The results of a comprehensive scan over the additional degrees of freedom of the sterile state (as described in Section 2.2) are displayed in Fig. 3 in which, for completeness, we also include the predictions of the simple "3+1 toy model" for the three body µ → eee decay. In grey we denote points violating at least one of the experimental constraints listed in Section 2.3, with the exception of the observables under study. Coloured points denote the two observables: in dark blue the predictions for CR(µ − e, Al), corresponding to the left y-axis; in cyan, the values of BR(µ → eee), displayed on the right y-axis. (Cosmological bounds would have typically disfavoured regions in parameter space for which m 4 0.1 GeV, and hence are not visible in Fig. 3.) Minimal extensions of the SM via sterile fermions -such as the simple "3+1 toy model"can easily account for sizable contributions to both CR(µ − e, Al) and BR(µ → eee), even above the current experimental bounds (respectively depicted by thick and thin solid horizontal lines in Fig. 3). This is in good agreement with previous analyses carried for specific models (e.g. low-scale type I seesaw [79], ISS [80]). For the mass regime of the mostly sterile state around the EW scale (m 4 ∼ 10 2 GeV), the leading contributions arise from Z-penguin diagrams; below the EW scale, box diagrams become increasingly important, and dominate the total width below a few GeV. (For the low mass regime, points are excluded as neutrino data cannot be accommodated, mostly due to an excessive departure from unitarity of theŨ PMNS ; bounds from neutrinoless double beta decays are also important in this mass regime.) For very large values of the mostly sterile heavy neutrino mass, m 4 10 TeV, the bound arising from CR(µ−e, Au) becomes the most constraining one.
As mentioned before, the scan leading to the results displayed in Fig. 3 explores all the new degrees of freedom of the additional sterile state, in particular the possible extra CP violating phases. We have verified that, as expected, these do not play a significant rôle for this type of observables. We have also considered an IH spectrum for the light neutrino spectrum, finding that this leads to similar results concerning the cLFV observables (the only significant difference being that bounds from neutrinoless double beta decay now exclude more important regions of the parameter space). The former is displayed in dark blue (left axis), while the latter is depicted in cyan (right axis). Grey points correspond to the violation of at least one experimental bound (other than those arising from CR(µ − e, Au) and BR(µ → eee)). A thick (thin) solid horizontal line denotes the current experimental bound on the CR(µ − e, Au) [4] (µ → eee decays [85]), while dashed lines correspond to future sensitivities to CR(µ − e, N) [7,8], see Tables 2 and 3.
The coherent µ − e conversion in muonic atoms, induced by an additional sterile neutrino, could certainly be probed in near future experiments, as Mu2e or COMET (both with Aluminium targets). In fact, for masses of the heavy (mostly) sterile state above the EW scale, the predictions of this simple model are well within reach of the first COMET phase -whose sensitivity corresponds to the upper horizontal dashed line in Fig. 3 (see Table 3) -or even DeeMe, which is expected to reach an O(10 −14 ) sensitivity, albeit for a silicon-carbide target.

Decay of muonic atoms to e − e − pairs
Despite the apparent similarity to the µ → eee decay, the (Coulomb enhanced) decay of a muonic atom to a pair of electrons substantially differs from the 3-body decay, both at the theoretical and at the experimental level.
Firstly, and as emphasised in the discussion of Section 3.2, the process' rate can be significantly enhanced in large Z atoms (in particular the contributions from contact interactions [11]). We thus begin the numerical analysis by comparing the prospects for two different nuclei; this is illustrated in Fig. 4 for Aluminium (Z = 13, dark blue) and Uranium (Z = 92, cyan). Grey points correspond to the violation of at least one experimental bound: the most stringent constraints arise, as expected, from µ → eee (and also from CR(µ − e, Au)).
The Coulomb enhancement is clearly visible: should this process be included in COMET's physics programme, the cLFV muonic atom decay should be within reach of COMET's Phase II, even for light nuclei, such as Aluminium (in the regime m 4 200 GeV); for heavier atoms, such as Uranium, branching ratios above 10 −15 render this process experimentally accessible (a similar situation occurs for Lead nuclei -albeit suppressed by a factor ∼ 7/9 when compared to Uranium [11]).
As mentioned in Section 3.2, in the absence of a complete estimation of the nuclear effects regarding the long-range photon cLFV interaction, we only considered contributions from contact interactions. It is possible that the additional dipole interactions further increase the total BR(µ − e − → e − e − , N) [11,99].
Should the decay of the muonic atom (µ − e − → e − e − ) be included in the physics programme of dedicated high-intensity facilities (as COMET), then it is only natural to investigate to which extent it can become a powerful probe of cLFV, and how it can complement information obtained from other nuclear assisted processes, such as µ − e conversion. We thus display in Fig. 5 the expectations for both observables, BR(µ − e − → e − e − , Al) and CR(µ − e, Al) in the simple "3+1" toy model. In general, the coherent conversion appears to have a stronger experimental potential, with contributions well within COMET reach for sterile masses above a few tenths of GeV. Interestingly, for the low-mass regime, the sterile fermion contributions to the muonic atom decay (via the box diagrams, as already discussed before) could be much larger, but this regime is heavily constrained on theoretical and experimental arguments. Heavier atoms, such as Lead, would further enhance these contributions.
To conclude the discussion of this observable, we present in Fig. 6 the predictions for the logarithm of the BR(µ − e − → e − e − , Al) in the (|U µ4 | 2 , m 4 ) parameter space of the simple "toy model". Shaded surfaces reflect the violation of at least one phenomenological or cosmological bound (see Section 2.3, as well as [108], [109]). The different solid lines correspond to the reach of future facilities: the projected exclusion limit from the LHC (14 TeV centre of mass energy, with 300 fb −1 data [109]); the expected sensitivity of FCC-ee regarding the production of heavy (RH) sterile neutrinos from Z → ν ℓ ν s (estimated 10 12 Z decays, for a 10-100 cm decay length [110]); DUNE (former LBNE, a beam dump experiment, searching for the decay products of sterile neutrinos produced in charmed meson decays) [111]; SHiP (a fixed-target experiment using highintensity proton beams at the CERN SPS [112,113]). As can be seen, there is little overlap between points associated with BR(µ − e − → e − e − , N) within COMET sensitivity and the reach of the above mentioned facilities. This is mostly a consequence of the large mass regime responsible for the latter contributions, which precludes the (direct) production of the sterile states. The potential sensitivity of COMET concerning this observable (which explores cLFV in the µ − e sector) would be complementary to that of cLFV in the µ − τ sector, in particular concerning high-energy observables such as the decay Z → µτ , which could be studied at FCC-ee (TLEP), running in e + e − mode close to the Z mass threshold. In fact, most of the points within COMET sensitivity are predicted to account for a BR(Z → µτ ) lying within FCC-ee reach ( 3 × 10 −13 ) [37]. This implies that this simple "3+1" toy model can be probed via two fully independentand yet strongly complementary -cLFV experimental approaches.

Muonium oscillation and decays
As discussed in Section 3.3, Muonium consists of a hydrogen-like atom (although free of hadronic interactions). We proceed to summarise the prospects regarding the sterile neutrino contributions to Muonium-antimuonium transitions and Muonium decay. Figure 7 displays the expected contributions for the simple "3+1" toy model, both to the Mu -Mu conversion (left) and to the cLFV Muonium decay Mu → e + e − (right). For completeness, the colour scheme of this figure illustrates the regimes which would be excluded from violation of cosmological bounds (red).
Although the sterile contribution could potentially account for values of the effective coupling constant not too far from the most recent measurement, these points are excluded as they would be associated to excessive CR(µ−e, Au) and BR(µ → eee). The phenomenologically viable parameter    Table 5.

Low-energy seesaw models
To conclude our discussion, we illustrate the contributions to the "nuclear assisted" cLFV observables arising in the framework of two well-motivated low-energy seesaw models.

The (3,3) Inverse Seesaw realisation
We first consider the ISS realisation in which three RH neutrinos and three additional steriles are added to the SM content. In Fig. 8, we display the predictions of the ISS concerning the Coulomb enhanced BR(µ − e − → e − e − , N), comparing it with the CR(µ − e, N). We consider again the case of Aluminium targets, and display the results as a function of the average mass of the heavier (mostly sterile) states, < m 4−9 > = 1 6 i=4...9 |m i | .
As can be seen, both observables are well within experimental reach for spectra containing (at least) one pseudo-Dirac pair even below the EW scale (concerning the coherent µ − e conversion) and above the TeV scale (for the muonic atom decay rate). Due to the contributions of the six sterile states, one finds predictions for both observables as large as those arising in the case of the simple "3+1" toy model; in particular, BR(µ − e − → e − e − , Al) can reach values above 10 −16 , thus within reach of COMET.
Remarkably, and despite the intrinsic constraints on its Yukawa couplings, this (3,3) ISS realisation offers the following possibilities: (i) a signal of CR(µ − e, Al) observable at COMET's Phase I; (ii) BR(µ − e − → e − e − , N) within the sensitivity of COMET's Phase II even for an Aluminium target. The ISS mass regions leading to the above signals can be further probed via complementarity studies, if one again considers cLFV from the µ − τ sector, and its potential observation at a future collider. As was already the case for the simple "3+1" toy model, the (3,3) ISS regions with a BR(µ − e − → e − e − , Al) 10 −17 would also induce BR(Z → µτ ) within reach of the FCC-ee [37].
A general overview of the (3,3) ISS realisation here studied is given in Fig. 9, in which we display the logarithm of the BR(µ − e − → e − e − , Al) in the (|U ℓi | 2 , m i ) parameter space of two states: one belonging to the lighest pseudo-Dirac pair (ν 5 ), and another to the heaviest pair, ν 9 . In addition to the expected sensitivity reach of the experiments already mentioned in Section 4.1.2, two additional curves are depicted in the panels of Fig. 9: the future linear collider (ILC) expected sensitivity (solid blue) in (|U ei | 2 , m i ) panels, which has been obtained assuming a √ s = 500 GeV and luminosity of 500 fb −1 [109,114]; the conservative and optimistic (solid red and dashed red respectively) projected limits at 90% C.L. in (|U τ i | 2 , m i ) panels, from semileptonic tau decays at a future B-factory [115].
As can be seen on the different panels of the left column, the states of the lightest pseudo-Dirac pair (ν 4,5 ), belonging to a sterile spectrum responsible for sizable BR(µ − e − → e − e − , Al), are well within reach of the different future facilities which are directly looking for these heavy sterile states. On the right column, one finds the summary of the heavier state's properties (mass and couplings -displaying U 2 ℓ9 on the different panels); this confirms the information of Fig. 8 -one finds diagonal bands corresponding to increasing values of the BR(µ − e − → e − e − , Al), for larger values of the masses and couplings. In this case, only a small subset of the heavier states could be directly probed at FCC-ee through Z → ν ℓ ν 9 decays. However, and as previously stressed, many of these states can be indirectly probed at the FCC-ee through the LFV Z → µτ decays. For the next-to heaviest states, ν 6,7 , one encounters an intermediate scenario, with a significant subset of the points within reach of future facilities. Fig. 10, the prospects for the observation of the |∆L e,µ | = 2 Muonium conversion (left) and of the cLFV Muonium decay (right) are similar to those identified in the simple "3+1" model; the phenomenologically allowed regions of the parameter space would typically lead to Re G MM 10 −13 and to BR(Mu → e + e − ) ∼ 10 −25 . Figure 11 summarises the prospects of the νMSM concerning the decay of muonic atoms to e − e − pairs. The results are displayed in the (U 2 µ , M ) parameter space, generated by the mass of the lightest (mostly sterile) state, M ≡ m 4 , and the mixings of the additional sterile states to the muon-neutrino,

The νMSM
At most, and in the case of Aluminium targets, one can expect BR(µ − e − → e − e − , Al) ∼ 10 −25 for a small subset of the investigated points; for heavier targets, such as Uranium, corresponding to   the results displayed in Fig. 11, values of BR(µ − e − → e − e − , U) ∼ 10 −22 can be reached. Although not displayed here, the rates for the coherent µ − e conversion were also found to be very small. This implies that the νMSM's parameter space is beyond experimental sensitivity concerning this class of nucleus-assisted cLFV observables. In fact, the smallness of the contributions to cLFV observables (among them those here studied) is characteristic of the νMSM given the smallness of the associated neutrino Yukawa couplings. Concerning Mu-Mu oscillations and cLFV Muonium decay, the νMSM's predictions are also extremely tiny: Re G MM 10 −26 and BR(Mu → e + e − ) ∼ 10 −35 .

Conclusions
In this work we have investigated the impact of sterile fermions on cLFV observables which occur in the presence of "muonic atoms": coherent µ − e conversion, the (Coulomb enhanced) decay of muonic atoms into e − e − pairs as well as Mu-Mu oscillations and cLFV Muonium decays. We have considered extensions of the SM which add one or more sterile neutrinos to its particle content: these include a simple "3+1" toy model or well motivated frameworks for the generation of neutrino masses, for instance low-scale seesaw mechanisms like the ISS and the νMSM. Due to their non-negligible mixing with the lighter (mostly active) neutrinos -which induces a departure from unitarity of theŨ PMNS -, the sterile states can have a significant phenomenological impact. In particular, they can provide non-negligible contributions to a number of cLFV observables which are being actively searched for. In turn, the new states are subject to a vast array of observational, experimental and theoretical constraints, which must be taken into account.
Our analysis confirmed that minimal extensions of the SM, such as the "3+1" toy model, can account for sizeable contributions to CR(µ − e, N) and also to muon three body decays, µ → eee (which is in general correlated with µ − e conversion). In fact, important regions of the parameter space can be probed in near future experiments such as Mu2e, DeeMe and COMET, while other regions are actually already excluded by current experimental limits on CR(µ − e, Au) and BR(µ → eee). As expected, a similar scenario arises in the context of the ISS, which also predicts that a signal of CR(µ − e, Al) could be potentially seen at COMET in its first phase.
The decay of a muonic atom into a pair of electrons has proven to be another powerful probe for cLFV, particularly sensitive to the presence of sterile neutrinos. While similar to µ → eee from the point of view of the elementary flavour violating transitions, this process has its rate strongly enhanced by Coulomb interactions between the charged leptons and the electromagnetic field of the nucleus (augmenting with increasing Z). The experimental relevance of this observable is manifest even for the simple "3+1" toy model: sterile neutrinos with masses m 4 > ∼ 800 GeV, lead to BR(µ − e − → e − e − , Al) within the reach of COMET, and the contributions would be further enhanced for heavier atoms, such as Lead or Uranium, thus improving the experimental potential.
The ISS also offers an interesting scenario concerning this observable. Not only the added contributions of the extra six sterile neutrinos could give rise to large BR(µ − e − → e − e − , Al) -as large as 10 −16 -but for an ISS spectrum containing at least one heavy pseudo-Dirac pair (with a mass above the TeV scale), one should have signals for both BR(µ − e − → e − e − , Al) and CR(µ − e, Al), observable at COMET's Phase II. This ISS mass region can be (complementary) probed via the search for cLFV Z → µτ decays, which are expected to be within reach of the future FCC-ee. Moreover, the lighter pseudo-Dirac pairs, belonging to a spectrum responsible for the large values of the cLFV observables, can also be directly searched for at facilities such as SHiP, DUNE and FCC-ee.
Sterile neutrinos could also contribute to Mu-Mu conversion, and to its cLFV decay Mu → e + e − . Although one could in principle have values of the effective coupling constant which would not be far from the most recent bounds, other cLFV observables, such as CR(µ − e, Au) and BR(µ → eee) preclude this possibility, and one is led to maximal values Re G MM ∼ 10 −13 in the "3+1" toy model and in the ISS, far from current experimental sensitivity. For completeness, we have also provided the expectations of these SM extensions regarding the Muonium decay. Although the experimental roadmap is not clear at present, it is possible that the cLFV Muonium decay will also be included in COMET's physics programme.
Concerning the νMSM, it was found that the predicted rates for the coherent µ − e conversion in Aluminium remain in general short of the expected future experimental sensitivity. Likewise, the expected contributions to the BR(µ − e − → e − e − , N) lie beyond experimental reach.
To summarise, our study has shown that minimal extensions of the SM with sterile neutrinos allow to explain future NP signals at high-intensity facilities dedicated to search for nuclei-assisted cLFV observables. In particular, that will be the case of CR(µ − e, N), BR(µ − e − → e − e − , N), the latter proving to be a very interesting probe of these SM extensions, in particular concerning the inverse seesaw model. On the contrary, the potential observation of Mu-Mu conversion or cLFV decay in the near future would disfavour sterile states (as those here considered) as the unique source of lepton flavour violation responsible for the latter cLFV observables.

A Appendix: Form factors for the cLFV decays
In this Appendix we provide the complete analytical expressions for the form factors and the loop functions entering the amplitudes of the cLFV processes described in Section 3.