LFV Higgs and Z-boson decays: leptonic CPV phases and CP asymmetries

Heavy neutral leptons are motivated by several extensions of the Standard Model and their presence induces modifications in the lepton mixing matrix, including new Dirac and Majorana CP violating phases. It has been recently shown that these phases play an important role in lepton number and lepton flavour violating decays and transitions, with a striking impact for the predicted rates of certain observables. In this work, we now consider the potential role of the heavy neutral fermions and of the new CP violating phases in Higgs and Z-boson lepton flavour violating decays H,Z→ℓα±ℓβ∓\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H,Z\rightarrow \ell _\alpha ^\pm \ell _\beta ^\mp $$\end{document}, as well as in the corresponding CP asymmetries. In order to allow for a thorough analysis, we derive the full analytical expressions of the latter observables, taking into account the effect of the (final state) charged lepton masses. A comprehensive exploration of lepton flavour violating Z and Higgs decays confirms that these are very sensitive to the presence of the additional heavy neutral leptons. Moreover, the new CPV phases have a clear impact on the decay rates, leading to both destructive and constructive interferences. Interestingly, in the μτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu \tau $$\end{document} sector, the Z→ℓα±ℓβ∓\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z\rightarrow \ell _\alpha ^\pm \ell _\beta ^\mp $$\end{document} rates are within reach of FCC-ee, with associated CP asymmetries that can potentially reach up to 20–30%.


Introduction
Numerous New Physics (NP) models have been explored in recent years in order to explain and understand the phenomenology of neutrino oscillations.Extensions of the Standard Model (SM) accounting for neutrino data open the door to new leptonic mixings -possibly going beyond the minimal Pontecorvo-Maki-Nagakawa-Sakata (PMNS) left-handed leptonic mixing matrix U PMNS -as well as to new sources of CP violation, through the presence of Dirac (and possibly Majorana) phases.Notice that unless CP conservation is enforced, CP violating (CPV) phases are in general associated with multi-generational fermion mixings.Leptonic sources of CP violation are a key ingredient to constructions aiming at explaining the observed baryon asymmetry of the Universe via leptogenesis; moreover, they have been found to have a significant impact on a number of observables, which are currently being investigated at high-energy colliders or in dedicated high-intensity experiments.
Among the many NP constructions successfully addressing the problem of neutrino mass generation and the observed lepton mixing pattern, one of the most minimal (yet very appealing) possibilities is that of SM extensions via heavy neutral leptons (HNL), in general of Majorana nature as is the case of right-handed neutrinos.Present in numerous scenarios (such as the type I seesaw and its variants, or arising from complete constructions as SO (10) grand unified models, among other examples), sterile fermions can be at the source of extensive contributions to numerous phenomena beyond neutrino oscillations: depending on the mass of these new states (and their interactions with the active neutrinos), such phenomena range from flavour observables, to high-energy collider signals, with a non-negligible impact on astroparticle physics or even cosmology.
The new CPV phases (Dirac and/or Majorana) are also expected to lead to significant (and detectable) effects in several searches.In addition to CP-violating observables (as is the case of electric dipole moments of charged leptons [33][34][35][36]), the new phases are also at the source of interference effects in LNV (and cLFV) semileptonic meson and tau decays [37] -especially in the case in which the SM is extended by two (or more) heavy neutrino states.In particular, the role of a second HNL has been explored in studies concerning the effects of CP violation, be it in high-scale seesaw scenarios, in the context of renormalisation group running and thermal leptogenesis [38,39], in assessing the possibility of resonant CP violation [40], or upon considering possible effects in forward-backward asymmetries at an electron-positron collider [41].Other studies have compared the expected number of events associated with same-sign and opposite-sign dileptons at colliders in the framework of leftright symmetric models [41][42][43].
Recently, a study was carried out, dedicated to investigating the impact of leptonic CPV phases on numerous cLFV transitions and decays [44]; working in the framework of a minimal low-energy extension of the SM via 2 Majorana heavy sterile states, the analysis revealed that the new phases could have a significant impact for the predicted cLFV rates, leading to a possible loss of correlation between cLFV observables (which would be otherwise manifest in the CP conserving case).Building upon this first thorough analysis of the effects of CPV phases regarding leptonic cLFV transitions and decays, we now address in detail the prospects in what concerns SM neutral boson decays.The present work is thus dedicated to a detailed investigation of the role of the CPV phases in cLFV Higgs and Z-boson decays (H, Z → ± α ∓ β ).A first part of the analysis concerns an in-depth study of the decay rates, revisiting the full computation of the amplitudes.In particular, we derive the full effective cLFV Z and Higgs vertices in the presence of heavy neutral leptons without any simplifying assumptions, finding good agreement with previous calculations [17,19,28,32,[45][46][47][48][49][50][51][52].
As argued in [44], the possible effects of CPV phases on the cLFV observables can lead to a strong loss of correlation between them, and thus smear out otherwise clean flavour patterns.A similar situation can be encountered in leptoquark models (albeit not necessarily due to CP violating interactions), such that an identification of specific leptoquark states would be often unfeasible just by taking into account future data on cLFV (meson) decays.As suggested in the literature, these difficulties can be overcome (or at least reduced) by comparing different CP final states, rather than taking the average of the observables (see e.g.[53][54][55]).In specific cLFV meson decays (e.g.B → Kµ ± τ ∓ ), measurements of both possible final state charge assignments have been carried out [56,57], thus allowing to constrain different flavour structures (as for example in the case of certain leptoquark models).
Similarly, in heavy neutral lepton models, it can be difficult to identify the presence of new CPV phases at the origin of peculiar flavour patterns.Henceforth, it could be interesting to compare the different CP final states of cLFV transitions: as a (more) direct manifestation of leptonic CP violation, one can then consider CP asymmetries in cLFV Z (and Higgs) decays.The latter observables could offer further hints not only on the underlying source of flavour violation, but also on its CP structure (possibly in connection with leptogenesis [58]).We thus carefully examine the impact of Dirac and Majorana CPV phases on such observables, and the results of our phenomenological analysis suggest that the CP asymmetries can reach up to 20 − 30%, in particular for cLFV Z → µτ decay rates, which can be within the reach of FCC-ee [59].
As we argue throughout this study, CP violating leptonic phases clearly have extensive effects on cLFV transitions and decays, and their presence must be systematically taken into account to fully explore HNL extensions of the SM.
The manuscript is organised as follows: we begin by briefly describing the features of this simple SM extension in Section 2, after which we revisit cLFV Z and Higgs boson decays, presenting full analytical expressions in Section 3. Our results concerning the role of the Dirac and Majorana CPV phases are collected in Section 4, which is followed by a dedicated discussion of CP asymmetries in cLFV neutral boson decays (Section 5).Finally, we summarise our main findings in the Conclusions.Further information on the analytical derivation and analysis is collected in the Appendices.

A minimal HNL extension of the SM
As done in a previous analysis [44], we follow a phenomenological approach, relying on a minimal SM extension which allows studying the effects of CPV phase-induced interferences in the lepton sector.We thus consider the SM extended by 2 sterile Majorana states; the neutral lepton sector is thus composed of 5 massive states (with masses m i , i = 1, ..., 5), which include 3 light (mostly active) neutrinos and two heavier states.From a low-energy perspective -and independently of the mechanism at the origin of neutrino mass generation -, the leptonic mixings are impacted by the presence of the new heavy states.Instead of the 3 × 3 U PMNS mixing matrix, a 5 × 5 unitary matrix U now parametrises mixings in the lepton sector; notice that the left-handed leptonic mixings are now encoded in its 3 × 3 upper left block (the would-be PMNS matrix, ŨPMNS ); as a consequence, the mixing in charged current interactions will be parametrised via a rectangular 3 × 5mixing matrix.The deviations of ŨPMNS from unitarity [21,60] can be conveniently cast in terms of a matrix η as follows [61] U In agreement with what was done in [44], U can be parametrised through a series of ten (complex) rotations R ij (with i = j), and a diagonal matrix including the four physical Majorana phases, ϕ i [33] ordering of the light neutrino spectrum; notice that the predictions for cLFV transitions are mostly independent of the ordering of the spectrum).These are collected in Appendix A, where we also discuss other relevant constraints on SM extensions with sterile fermions.
As already stated, due to the presence of the additional sterile fermion states, the would-be PMNS matrix is no longer unitary, thus leading to modified charged and neutral lepton currents.In the physical basis, and for the generic case of a number n S of sterile fermions, the relevant terms of the interaction Lagrangian can be cast as in which α, ρ = 1, . . ., 3 denote the flavour of the charged leptons, while i, j = 1, . . ., 3 + n S correspond to the physical (massive) neutrino states; P L,R = (1 ∓ γ 5 )/2, g w is the weak coupling constant, and Regarding the interaction of Z bosons with charged leptons, the SM vector and axial-vector currents can be cast in terms of the C V and C A coefficients, respectively given by C V = − 1 2 + 2 sin 2 θ w and C A = − 1 2 .Moreover, the coefficients C ij are defined as: For completeness, the new interaction vertices (Feynman rules) are presented in Appendix B.

