Measurements of $\mu\to 3e$ Decay with Polarised Muons as a Probe of New Physics

Working within the Standard Model effective field theory approach, we examine the possibility to test charged lepton flavour violating (cLFV) new physics (NP) by angular measurements of the outgoing electron and positrons in the polarised $\mu^{+}\to e^{+}e^{+}e^{-}$ decay. This decay is planned to be studied with a sensitivity to the branching ratio of $B_{\mu\to3e}\sim 10^{-15}$ by the Mu3e experiment at PSI. To illustrate the potential of these measurements, we consider a set of eight effective operators generating the decay at the tree level, including dimension-five dipole operators and dimension-six scalar and vector four-fermion operators of different chiralities. We show that if the polarised $\mu^{+}\to e^{+}e^{+}e^{-}$ decay is observed and is induced by a single operator, data from three angular observables - two P-odd asymmetries and one T-odd asymmetry - can allow to discriminate between all considered operators with the exception of scalar and vector operators of opposite chirality. For these two types of operators, the P-odd asymmetry is maximal in magnitude; they can be distinguished, in principle, by measuring the helicities of the outgoing positrons. The observation of a non-zero T-odd asymmetry would indicate the presence in the decay rate of a dipole and vector operator interference term that is T (CP) violating.

processes have also been considered in models with supersymmetry, grand unification and continuous and discrete flavour symmetries [38][39][40][41][42][43][44][45][46]. Given the number and variety of possible NP models in which the cLFV processes are allowed and can proceed with rates close to the current upper bounds, one can instead study the impact of NP at low energies in the model-independent SM effective field theory. Below the scale of NP, Λ NP , but above the electroweak symmetry breaking (EWSB) scale, one adds higher-dimensional operators to the SM Lagrangian built from SM fields and respecting the SU(3) c × SU(2) L × U(1) Y gauge symmetry, where L (d) contains effective operators O i . An operator of dimension d is suppressed by (d − 4) powers of Λ NP . The only operator at dimension-five is the |∆L| = 2 Weinberg operator; additional (total) lepton number violating operators only arise at odd higher dimensions. Lepton flavour can be broken by operators of dimension-five and higher. At energies below the EWSB scale, one must match the operators above onto a set of dimensiond operators constructed from degrees of freedom lighter than M W which respect the broken-phase SU(3) c × U(1) em gauge symmetry. The Weinberg and higher-dimensional |∆L| = 2 operators, for example, generate Majorana neutrino masses, i.e. L (3) = − 1 2 M νν c L ν L + h.c.. Also generated in the low-energy effective field theory are operators which induce cLFV processes. For example, there are 42 (48) operators with three or four legs that trigger the µ ↔ e cLFV processes in Table 1 at tree (one-loop) level [58]. If any of these processes are observed, the question about which operator is generating the process will inevitably arise. The complementarity of µ ↔ e cLFV probes in identifying and constraining the possible operators was investigated in Ref. [58][59][60]. Taken into account was the renormalisation group (RG) running of the operator coefficients from the muon mass m µ up to the electroweak scale, M W . It was found that upper bounds on µ → eγ, µ → eγγ, µ → 3e decays and µ − e conversion in nuclei can probe orthogonal directions in the parameter space at a given scale. A higher degree of orthogonality can help discriminate a particular model, in which the indicated processes are triggered by specific combinations of operators whose coefficients at Λ NP are predicted by the model.
The quoted analyses only use upper limits on the branching ratios and capture rates of µ ↔ e flavour changing processes. If one or more of these processes is detected, additional information can be obtained from the angular distributions and helicities of the final state particles. For example, in Refs. [61], the full differential rate for µ + → e + e + e − decay with a polarised muon was used to construct three angular distribution asymmetries -two P-odd and one T-odd -sensitive to different combinations of operators (contributing to the process at tree level). The asymmetries were calculated in specific versions of SU(5) and SO(10) SUSY grand unified theories and it was found that the sign of the P-odd asymmetry can be different in the two models. The consideration of angular observables is particularly relevant as the next generation of µ ↔ e cLFV search experiments are likely to use a polarised muon beam. For example, the high-intensity muon beams (HIMB) proposal at PSI aims to provide a flux of O(10 9 ) longitudinallypolarised muons per second to the Mu3e and MEG II experiments [62], with an estimated average µ + polarisation of 93%.
In the present article, following a purely phenomenological approach, we investigate the possibility to discriminate between different effective operators contributing to the amplitude of the polarised µ + → e + e + e − decay by using data on the three angular observables -two P-odd and one T-odd asymmetries. Our aim is not to perform a comprehensive analysis of all operators that can induce the µ + → e + e + e − decay, but to illustrate the potential of the method on the example of a subset of eight commonly-used effective operators that generate the decay at tree level. Assuming the scale of NP Λ NP to be around 1 TeV, one can neglect to a good approximation in the analysed P-odd and T-odd asymmetries the effects of the renormalisation group running of the considered operators above and below the electroweak scale. In addition to the case of having just one operator responsible for the decay, we obtain predictions for the asymmetries of interest also in the cases when two operators are triggering the decay.
This work is organised as follows. In Section 2, we review the effective operators that are commonly used to parametrise the impact of NP above the EW scale on cLFV observables, such as the µ + → e + γ and µ + → e + e + e − decay rates, and that we are going to use also in our analysis. We then examine the differential branching ratio for the µ + → e + e + e − decay as a function of four kinematic variables, in particular of the two angles defining the directions of the outgoing positrons and electron. We next derive the total branching ratio and single differential branching ratios in the two angles. In doing so, we define three asymmetries (two P-odd and one T-odd) and derive their dependence on the coefficients of the effective operators of interest. We determine further the values of these asymmetries when only one operator, or a combination of two operators, induces the decay. In Section 3, we explore how a future experiment such as Mu3e can measure these asymmetries and if so, how the NP contributing to the cLFV process can be identified. We identify some unique signatures and comment on how helicity measurements of the outgoing positrons and electron can complement the angular observables. We conclude in Section 4.
2 New Physics Contributions to µ + → e + e + e − Decay Rate In the low-energy effective field theory (LEFT) of the SM, the eight operators of lowest dimension that contribute at tree level to the cLFV µ + → e + e + e − decay rate can be written, without loss of generality (up to Fierz rearrangements), as (see, e.g., [50,58,63]) Here, the dimension-five operators with coefficients C D,L and C D,R are dipole operators and the dimension-six are four-fermion operators. The coefficients have been normalised to the Fermi constant G F ; an O(1) coefficient therefore corresponds to Λ NP ∼ m W .
In the presence of the cLFV operators in Eq. (2), the differential branching ratio for µ + (p) → e + (k 1 )e + (k 2 )e − (k 3 ) decay with an incoming positive muon of polarisation ε = (0, P ), written in the muon rest frame, can be expressed as where dΩ ε = d cos θ ε dφ ε , P = | P | is the magnitude of the muon polarisation vector, and we define x 1 ≡ 2E k 1 /m µ and x 2 ≡ 2E k 2 /m µ , with E k 1 and E k 2 the energies of the outgoing positrons in the muon rest frame (the former taken to have larger outgoing energy, E k 1 ≥ E k 2 ). For convenience, the z-axis is set to be along the outgoing electron direction k 3 and the decay plane is fixed to be the x-z plane. The angle θ ε is defined to lie between P and the outgoing electron direction, while the angle φ ε is the azimuthal rotation in the x-y plane between P and the decay plane. In Eq. (3), the C i (x 1 , x 2 ) are functions of the coefficients in Eq. (2) and are given in Appendix A. These expressions are obtained in the limit m e → 0. The term in Eq. (6) proportional to cos θ ε is obtained from the dot product k 3 ·ε, while the sin θ ε cos φ ε term originates from the products k 1 · ε and k 2 · ε. On the other hand, the term proportional to sin θ ε sin φ ε is generated by ε µνρσ ε µ k ν 1 k ρ 2 k σ 3 , where ε µνρσ is the fourdimensional Levi-Civita symbol. Under a parity (P) transformation, the polarisation and momentum four-vectors transform as ε → ε = (0, P ) and k i → k i = (E k i , − k i ), respectively. Under a time (T) reversal, they instead transform as ε → ε = (0, − P ) and k i → k i = (E k i , − k i ). As a result, the dot products k i · ε and the terms proportional to cos θ ε and sin θ ε cos φ ε are P-odd, while ε µνρσ ε µ k ν 1 k ρ 2 k σ 3 and the term proportional to sin θ ε sin φ ε are T-odd. Assuming that NP contributing to the cLFV process is invariant under CPT, the non-conservation of T implies the violation of CP.
To obtain a double differential branching ratio in the angular variables cos θ ε and φ ε , it is possible to integrate Eq. (3) over the variables x 1 and x 2 . The kinematic limits of the three-body phase space in x 1 and x 2 can be found by first examining those in the kinematic variables s 1 = (p−k 1 ) 2 = m 2 µ (1−x 1 )+m 2 e and s 2 = (p−k 2 ) 2 = m 2 µ (1−x 2 )+m 2 e . These are where λ(x, y) = (1 − x − y) 2 − 4xy and the condition s 1 ≤ s 2 follows from Neglecting the electron mass, one obtains the integration ranges [0, 1] and [1 − x 1 , x 1 ] for x 1 and x 2 , respectively. However, given the choice can be taken as the final integration limits. These limits lead to logarithmic divergences since the functions C i (x 1 , x 2 ), more specifically, the terms in C i (x 1 , x 2 ) proportional to the dipole coefficients C D,L and C D,R , which correspond to photon-penguin diagrams, tend to infinity as x 1 → 1 (or equivalently, as s 1 → 0).
In Ref. [61] the collinear singularities were avoided by introducing the cut-off parameter δ in the x 1 integration as . For all operators in Eq. (3) other than the dipole operators, the δ → 0 limit can be taken without leading to a logarithmic divergence. Non-zero δ therefore gives small corrections at higher orders in an expansion in δ. On the other hand, the lowest order contribution for the dipole operators is proportional to the logarithm of δ. In their numerical estimates, the authors of Ref. [61] fixed δ = 0.02. To understand more exactly the dependence on m e , we instead perform the phase space integration for the dipole operators over the variables s 1 and s 2 , with the limits in Eqs. (4) and (5) and without neglecting m e . This allowed us to avoid introducing an undefined cut-off parameter in the calculations. Finally, we expand the final result in the small ratiom e ≡ me mµ . We obtain the double differential branching ratio as The total branching ratio, which can be found by integrating Eq. (6) over cos θ ε and φ ε with the ranges [−1, 1] and [0, 2π], respectively, can be written as where L 1 (x) = 17 + 12 ln 4x. Firstly, we note that the terms in the total branching ratio proportional to |C D,L | 2 and |C D,R | 2 are practically equal to previous results in the literature, cf. Refs [11][12][13]41]. Taking only the contribution from C D,L(R) , we obtain Table 2: Values of the asymmetries A P µ→3e and A P µ→3e when each of the coefficients of the operators in Eq. (2) are taken to be non-zero at a time, with the equivalent coefficients of Refs. [58,61] shown for clarity.
The differential branching ratio in cos θ ε can be found by integrating Eq. (6) over the azimuthal angle. This gives The P-odd asymmetry term A P µ→3e that remains is where L 2 (x) = 170 − 171 ln 2 + 12 ln 4x. In Table 2, we show the values of the asymmetry A P µ→3e when each of the coefficients in Eq. (2) is taken to be non-zero at a time. The operators with coefficients C ee S,LL , C ee S,RR , C ee V,LL and C ee V,RR induce maximal asymmetries, A P µ→3e = ±1, while those with coefficients C D,L , C D,R , C ee V,LR and C ee V,RL have |A P µ→3e | < 1. As is evident from Eq. (10), the value of A P µ→3e changes sign when the replacement L ↔ R is made. We note that the only operators that interfere at lowest order inm e , among those considered by us, are the dipole and vector operators. Interference between the other operators can therefore be neglected.
One can obtain the differential branching ratio in the azimuthal angle φ ε by integrating Eq. (6) over cos θ ε , giving with Here, the P-odd asymmetry A P µ→3e is given by and the T-odd asymmetry A T µ→3e by Fewer operators in Eq. (2) contribute to these asymmetries at leading order in m e /m µ .
In Table 2, we show the values of the asymmetry A P µ→3e when each of the coefficients in Eq. (2) is taken to be non-zero at a time. In isolation, the scalar and vector operators with coefficients C ee S,LL , C ee S,RR , C ee V,LL and C ee V,RR do not contribute to either A P µ→3e or A T µ→3e , while only the imaginary part of the interference between dipole and vector operators contributes to A T µ→3e . The general expressions for the three asymmetries are compatible with those derived in Ref. [61]. In contrast, our asymmetry A P µ→3e does not depend on the cut-off parameter δ which was set to δ = 0.02 in the numerical estimates in [61]. In our case, the role of this parameter is played by 4m 2 e /m 2 µ . In addition, the expressions for our asymmetries are a factor of two larger than those reported in Ref. [61] as a consequence of the way we define them. Further, our T-odd asymmetry is of opposite sign.
In order to observe CP violation in the µ + → e + e + e − decay, both the dipole and vector operators must be present; additionally, the coefficients of these operators must have a non-zero relative phase. This is possible, for example, in the canonical type I seesaw [12][13][14] and general SUSY [64,65] extensions. In the former, the dipole coefficient C D,L is generated by a photon-penguin diagram with a heavy neutral lepton and W ± boson in the loop. The vector coefficient C ee V,LL , on the other hand, gets contributions from box diagrams containing two heavy neutral leptons and W ± bosons. The penguin and box diagram amplitudes are proportional to U ei U * µj and U ei U * µj U ei U * ej , respectively, where U i is the mixing between the charged lepton and neutrino mass eigenstate N i (i = 1, · · · , 3 + n, with n the number of heavy neutral leptons). The dependence of the coefficients C D,L and C ee V,LL on additional Dirac CP phases in the extended mixing matrix can therefore be different. In SUSY extensions, the coefficients can also depend differently on phases in the soft SUSY-breaking sector.

Phenomenological Analysis
We now discuss the implications of a detection of the µ + → e + e + e − decay of a polarised µ + . We consider, in particular, how angular and helicity measurements of the outgoing positrons and electron can help to identify the NP contributing to the cLFV process.
In future, the µ + → e + e + e − decay will be searched for by the Mu3e experiment at PSI. Phase I of the experiment will receive O(10 8 ) positive muons per second with average momenta p ∼ 29.8 MeV from pion decays at rest near the surface of the production target (the πE5 compact muon beamline). The muons produced in this way are fully polarised; the beam is thus estimated to have an average polarisation of P ∼ 93%. This phase of the experiment is expected to probe down to B µ→3e ∼ 2 × 10 −15 , in contrast to the current limit of B µ→3e < 1.0 × 10 −12 from the SINDRUM experiment [49]. This is a much larger improvement on the current experimental bound compared to that expected with the planned search for the µ + → e + γ decay. Using the πE5 beamline, the MEG II experiment aims to reach B µ→eγ ∼ 6 × 10 −14 compared the current bound of B µ→eγ < 4.2 × 10 −13 from the original MEG experiment. In the longer term, the HIMB proposal aims to provide O(10 9 ) muons per second to Phase II of the Mu3e experiment, allowing to reach even higher sensitivity to B µ→3e . The Mu3e detector is designed to optimally track the positrons and electron from decays of muons at rest, and therefore to fully reconstruct the event kinematics. To cover a large portion of the final state phase space, the experiment aims to reconstruct tracks over a large solid angle (limited by the beam entry and exit points) with momenta ranging from ∼ 1 MeV to half the muon mass [50]. The experiment must be able to discriminate the µ + → e + e + e − decay from backgrounds such as µ + → e + e + e − νν and a Michel positron accompanied by a e + e − pair from photon conversion or Bhabha scattering. The experiment can nevertheless exploit the kinematics of the µ + → e + e + e − decay at rest, | i k i | = 0 and i E k i = m µ , to mitigate these backgrounds.
To understand how the asymmetries can be measured by an experiment such as Mu3e (in the scenario where µ + → e + e + e − is observed), we note that they can be extracted from Eq. (6) as follows: The P-odd asymmetry A P µ→3e can be determined by counting the number of electrons with outgoing angles θ ε in the range [0, π 2 ] minus those in the range [ π 2 , π]. Other than the outgoing electron direction, this only requires knowledge of the magnitude and direction of the muon polarisation P . To instead determine the P-odd and T-odd asymmetries A P µ→3e and A T µ→3e , it is also necessary to measure the directions of the outgoing positrons to identify the decay plane and the azimuthal angle φ ε . The asymmetry A P µ→3e can be ascertained by counting the number of decay topologies with angles φ ε in the range [− π 2 , π 2 ] minus those in the range [ π 2 , 3π 2 ]. Conversely, for the asymmetry A T µ→3e one counts the number of decays with angles φ ε in the range [0, π] minus those in the range [−π, 0].
In Table 2 we have already shown the values of the asymmetries A P µ→3e and A P µ→3e when each effective coefficient is taken to be non-zero at a time. It follows from the results in Table 2 that the measurement of A P µ→3e can allow to discriminate between the dipole, scalar RR (scalar LL), vector LL (vector RR), vector LR and vector RL operators, while it will be impossible to distinguish between scalar RR (scalar LL) and vector RR (vector LL) operators. With the measurement of A P µ→3e alone it may be possible to discriminate between the dipole, scalar or vector LL or RR, and vector LR or RL operators. At the same time, the difference between the values of A P µ→3e due to the dipole R (dipole L) and vector RL (vector LR) operators is very small and thus distinguishing between these operators requires an extremely precise measurement of A P µ→3e . If more than one operator is present, the asymmetries instead take a range of values depending on the relative sizes of the respective coefficients. For example, in the presence of the vector coefficients C ee V,LL and C ee V,LR (with all other coefficients set to zero), we have which clearly range from the minimum values A P− µ→3e = −1 and A P − µ→3e = 0 when x 1 and the maximum values A P+ µ→3e = 1 3 and A P + µ→3e = 32 105 when x 1. Actually, the values of the asymmetries A P µ→3e and A P µ→3e are correlated: A P µ→3e = 8 35 (1 + A P µ→3e ). Observing this correlation would strongly imply that the combination of the vector LL and LR operators triggers the µ + → e + e + e − decay.
A measurement of the maximal asymmetry A P µ→3e = 1 (−1) would imply that one of the coefficients C ee S,LL or C ee V,RR (C ee S,RR or C ee V,LL ) dominates over the others. The other P-odd asymmetry A P µ→3e is approximately zero when both C ee S,LL and C ee V,RR dominate, while the T-odd asymmetry (which relies on the interference between the dipole and vector operators) is also suppressed; neither can therefore be used to distinguish the two contributions. However, the scalar and vector operators can be distinguished with an additional measurement of the helicity of the outgoing positrons. In Table 3, we show the helicities of the final states for each of the operators in Eq. (2). For both the scalar and vector operators (with the coefficients C ee S,LL and C ee V,RR , respectively) the helicity of the outgoing electron is positive. The outgoing positrons, on the other hand, have positive helicities for the scalar operator and negative helicities for the vector operator. A measurement of these helicities can therefore be used to distinguish the two scenarios.
In the presence of the dipole and vector operators, which interfere at leading order in m e , the asymmetry does not simply range between the values in Table 2 Figure 1: Ranges of the asymmetries A P µ→3e (blue), A P µ→3e (green), and A T µ→3e (red), with solid lines indicating the maximum asymmetries A X+ µ→3e and the dashed lines the minimum asymmetries A X− µ→3e . In the left panel, only the coefficients C D,L and C ee V,LL are taken to be non-zero, while in the right panel only C D,L and C ee V,LR are considered. The ranges are shown as a function of the phase difference ∆φ between the coefficients. on the relative phase ∆φ between the coefficients. For the coefficients C D,L and C ee V,LL , the minimum and maximum values of A P µ→3e depend on the phase as .
In the left panel of Fig. 1, we depict the ranges of the asymmetries A P µ→3e (blue), A P µ→3e (green), and A T µ→3e (red) as a function of the relative phase ∆φ when only C D,L and C ee V,LL are taken to be non-zero. The solid and dashed lines indicate the maximum and minimum asymmetries A X+ µ→3e and A X− µ→3e , respectively. It can be seen that for ∆φ = 0, the P-odd asymmetry has the wider range of −1 ≤ A P µ→3e ≤ 0.87 compared to −1 ≤ A P µ→3e ≤ 0.63 for ∆φ = π 2 . In Fig. 1, right panel, we also show the ranges when C D,L and C ee V,LR are taken to be non-zero.
If a future experiment does not observe µ + → e + e + e − decay, this sets an upper limit on the branching ratio B µ→eγ in Eq. (7) and therefore on each of the coefficients in Eq. (2) if the other coefficients are set to zero. A positive detection of µ + → e + e + e − decay instead sets an allowed range for each coefficient. For example, an observed branching Figure 2: Assuming that the cLFV decay µ + → e + e + e − is observed with a branching ratio of B µ→3e = (2.0 ± 0.5) × 10 −15 , we depict the allowed regions in the (C D,L , C ee V,LL ) parameter space for a relative phase difference between the coefficients of ∆φ = 0 (top) and ∆φ = π 2 (bottom). We also show the size of the asymmetries A P µ→3e (left), A P µ→3e (middle), and A T µ→3e (right). Measurements of these asymmetries can further constrain the parameter space, or indicate additional non-zero effective coefficients. ratio of B µ→3e = (2.0 ± 0.5) × 10 −15 constrains the dipole coefficients to lie in the range 8.1 × 10 −9 ≤ |C D,L(R) | ≤ 1.0 × 10 −8 . If two coefficients are allowed to be non-zero, the favoured region becomes an circle or ellipse in the parameter space of the coefficients. If the dipole and vector coefficients C D,L and C ee V,LL (or C ee V,LR ) are taken to be non-zero, the alignment of the ellipse is controlled by the phase difference ∆φ between the coefficients. In Fig. 2, we plot the allowed regions in the (C D,L , C ee V,LL ) plane if a future experiment measures B µ→3e = (2.0 ± 0.5) × 10 −15 , for ∆φ = 0 (top panels) and ∆φ = π 2 (bottom panels). For ∆φ = 0, the real part of C D,L C ee * V,LL in Eq. (7) is non-zero and therefore the ellipse is not aligned along the C D,L and C ee V,LL axes. For ∆φ = π 2 , the real part is indeed zero and the ellipse is aligned. In Fig. 2, we allow for negative values of the coefficients C D,L and C ee V,LL . The allowed region for the phase difference ∆φ = π is therefore the mirror of the ∆φ = 0 ellipse along either the C D,L or C ee V,LL axis. In the left, centre and right panels of Fig. 2 we also show the magnitudes of the asymmetries A P µ→3e (blue), A P µ→3e (green), and A T µ→3e (red), respectively. The ranges of these asymmetries are consistent with the minimum and maximum values in Eqs. (19), (20) and (21).
We would finally like to note that a measurement of the process µ − e − → e − e − in a muonic atom, first suggested in Ref. [66], can provide additional information on the effective coefficients of Eq. (2). Assuming a non-relativistic bound state for the muon, a 1S electron, a point nuclear charge and the plane wave approximation for the final state electrons, one obtains an analytical result for the rate Γ(µ − e − → e − e − ) with interference terms between the four-fermion operators [67]. The dependence of the rate on the dipole coefficients was studied separately in Ref. [68]; however, in principle there are also interference terms between the dipole and four-fermion operators. In a similar manner to the µ + → e + e + e − decay, the angular distribution of the outgoing electrons for a polarised muon can further identify the contribution of each operator [69]. In future, it is possible for Phase I of the COMET experiment to search for the clean signal of two back-to-back electrons with a total energy of the muon mass minus the muon binding energy [53]. A complete study of the complementarity of the cLFV processes in Table 1 (including measurements of angular distributions and helicities) is beyond the scope of this work.