Revisiting cLFV neutral boson decays: analytic expressions
In this section we present expressions for the relevant cLFV decays of Z and Higgs bosons in the presence of n S singlet fermions.In general, the cLFV boson decays receive contributions at the oneloop level from boson vertex corrections, as well as from charged lepton self-energy diagrams.In the following, no approximation will be made leading to the full amplitudes.In this section, the results are obtained working in unitary gauge; for completeness, we also present the full expressions derived in Feynman-'t Hooft gauge in Appendix C.

Z-boson decays
For on-shell charged leptons, the cLFV Z vertex can in general be decomposed as in which q µ = p µ α + p µ β is the momentum of the Z boson and σ µν = i 2 [γ µ , γ ν ]; F L,R S,V,T are scalar, vector and tensor form factors.For on-shell Z-boson decays, the scalar amplitudes do not contribute due to the Ward identity q µ * µ (q) = 0, but they can nevertheless appear if the Z boson is an offshell intermediate state.Furthermore, all amplitudes but F L V vanish in the limit of massless charged leptons.We emphasise that in the calculation of the contributing diagrams, which are shown in Fig. 1, we do not neglect the final state charged lepton masses.Although the assumption of vanishing charged lepton masses is indeed justified for the cLFV decay rates, as we will discuss later on, the behaviour of the CP asymmetries significantly depends on the charged lepton masses, thus demanding a full computation 1 .
We separate the contributions to the invariant amplitudes according to the different topologies (the superscript refers to the diagrams of Fig. 1) as: The different contributions are given by 2 with the Passarino-Veltman functions [63] following the LoopTools convention [64]: , and W , m 2 i , m 2 j ) (with rs = 0, 1, 2, 11, 22, 12, 00); we integrate in D = 4 − 2ε dimensions.
Note that in all above contributions, the terms proportional to A 0 (M 2 W ) and to B 0,00 (q 2 , M 2 W , M 2 W ) do not contribute to the amplitude since they do not depend on the internal neutrino masses (and therefore cancel due to the (semi-)unitarity3 of U).Among the distinct contributions listed above, terms proportional to the B 0,1 and C 00 functions contain ultraviolet (UV) divergences which are regulated dimensionally (we integrate in (4 − 2ε) dimensions).In the case of on-shell Z-boson decays (q 2 = M 2 Z ), the sum over all contributions is, however, manifestly finite as shown in Appendix C. We thus set q 2 = M 2 Z and safely take the limit D → 4. A few additional comments are in order: firstly, we note here that all contributions, except F L V , are generically suppressed by at least one power of the charged lepton masses, and can be therefore safely neglected in what regards their contribution to the decay rate.Secondly, note that functions of the form B 0,1 (m 2 α,β , M 2 W , m 2 i ) (or similar) contain potentially problematic large logarithms ∼ log(m 2 i /M 2 W ) be it in the limit of very small or very large neutrino masses.However, in the sum over all contributions (and internal states) these logarithms cancel and, as expected, the light neutrino contribution to the amplitude is negligible.Finally, we have numerically checked that our results agree with those of [17,45,52], which were obtained in the limit of vanishing charged lepton masses.Neglecting suppressed contributions, the decay rate can be written as4 We notice that concerning the evaluation of the Passarino-Veltman functions we use the public Fortran code LoopTools [64] (wrapped into a dedicated python code).From a numerical point of view, the observable evaluated is the CP averaged decay rate as current search results do not distinguish the charges of the final state leptons.