Conclusions
In the present article, following a purely phenomenological approach, we have investigated the possibility to discriminate between different effective operators contributing at tree level to the amplitude of the polarised µ + → e + e + e − decay by using data on the three angular observables -two P-odd and one T-odd asymmetries, A P µ→3e , A P µ→3e and A T µ→3e , defined in Eqs. (15), (16) and (17). The consideration of these angular observables is particularly relevant as the next generation of experiments on µ ↔ e cLFV processes are likely to use polarised muon beams. In particular, the high-intensity muon beams (HIMB) proposal at PSI aims to provide a flux of O(10 9 ) longitudinally-polarised muons per second to the Mu3e and MEG II experiments [62] on µ + → e + e + e − and µ + → e + γ decays, with the average µ + polarisation estimated to be approximately of 93%. The measurement of the P-odd asymmetry A P µ→3e in the decay µ + → e + e + e − requires the knowledge only of the direction of the outgoing electron and of the magnitude and the direction of the µ + polarisation. To determine the P-odd and T-odd asymmetries A P µ→3e and A T µ→3e , it is also necessary to measure the directions of the outgoing positrons to identify the decay plane and the azimuthal angle φ ε between the µ + polarisation vector and the decay plane. The effective operators considered by us respect the CPT symmetry and thus an observation of a non-zero T-odd asymmetry would imply violation of the CP symmetry. This is only possible if at least two operators, which have different phases and whose interference term in the decay rate is not zero (i.e., is not strongly suppressed), generate the µ + → e + e + e − decay.
The aim of our study was not to perform a comprehensive analysis of all operators that can induce the µ + → e + e + e − decay, but to illustrate the effectiveness of the method on the example of the subset of eight commonly-used effective operators that generate the decay at tree level, given in Eq. (2). This subset includes dimension-five left-handed (L) and right-handed (R) dipole operators and dimension-six four-fermion operators involving products of scalar and vector currents of different chiralities (left-left (LL), right-right (RR), left-right (LR) and right-left (RL)). The only operators, among those considered by us, that were found to have significant interference terms in the µ + → e + e + e − decay rate at lowest order in m 2 e /m 2 µ , are the dipole and vector operators. The interference between the other operators was therefore neglected. Assuming the scale of NP to be around Λ NP ∼ 1 TeV, one can neglect to a good approximation in the analysed P-odd and T-odd asymmetries the effects of the renormalisation group running of the considered operators above and below the electroweak scale.
It follows from the results of our study that if only one of the set of considered operators triggers the µ + → e + e + e − decay, the measurement of the P-odd asymmetry A P µ→3e can allow to discriminate between the dipole, scalar RR (scalar LL), vector LL (vector RR), vector LR and vector RL operators, while it will be impossible to distinguish the scalar RR (scalar LL) and vector RR (vector LL) operators (cf. Table 2). With the measurement of A P µ→3e alone it may be possible to discriminate between the dipole, scalar or vector LL or RR, and vector LR or RL operators. At the same time, the difference between the values of A P µ→3e due to the dipole R (dipole L) and vector RL (vector LR) operators is very small ( Table 2) and thus distinguishing between these operators requires an extremely high precision in the measurement of A P µ→3e . In addition to the case of having just one operator responsible for the decay, we have derived predictions for the asymmetries of interest also in the cases when two operators are triggering the decay (Section 3). Some of these predictions are illustrated in Figs. 1 and 2. The observation of a non-zero T-odd asymmetry in the polarised µ + → e + e + e − decay, for example, would imply that, within the considered set, dipole and vector operators with different phases are generating the decay.
Our study has demonstrated that if the µ ↔ e cLFV processes are observed, measuring the angular distributions and the corresponding P-odd and T-odd asymmetries of the µ + → e + e + e − decay with a polarised µ + would provide an additional effective method of identifying the beyond the SM physics operators that trigger these processes.