Higgs decays:
In the presence of n S sterile fermions, one finds new contributions leading to cLFV H → ∓ α ± β (α = β).Working in unitary gauge, one can identify three types of diagrams as depicted in Fig. 2: corrections to the H vertex due to two neutrinos (cf.diagram (a)), two W ± bosons and one neutrino (diagram (b)), as well as corrections to the fermion legs (diagrams (c) and (d)).
We have carried out the evaluation of the new contributions, leading to the following expression for the branching ratio: where Γ H denotes the total Higgs width.The cLFV width can be cast as with M H the SM Higgs mass, and as illustrated in Fig. 2. Below, we collect the expressions for each of the form factors, F i L(R) , cast in terms of the relevant entries of the leptonic mixing matrix U αi , the masses of the external leptons, m α(β) , and finally of the physical masses of the neutral fermions running in the loop, m i(j) (i, j = 1, ..., 3 + n S ).As before, our results are given in terms of Passarino-Veltman functions following LoopTools conventions; for previous works -relying on different Passarino-Veltman basis -, see [32,48,49,51].
with the Passarino-Veltman functions where The following section will be devoted to numerically studying the rates for the cLFV Higgs and Z decays, in particular emphasising the role of the leptonic CPV Dirac and Majorana phases.

cLFV neutral boson decays: the impact of leptonic CP phases
We now proceed to address how the presence of the Dirac and Majorana CP violating phases can modify the predicted rates for cLFV Z and Higgs decays.We first discuss "simplified scenarios", without aiming at a complete phenomenological discussion; a thorough phenomenological study will be carried out in subsequent sections.