A Relevant Functions
The functions C i (x 1 , x 2 ) appearing in the differential branching ratio for the µ + → e + e + e − decay in Eq. (3), where x 1 ≡ 2E k 1 /m µ and x 2 ≡ 2E k 2 /m µ and E k 1 ≥ E k 1 are the outgoing positron energies, are of the form Here, the functions F i (x 1 , x 2 ) are given by the functions G i (x 1 , x 2 ) by and finally the H i (x 1 , x 2 ) functions by , , In these expressions we have neglected the electron mass m e . These expressions are in agreement with those in Ref. [61]. It is straightforward to find the single differential branching ratio in the parameter where the F i (x 1 ) functions are computed as As mentioned in the main text, for the dipole operators one encounters a logarithmic divergence when neglecting the electron mass and performing the phase-space integral over x 1 . The approach of Ref. [61] was to introduce a cut-off parameter δ in the upper integration limit of x 1 . However, the physical interpretation of δ is unclear, and for numerical estimates it must be set arbitrarily to some value. In order to find the exact dependence of the total branching ratio B µ→3e and the P-odd asymmetry A P µ→3e on the electron mass m e , we instead express the differential branching ratio for µ + → e + e + e − decay as a function of the kinematic variables s 1 = (p − k 1 ) 2 = m 2 µ (1 − x 1 ) + m 2 e and s 2 = (p − k 2 ) 2 = m 2 µ (1 − x 2 ) + m 2 e . We find (considering only the dipole operators) dB D,L(R) µ→3e ds 1 ds 2 dΩ ε = 3 2π C 1 (s 1 , s 2 ) + C 2 (s 1 , s 2 )P cos θ ε + C 3 (s 1 , s 2 )P sin θ ε cos φ ε , where the functions C i (s 1 , s 2 ) are given by with X ≡ |eC D,L | 2 + |eC D,R | 2 and Y ≡ 1 where we have used λ(x, y) = (1 − x − y) 2 − 4xy and s = m 2 µ + 3m 2 e − s 1 − s 2 . Using the phase-space integration limits of s 1 and s 2 in Eqs. (4) and (5), respectively, we then obtain the total branching ratio as where L 1 (x) = 17 + 12 ln 4x and we have only retained the terms at lowest order in an expansion inm e = me mµ . Similarly, the P-odd asymmetry is found to be A P where L 2 (x) = 170 − 171 ln 2 + 12 ln 4x.