Simplified scenarios
We begin by considering some illustrative scenarios, whose aim is to offer a first insight into the role of phases in neutral boson decays, and a direct comparison of the CP conserving limit with the CPV case (non-vanishing Dirac and/or Majorana phases).Notice that at this stage we do not apply any constraints, so that one might encounter scenarios which are already experimentally excluded.This first exploratory study relies on simple assumptions: we set θ α4 = θ α5 and, concerning the new sterile states, we assume that they are heavy and mass-degenerate, above the electroweak (EW) scale, typically of the order of a few TeV, m 4 = m 5 Λ EW .In particular, we consider a simple benchmark choice for the values of the mixing angles, sin θ 14 = sin θ 15 = 10 −3 , sin θ 24 = sin θ 25 = 0.01 and sin θ 34 = sin θ 35 = 0.1, as well as three representative values of the (degenerate) heavy fermion masses m 4 = m 5 = 1, 5, 10 TeV.In Fig. 3 we display the cLFV rates for the B → α β decay as a function of the degenerate masses of the heavy sterile states (m 4,5 ), in the CP conserving limit with all CPV phases set to zero.Notice that in this figure (and generically in what follows), the observable associated with the cLFV decays of a B-boson refers to the average of the two combinations of final state charges, with B = H, Z.In our analysis, we further take the following values for the neutral boson decay widths, Γ Z = 2.4952 GeV and Γ H = 4.115 MeV , with Γ Z as given in [65], and the SM prediction for Γ H from the CERN Yellow Report [66] (for As seen from Fig. 3, the sterile states lead to significant contributions to the cLFV decay rates, both for the Z and the Higgs; the dependence of the rates on the heavy states' mass confirms the findings of previous studies (see, e.g.[17,32,51]) 5 .In particular, and although in the present study no approximation is made regarding the mass of the final state leptons, the predictions for the cLFV Z decay rates are in excellent agreement with the results of [17,45,52], obtained in the limit in which m α,β → 0.
Although the decays are not formally identical (scalar vs. vector boson decays), the common topology of the one-loop contributions leads to similar behaviours for the observables.In both cases, and in the limit of large m 4,5 , the diagrams with two virtual neutrinos (diagrams (a) of Figs. 1 and 2), provide the dominant contribution.In Fig. 4 we now display the (individual) effects of Dirac and Majorana CPV phases for the distinct cLFV rates H → α β .As mentioned before, we also consider three regimes for the degenerate heavy fermion masses, m 4 = m 5 = 1, 5, 10 TeV.While the effects of the Dirac phases δ 14 (δ 24 ) are absent for µτ (eτ ) final states -a consequence of the very simplifying hypotheses for the mixing anglesall cLFV Higgs decays exhibit a clear dependence on δ 34 .This is a consequence of the contributions proportional to C ij m 2 i , which introduce s 34,35 -enhanced terms, see Eqs. (23,24).In all cases, Dirac phases open the possibility of a strong suppression in the decay rates for values ∼ π (albeit for different cLFV decays, this behaviour had already been identified in [44]).The cLFV Higgs decay rates are also sensitive to Majorana phases, which contribute via the terms proportional to C * ij , see Eqs. (23,24).Although not displayed here, the behaviour of the BR(H → α β ) with respect to other phases (δ i5 and ϕ 5 ) is identical to what is illustrated in Fig. 4.
As already noticed in [44], cLFV Z decays exhibit a clear dependence on both Dirac and Majorana CP violating phases.For completeness we display in Fig. 5 a study analogous to that of Fig. 4, considering the impact of Dirac and Majorana phases on BR(Z → α β ), for distinct masses of the heavy mediators.Moreover, these decays are closely related to several other cLFV leptonic decays, in particular to processes receiving contributions from Z-penguin topologies.While in the present study -and as already mentioned -no approximations were made leading to the computation of the cLFV Z decay rates, the findings are in very good agreement with the previous study of [44].Finally, notice that for both cases of Z and Higgs cLFV decays the effects of the phases are amplified for increasing masses of the sterile states.
In order to illustrate the joint effect of Dirac and Majorana phases, we display in Fig. 6 the iso-surfaces for BR(H → τ µ) as spanned by the CPV phases δ 34 and ϕ 4 , for fixed values of the mixing angles and sterile states' masses.As visible, different regimes for the CPV phases can account for variations on the predicted values for the rates by as much as 7 orders of magnitude, due to the impact of the (destructive) interference effects.The comparison of left and right panels also reveals the possibility of constructive interferences: in this case, effects already emerge for opposite-sign mixings (illustrated by s 24 = −s 25 ), leading to a significant increase in the branching ratio.

Phenomenological impact
In this section we complete and generalise the discussion of the previous one, carrying out a comprehensive phenomenological analysis to fully assess the effect of the CPV Dirac and Majorana phases on Z-and Higgs boson cLFV decays, in the framework of a minimal low-energy extension of the SM via 2 heavy neutral leptons.
Figure 4: Rates of cLFV Higgs decays as a function of the CP violating Dirac and Majorana phases.
From left to right and top to bottom, dependence on δ 14 , δ 24 , δ 34 and ϕ 4 (all phases set to zero in each case, except for the one displayed).The colour code denotes the flavour composition of the final state lepton pair: eµ (blue), eτ (orange) and µτ (green).In both panels, solid, dashed and dotted lines respectively correspond to the following heavy neutral fermion masses: m 4 = m 5 = 1, 5, 10 TeV.
We thus now take into account all available experimental constraints which include, among others, limits on the active-sterile mixings, negative results of direct and indirect searches for the sterile states, as well as EW precision tests; likewise, bounds on searches for other cLFV transitions -such as radiative and three-body leptonic decays as well as rare transitions in muonic atoms -are also taken into account (see Appendix A for the constraints on HNL extensions of the SM, and Appendix D for the relevant expressions for the cLFV observables, as well as the associated current bounds and future sensitivities.).The parameter space of the model is also now thoroughly explored: no assumptions are made on relations between mixing angles, nor on the heavy mediator masses.The former are now randomly and independently varied and drawing samples from log-uniform distributions, further randomly varying their signs, while the latter are no longer taken to be degenerate (only sufficiently close in mass to allow for interference effects 6 ).
In summary, the ranges of the parameters to be here explored are These ranges7 correspond to regimes for which the CP conserving cases comply with experimental bounds (cf.Appendix A).From the 10 4 points thus selected, and for each tuple of mixing angles, we   then vary all CPV phases associated with the sterile states, i.e. δ α4,5 , ϕ 4,5 ∈ [0, 2π] (with α = e , µ , τ ), drawing 100 values for each of the 8 from a uniform distribution.
were not considered as for the lower bounds these already correspond to predictions well beyond future experimental sensitivities.
In the following discussion, the obtained theoretical predictions are then compared with the existing current bounds (and prospects for future sensitivities) of the Z-and Higgs boson cLFV decays; these are summarised in Table 1.
A first overview of the prospects for cLFV Higgs decays is presented in Fig. 7, in which we display the rates of H → α τ (α = e, µ) vs. radiative and 3-body cLFV tau-lepton decays.We emphasise once again that now all active-sterile mixing angles, as well as Dirac and Majorana CP phases, are randomly varied, without any underlying simplifying assumptions.Notice that in what follows we will mostly focus on final state lepton pairs including one tau-lepton, as the decay for eµ final states is associated with comparatively smaller rates.
As expected from the underlying topologies of the different decays, H → α τ processes exhibit a significant correlation with τ → 3 α (notice that for heavy sterile states, the latter decays are in general dominated by Z-penguin contributions, which are also topologically similar to Higgs decays); this is in contrast to the radiative tau-lepton decays, τ → µγ.Although CPV phases (both Majorana and Dirac) lead to extensive interference effects, their impact is less striking than what had been found for purely leptonic processes (e.g.µ − e conversion vs. µ → 3e, see [44]).
Another conclusion to be drawn is that in the framework of this SM extension via 2 heavy neutral leptons (with masses around a few TeVs), cLFV Higgs decays are in all cases beyond future sensitivity; although larger rates could in principle be possible, the associated active-sterile mixing regimes are disfavoured from bounds on several electroweak precision observables (including invisible Z decays, among others) and cLFV τ → 3µ decays.In turn, this suggests that a near future observation of H → α τ decays would strongly disfavour this minimal class of HNL extensions of the SM.
A comparison between the prospects of cLFV Higgs and Z decays is presented in Fig. 8, in which we display H → α τ vs. Z → α τ , for α = e, µ.Both Z → eτ and Z → µτ decays are within future sensitivity of FCC-ee (running at the Z-pole), the latter offering the most promising prospects.As already anticipated, both observables exhibit a strong correlation in the CP conserving case (especially for µτ final states); such correlated behaviour is indeed affected by the presence of CP violating phases, which can induce both constructive and destructive interference.Nonetheless, and as mentioned before, the observed loss of correlation is milder than for purely leptonic processes [44].

CP asymmetries in cLFV boson decays
As mentioned in the Introduction, and extensively discussed in the present work, in models in which the SM is extended via heavy Majorana fermions, the "effective" Z and Higgs vertices to leptons are modified, allowing for both flavour and CP violation.In view of the capabilities and clean environment of a future FCC-ee (tera-Z factory), one can then consider to which extent such a minimal BSM construction could be at the source of non-vanishing contributions to CP-asymmetries, in particular concerning Z-boson decays, A CP (Z → α β ).The latter are defined as ).Dotted (dashed) lines denote current bounds (future sensitivity) as given in Table 1.(For additional information, see detailed description in the text.) As can be seen by the illustrative results of Fig. 9, obtained under the same simple assump- tions of the preliminary discussion of Section 4.1, one can indeed have non-negligible contributions to A CP (Z → α β ), induced by both Majorana and Dirac CPV phases.Individually, the former have a reduced impact, while the latter (especially δ 34 , through the C ij terms) can lead to very large asymmetries.As expected, and in general terms, the most promising channels appear to be those leading to final states containing one tau-lepton.It is also worth mentioning that the asymmetries are significantly dominated by the contributions to the cLFV Z decays stemming from the diagrams with 2 neutrinos in the loop (see Fig. 1 (a)).Also notice that the CP asymmetries computed in the limit of vanishing charged lepton masses are typically predicted to be larger; moreover their behaviour with respect to the CPV phases significantly differs.In Fig. 10, we explore the joint effects of Dirac and Majorana phases on A CP (Z → µτ ), still under simplifying assumptions for the mixing angles and masses of the heavy states.As can be seen, provided that Dirac phases are present, the Majorana phases have a significant impact (which was not the case should they be the only source of CPV, as seen from Fig. 9).
A thorough phenomenological study has been conducted for the CP asymmetries in the most promising channels of the cLFV Z decays, in particular those leading to final states containing one tau lepton (in view of the better prospects for observation of the decay itself).In Fig. 11 we display the CP asymmetries A CP (Z → α τ ) vs. the associated cLFV decay rates, BR(Z → α τ ), for α = e, µ. (See Section 4.2 for details on the underlying numerical scan of the parameter space.)Once all the relevant experimental and phenomenological constraints are taken into account, one can still have very large asymmetries (up to 100%) in both cases; notice however that such regimes are associated with rates for the cLFV process that are significantly beyond future sensitivities.For In the µτ sector, one can in fact simultaneously test the presence of additional heavy neutral leptons (and their CP violating phases) via several observables: this is displayed in Fig. 12, in which we emphasise the joint behaviour of three observables, A CP (Z → µτ ), BR(Z → µ ± τ ∓ ) and the pure leptonic decays BR(τ → µµµ).As can be seen from both panels (complementary views of the observables), for regimes leading to both BR(Z → µ ± τ ∓ ) and BR(τ → µµµ) within future sensitivities 8 , one can still have A CP (Z → µτ ) as large as 20% (cyan points).Although not a "smoking gun", the joint observation of the three observables could be interpreted as highly suggestive of such a extension of the SM via at least 2 heavy Majorana fermions (featuring new CPV phases).(Notice that for minimal "3+1" extensions via a single HNL, the cLFV rates are in general unaffected by the presence of CPV phases.)Finally, and to further highlight the impact of (potential) measurements of the CP asymmetries, we have selected the following CP conserving (P 1 ) and CP violating (P 2 ) benchmark points, specified as follows: Both benchmark points (CP conserving and CP violating) lead to common cLFV predictions as displayed in Table 2, thus rendering them indistinguishable if cLFV signals are observed in the future: in particular, notice that only the predictions for µ → eee, CR(µ − e, Al), τ → µµµ and Z → µτ lie within future sensitivities.Notice that for P 2 , which is associated with smaller mixing angles, CPV phases are at the source of constructive interferences, leading to predictions similar to P 1 (larger active-sterile mixings).From a phenomenological point of view, and since they share the same cLFV "profile", these two scenarios are in essence indistinguishable, and little can be learnt regarding the presence of leptonic  CP violating phases.However, the CP asymmetries in Z-boson decays offer a clear distinction between them: P 2 leads to the following predictions for the CP asymmetries In particular, A CP (Z → µτ ) is potentially observable, thus allowing to disentangle between CP conserving and CPV scenarios.Similarly to what is currently pursued for (semi-)leptonic LFV decays of mesons (e.g.B → K ( * ) ± α ∓ β ) -as motivated by leptoquark models -, we thus recommend a measurement of individual cLFV Z decays (i.e. with charge "tags" of the final state charged leptons).An independent measurement of both Z → µ − τ + and Z → µ + τ − rates would significantly help in disentangling CP conserving from CP violating scenarios of HNL models of new physics.
CP asymmetries in cLFV Higgs boson decays have also been considered.However, and in comparison with A CP (Z → α τ ), the scalar decays are far less promising.As discussed in the previous section, the associated branching ratios (for regimes fulfilling all imposed constraints) are typically beyond future experimental reach; the latter issue non-withstanding, one would nevertheless be led to maximal values of A CP (H → α τ ) much smaller than those found for Z decays, at most O(10 −12 ).Such a difference can be understood by considering the structure of the "effective Higgs vertex" (see Eq. ( 21)), where the only significant source of a CP asymmetry would stem 9 from the interference term (∝ Re(∆F L ∆F * R )); however this would be suppressed by the smallness of the light lepton masses.

Conclusions
Barring symmetries to enforce CP conservation, leptonic CPV phases are a generic feature of SM extensions via heavy neutral fermions.Arising in UV-complete models from various sources (among them Yukawa couplings), in the low-energy (effective) realisations both Dirac and Majorana phases can be present as intrinsic physical degrees of freedom of these SM extensions.Following a first assessment of the role of CP violating phases on cLFV observables [44], here we have carried out a thorough investigation of the impact of Dirac and Majorana phases on cLFV neutral boson decays.We have thus revisited Higgs and Z-boson decay rates, computing the widths for Z → α β and H → α β without any simplifying approximations.The full computations were done in both Feynman-'t Hooft and unitary gauges, and the complete expressions are given for future application in a user-friendly way.
Focusing on minimal SM extensions by heavy Majorana sterile fermions, the impact of the CPV phases on a number of cLFV decays was considered.In particular, we have illustrated the results for a minimal toy-model, in which 2 heavy Majorana sterile states are added to the SM content (without any assumption on the underlying mechanism of neutrino mass generation).Despite its minimality and simplicity, these constructions can be interpreted as a low-energy phenomenological limit of a complete high-energy model; the impact of the heavy states is a consequence of their masses, and of their mixings to the active states (including CP violating phases).Even if inferred for such a minimal SM extension, the conclusions of the present study are generic for heavy neutral lepton models (in which the new leptonic mixings are at the source of lepton flavour violation).In this sense, the results are comprehensive and generic, since the parameters are not related nor are certain regimes protected (or precluded) by any imposed flavour symmetry. 9Notice that, as shown by the expressions detailed in Section 3.2, one has F αβ L F βα R ; therefore any variation in FL is compensated by one in FR in the decay rate.No such effect is found for Z decays.
Confirming the findings of previous studies, cLFV Z and Higgs decays are very sensitive to the presence of the additional heavy neutral leptons; moreover, the CP violating phases have a clear impact on the decay rates, leading to both destructive and constructive interferences (again with the Dirac phases allowing for striking reductions -cancellations -of the widths).Following a first illustrative study, we have conducted a thorough phenomenological investigation, further emphasising the effect of CPV phases on (possibly) correlated behaviours between observables.
While cLFV branching ratios for Z → α β are potentially within reach of a future FCC-ee (especially for heavy states with masses O(1 TeV) or above), the prospects for the analogous Higgs decays are not very promising, as a consequence of the poorer future sensitivity.In the CP conserving limit, and as expected in view of the (leading) contributing topologies, one observes clearly correlated patterns between B → α β and α → 3 β (for B = Z, H).The presence of the CP phases (Dirac and Majorana) does affect the correlated patterns, although in a less striking manner than what occurs for the correlation of purely leptonic cLFV observables (cf.findings of [44]).
Motivated by the clean environment and future sensitivity of an FCC-ee (running at the Z-pole) for Z → α β transitions, we have also considered a distinct class of observables: the CP asymmetries associated with the latter cLFV boson decays.
Non-negligible -large -CP asymmetries in Z decays turn out to be a generic feature of HNL extensions encompassing at least 2 heavy states.For final states composed of a tau and a light charged lepton, with decay rates potentially within future sensitivity, one can have very large CP asymmetries; in particular, these can be as large as 20-30% for the case of A CP (Z → µτ ), interestingly in association with sizeable rates for τ → 3µ (also within future sensitivity).
As we have strongly emphasised, the presence of CP violating phases can open the possibility of having identical cLFV signatures in association with very distinct regimes of active-sterile mixings (a consequence of the interference effects due to the phases).Exploring the CP asymmetries in cLFV Z decays then offers a complementary probe to clearly assess the presence of CP violation.As highlighted, and despite not definitive, the joint observation of cLFV leptonic decays, together with cLFV Z decays and the associated CP asymmetry would strongly hint towards the presence of HNL-mediated flavour and CP violation.Moreover, and while a given cLFV observation can either stem from minimal extensions by only one sterile fermion ("3+1" model), the observation of a CP asymmetry is a clear signal of the presence of non-minimal HNL realisations (at least a "3+2"-like SM extension).
If on the one hand it is clear that CP violating phases should be in general taken into account upon comparison between prediction and observation in the context of cLFV HNL extensions of the SM, A CP (Z → α β ) might hold the key to clearly establishing the presence of leptonic CP violation; in turn, this might have strong implications regarding leptogenesis (relying on complete models including heavy sterile states).Whenever possible, data from individual channels (i.e.

A.1 Neutrino oscillation data
Below we collect the NuFIT 5.1 global fit results [62] for neutrino mixing data as used in our analysis.Although the predictions for cLFV transitions are mostly independent of the ordering of the spectrum and of the lightest neutrino mass, for concreteness we notice that in the numerical analysis we have assumed a normal ordered light neutrino spectrum, with the lightest neutrino mass in the range m 0 ∈ [10 −5 , 10 Table 3: Global fit results obtained by NuFIT 5.1 [62] for neutrino mixing data, not including experimental data from oscillation experiments measuring atmospheric neutrinos.In our numerical analysis we assume normal ordering of the light neutrino spectrum (∆m 3 = ∆m 31 > 0) and fix the neutrino mixing parameters to their central values.

A.2 Further phenomenological constraints
The mixing angles θ α4(5) between the active and heavy states can be constrained from precision observables which are insensitive to the new CPV-violating phases, and which can be divided into high-and low-energy ones.
At high energies, one can construct the lepton universality ratios of W boson decays [70] with α, β ∈ {e, µ, τ }.Moreover, the invisible decay width of the Z boson, Γ(Z → inv), plays a crucial role in constraining scenarios with heavy neutral leptons, as it is typically reduced when compared to the SM prediction.At low energies, several ratios from (semi-)leptonic decays of the τ -lepton and leptonic decays of light mesons can be constructed; these are sensitive to the modification of the W ν vertex in the presence of sterile fermions.In particular, one has [70,71] where R P ≡ Γ(P + → e ν) Γ(P + → µ ν) Additionally, we take into account upper bounds on the entries of η (see Eq. ( 1)) as derived in [72], taking into account indirect modifications of the Fermi constant G F , the weak mixing angle sin 2 θ w and the mass of the W boson, among others.For masses of the HNL at the TeV scale (as we are interested in our study), constraints from direct searches at colliders or from cosmology (such as big bang nucleosynthesis) are not competitive and will not be taken into account; nonetheless, we can derive theoretical constraints by imposing that decays of the new heavy states comply with perturbative unitarity [73][74][75][76][77], posing a direct bound on their decay width Γ ν i /m i < 1/2 for i ≥ 4. Since the dominant contribution to their decay stems from the W exchange, the bound can be written as where α w = g 2 w /(4π).Finally, the inclusion of sterile states leads to the modification of the prediction for the Majorana effective mass, m ee , to which the amplitude for neutrinoless double beta decay is proportional.In the presence of n S new states, m ee can be written as in which p 2 corresponds to the typical virtual momentum, p 2 −(100 MeV) 2 .We take into account the KamLAND-ZEN upper limit [78] m ee ≤ (61 ÷ 165) MeV.

B Feynman rules for SM extensions via Majorana sterile fermions
In Table 4, we summarise in a compact way the relevant interaction vertices, as inferred from the Lagrangian terms in Eq. ( 4).In the Z and Higgs vertices, the arrows denote the momentum flow.We note here that diagrams including at least one Zn i n j or Hn i n j vertex have to be symmetrised (factor 2) due to the Majorana nature of the physical neutrinos.
Table 4: Feynman rules for W , Z and Higgs interactions (and associated Goldstone bosons) in SM extensions via Majorana sterile fermions: n i denotes neutrino mass eigenstates (with i = 1 . . .5), while α corresponds to charged leptons.

C Additional information on leptonic Z and Higgs decays
In this appendix we discuss in detail different aspects of the computation of the cLFV Z and Higgs decays.We present sub-dominant terms not included in the main text, and we also include the results derived in the Feynman-'t Hooft gauge.
For the purpose of renormalisation, and in order to achieve finite and gauge-invariant results for the analytical calculation of the loop amplitudes, we work in the framework of low-scale type-I seesaw models [79][80][81][82][83].We thus hypothesise a neutrino Majorana mass matrix (complex and symmetric) of the form Notice that in the interaction basis the active neutrino mass vanishes due to gauge invariance.Light neutrino masses are generated by Majorana mass terms from the heavy sector (with n S HNL) and Yukawa couplings.The full (Majorana) neutrino mass matrix M ν can be diagonalised by a unitary rotation Due to vanishing active neutrino masses in the interaction basis, the mixing matrix U has certain properties pertinent for the renormalisability of the model (which will be subsequently detailed).The charged lepton masses are assumed to be diagonal, such that the leptonic mixing matrix is then given by U.
In what follows, we provide detailed expressions for the cLFV Z-boson decay rates, also discussing several points related with the renormalisation of the contributions.

C.1 Full cLFV Z-boson decay rate
The full decay rate (adding the sub-leading contributions to the width given in Eq. ( 18)) is

C.2 Divergences in cLFV Z-boson decays
As mentioned in Section 3.1, terms proportional to the B 0,1 and C 00 functions contain UV-divergences which are regulated dimensionally (via an integration in D = 4 − 2ε dimensions).The UV-divergent pieces of the different contributions are given by (the superscripts associated with the distinct terms refer to the diagrams of Fig. 1): in which γ E is the Euler-Mascheroni constant (all other contributions are finite).Being a (semi-)unitary matrix, U fulfils the following identities: i,j i,j In renormalisable extensions of the SM with fermion singlets, such as the type-I seesaw (and its variants), the following additional identities hold: in which Mαβ (with α, β = 1, 2, 3) denotes the mass matrix of the active neutrinos in the interaction basis, which generally contains only vanishing elements.Using the identities in Eqs.(48)(49)(50)(51)(52)(53) and the tree-level definition of cos 2 θ w = M 2 W /M 2 Z , it can be easily shown that the sum of diagrams (a)-(d) is (only) finite if the Z boson is on-shell (i.e.q 2 = M 2 Z ).Thus, for the purpose of on-shell Z-boson decays, we can safely consider the limit D → 4 in Eqs.(8)(9)(10)(11)(12)(13)(14)(15)(16), and the (partial) amplitude is manifestly finite and independent of the 't Hooft renormalisation scale 10 .For the off-shell cLFV vertex, the corresponding box-diagrams (that contribute to the same process as the Z-penguin) have to be taken into account in order to yield a finite and gauge-independent amplitude.

C.3 Additional amplitudes for the cLFV Z-vertex
Here we present the additional "scalar" amplitudes that appear in the effective cLFV Z-vertex, but which do not contribute to the decay rates since they vanish due to the Ward identity.The contributions are given by (see Eq. ( 6) for their definition) with the Passarino-Veltman functions defined as in Eqs.(8)(9)(10)(11).The remaining (non-vanishing) contributions are given by with the Passarino-Veltman functions as defined in Eqs.(12)(13)(14)(15)(16)(17).Note that the "scalar" amplitudes are neither finite nor gauge invariant on their own; UV-divergences and gauge-dependent terms are however cancelled once the corresponding box-diagrams (that contribute to the same process as the off-shell Z-penguin) are taken into account.

C.4 Amplitudes for cLFV Z-boson decays in the Feynman-'t Hooft gauge
For completeness, we further include the amplitudes for the cLFV Z decays in the Feynman-'t Hooft gauge; the labels (a)-(d) refer to the topologies presented in Fig. 1, which now also include additional diagrams with Goldstone bosons. with in which where B α,β 0,1 = B 0,1 m 2 α,β , m 2 i , M 2 W .We can now compare the analytic expressions in both unitary (UG) and Feynman-'t Hooft (FG) gauges.After a substantial amount of algebra one can show that all differences in F R V , F L T and F R T between the two gauges (diagram by diagram) are proportional to a common factor of (D − 4) and are thus trivially 0 in D = 4 dimensions.For F L V the situation is slightly more involved due to a permutation of arguments in the two-point Passarino-Veltman functions.Neglecting terms independent of the internal fermion masses, and setting which can be shown to vanish relying on the identities in Eqs.(48)(49)(50)(51)(52)(53), and by using the following property of the two-point functions

C.5 Divergences in cLFV Higgs decays
As done for the Z-boson decays, we now briefly discuss divergences emerging in association with the cLFV Higgs decays diagrams: In a similar fashion to what was explained for the cancellation of the divergences in the Z-boson decay, using the identities in Eqs.(48)(49)(50)(51)(52)(53), it can be easily shown that the sum of diagrams (a)-(d) from Fig. 2 is finite.

C.6 Amplitudes for cLFV Higgs decays in the Feynman-'t Hooft gauge
As done before, and for completeness, below we collect the results for the cLFV Higgs decays in the Feynman-'t Hooft gauge (the labels (a)-(d) refer to the topologies presented in Fig. 2, which now also include additional diagrams with Goldstone bosons): with where with Comparing the analytic expressions in unitary and Feynman-'t Hooft gauges, one has, after setting q 2 = m 2 h and neglecting terms independent of the internal neutrino masses, which can be shown to vanish using the identities in Eqs.(48)(49)(50)(51)(52)(53)71).

D Leptonic cLFV observables
In this Appendix we collect several expressions concerning the decay widths of numerous cLFV leptonic decays and rare transitions, which are relevant for our numerical studies (constraints and implications for probing the model).Likewise, in Table 5, we summarise the bounds and future sensitivities for these cLFV processes.Table 5: Current experimental bounds and future sensitivities on cLFV observables here considered.The quoted limits are given at 90% C.L. (Belle II sensitivities correspond to an integrated luminosity of 50 ab −1 ).

D.1 Radiative and three-body decays
Below we summarise the expressions for the radiative and β → 3 α decays (the full expression for the most general case of the 3-body decay can be found in [19]), closely following the notation of [44].
The rates for the radiative and three-body decays in the SM extended via n S heavy sterile fermions, are given by [18] BR( As before, M W is the W boson mass, α w = g 2 w /4π denotes the weak coupling, s w the sine of the weak mixing angle, and m β (Γ β ) the mass (total width) of the decaying charged lepton of flavour β.The form factors G βα γ , F βα γ , F β3α box are given by [18,19] In the above expressions, the sums are understood to be taken over all neutral mass eigenstates (i, j = 1, ..., 3 + n S ).We recall that the C ij function is defined as The distinct loop functions (with arguments defined as x i = m 2 i /M 2 W ) are summarised in Appendix D.3, with the corresponding arguments defined as x i = m 2 i /M 2 W .

D.2 cLFV in muonic atoms
Although other cLFV transitions and decays can also occur in the presence of muonic atoms (as for example Muonium oscillations and decays, or the Coulomb enhanced decay µe → ee), here we collect the expressions for the coherent conversion rate in the presence of a nuclei (N), as the latter turns out to be among the most relevant cLFV observables due to its constraining power.The neutrinoless conversion rate is given by [18] CR(µ − e, N) = 2G in which Γ capt.denotes the capture rate for the nucleus N, with D, V (p) and V (n) corresponding to nuclear form factors (see [97]), and e is the unit electric charge.The above form factors are given by [18,19] to which one must add Here, x q = m 2 q /M 2 W and V is the Cabibbo-Kobayashi-Maskawa (CKM) quark mixing matrix.All relevant loop functions can be found in Appendix D.3.

D.3 Loop functions
For completeness we collect here the loop functions for the purely leptonic cLFV decays and transitions discussed above, as well as some relevant limits (as presented in [18,19]).

Photon dipole and anapole functions
Z-penguin: two-and three-point functions x − y (4 + xy) In the above, we highlight that F Xbox corresponds to F box in [19]; moreover it has an opposite global sign in comparison to [18], which further impacts F β3α box (see [18,19]).

Figure 5 :Figure 6 :
Figure 5: Rates of cLFV Z-boson decays as a function of the CP violating Dirac and Majorana phases.From left to right and top to bottom, dependence on δ 14 , δ 24 , δ 34 and ϕ 4 (all other phases set to zero in each case).The colour code denotes the flavour composition of the final state lepton pair: eµ (blue), eτ (orange) and µτ (green).In all panels, solid, dashed and dotted lines respectively correspond to the following heavy fermion masses: m 4 = m 5 = 1, 5, 10 TeV.

Table 2 :
cLFV predictions (cLFV "profile") associated with benchmark points P 1 and P 2 (see text for details); observables not displayed are associated with rates lying beyond future sensitivity.