Probing Leptogenesis at Future Colliders

We investigate the question whether leptogenesis, as a mechanism for explaining the baryon asymmetry of the universe, can be tested at future colliders. Focusing on the minimal scenario of two right-handed neutrinos, we identify the allowed parameter space for successful leptogenesis in the heavy neutrino mass range between $5$ and $50$ GeV. Our calculation includes the lepton flavour violating contribution from heavy neutrino oscillations as well as the lepton number violating contribution from Higgs decays to the baryon asymmetry of the universe. We confront this parameter space region with the discovery potential for heavy neutrinos at future lepton colliders, which can be very sensitive in this mass range via displaced vertex searches. Beyond the discovery of heavy neutrinos, we study the precision at which the flavour-dependent active-sterile mixing angles can be measured. The measurement of these mixing angles at future colliders can test whether a minimal type I seesaw mechanism is the origin of the light neutrino masses, and it can be a first step towards probing leptogenesis as the mechanism of baryogenesis. We discuss how a stronger test could be achieved with an additional measurement of the heavy neutrino mass difference.


Introduction
Motivation. The Standard Model (SM) of particle physics allows to describe almost all phenomena in nature at a fundamental level [1]. Neutrino flavour oscillations, which clearly indicate the existence of neutrino masses, are the only phenomenon observed in the laboratory that points unambiguously towards the existence of new states beyond the SM. If the neutrino masses are at least partly generated by the Higgs mechanism in the same way as the masses of all other fermions, then this necessarily implies the existence of right handed neutrinos ν R . Since the ν R are singlets under the gauge symmetries of the SM, they are allowed to have a Majorana mass term − 1 2 ν c R M ν R in addition to the usual Dirac mass term that is generated by the Higgs mechanism. If there are n s right-handed neutrinos ν Ri , then M is an n s × n s flavour matrix with eigenvalues M i . These Majorana masses M i are free parameters. Their magnitude cannot be fixed by neutrino oscillation data alone, since light neutrino oscillations are only sensitive to a particular combination of M i and the right-handed neutrinos' Yukawa couplings Y ia to SM leptons of flavour a, cf. eq. (2.5). The range of masses allowed by neutrino oscillation data is in principle very large; even in the minimal model with n s = 2 it reaches from the eV scale [2] up to values 1-2 orders of magnitude below the Planck mass [3]. The implications of the existence of right-handed neutrinos for particle physics and cosmology strongly depend on the magnitude of M , see e.g. ref. [4] for a review. Traditionally it is often assumed that the M i are much larger than the electroweak scale. In this case the smallness of the observed neutrino masses can be explained by the usual seesaw mechanism [5][6][7][8][9][10], i.e., the smallness of v/M i , where v = 174 GeV is the vacuum expectation value of the Higgs field. However, this choice is entirely based on theoretical arguments that e.g. relate M i to the scale of grand unification.
Alternatively one could e.g. argue that the M i and electroweak scale have a common origin [11][12][13][14]. A value of M i at (comparably) low scales is technically natural because B−L becomes a symmetry of the model in the limit where all M i vanish, and large radiative corrections to the Higgs mass that plague high scale seesaw models can be avoided. Many low scale seesaw models indeed involve an approximate "lepton number"-like symmetry. The smallness of the light neutrino masses in these symmetry protected seesaw scenarios does not primarily come from the seesaw suppression by the parameters v/M i (which are not small), but is due to the smallness of the symmetry violation, cf. eqs. (2.21). 1 A benchmark scenario can be found in [15].
The probably most studied model that involves a (sub-) electroweak scale seesaw is the Neutrino Minimal Standard Model (νMSM) [16,17]. This was indeed the model in which it was first discovered that leptogenesis is feasible for M i < v in the minimal scenario with n s = 2 [17], which is the setup that we are concerned with in the present paper. The νMSM realises the an approximate B − L symmetry in part of its parameter space [18]. Other popular models of this type are often referred to as "inverse seesaw models" [19][20][21][22], "linear seesaw models" [23][24][25][26] (see also [27][28][29][30][31]) and "minimal flavour violation" [32,33]. In the present paper we take an agnostic approach to the magnitude of the M i and focus on the mass range that is accessible to near future experiments.
The Yukawa couplings Y of the right-handed neutrinos ν R in general violate CP, while M violates lepton number L. Since lepton number L and baryon number B can be converted into each other by electroweak sphaleron processes in the early universe [34], this opens up the possibility that they are the origin of the baryon asymmetry of the universe (BAU) [35] and thereby solve one of the great mysteries in cosmology that cannot be explained within the SM. 2 The experimentally observed BAU can be quantified by the baryon-to-entropy ratio [37] Y B obs = (8.6 ± 0.1) × 10 −11 . (1.1) While leptogenesis generically requires rather large M i [38], an approximate B − L symmetry can alleviate this requirement in different ways [17,[39][40][41]. For M i below the TeV scale, and within the minimal seesaw model, there are three different mechanisms that generate lepton asymmetries. For M i above the electroweak scale, the baryon asymmetry can be generated in heavy neutrino decays via resonant leptogenesis [39]. The lower bound on the mass comes from the requirement that the heavy neutrinos freeze out and decay before sphalerons freeze out at T sph ∼ 130 GeV [42]. For masses M i below the electroweak scale, the baryon asymmetry can be produced via CP-violating oscillations of the heavy neutrinos during their production (instead of their decay) [17,43]. This mechanism is also known as baryo-or leptogenesis from neutrino oscillations. Finally, there is a contribution to the asymmetries from the lepton number violating (LNV) decays of Higgs quasiparticles with large thermal masses into ν Ri and SM leptons [44,45].
Goals of this work. In the present paper we investigate the perspectives to probe low scale leptogenesis in the minimal seesaw model with n s = 2 at the proposed future lepton colliders, the electron-positron mode of the Future Circular Collider (FCC-ee) at CERN [46], the Circular Electron Positron Collider (CEPC) in China [47] and the International Linear Collider (ILC) in Japan [48,49]. We focus on the mass range 5 GeV < M i < 80 GeV, i.e., on heavy neutrinos that are heavier than b-mesons and lighter than W bosons. In this regime, large numbers of heavy neutrinos can be produced (see e.g. [15,[50][51][52][53][54][55][56][57][58]) from on-shell Z bosons at the so-called Z pole run of a lepton collider, or from W boson exchange at higher center-of-mass energies.
• The interaction strength of the heavy neutrinos with SM leptons of flavour a can be characterised by an active-sterile mixing angle U 2 a . In previous studies, the potential of future experiments to discover leptogenesis has been estimated by comparing the projection of the viable leptogenesis parameter region on the M i − U 2 a planes to the 3 Previous studies suggest that the LHC cannot probe the parameter region where leptogenesis is possible in the minimal model [65][66][67][68][69], or may just touch it [62]. The potential to search for heavy neutrinos with larger masses has recently e.g. been studied for both, hadron [50,67,[70][71][72][73][74][75][76][77] and lepton colliders [15, 50, 52-55, 76, 78].
projection of the experimental sensitivity on these planes. This comparison of projections is not fully consistent because both, low scale leptogenesis and the experimental sensitivities, depend not only on the total interaction strength U 2 = a U 2 a , but also on the relative size of U 2 e , U 2 µ and U 2 τ . 4 For instance, the fact that a given combination of M i and U 2 µ is consistent with both (successful leptogenesis and an observable event rate at a collider) does not necessarily imply that those heavy neutrinos that are able to generate the BAU can actually be discovered at the collider because the parameter points for which one or the other can be realised may correspond to very different U 2 e . Using measurements of the U 2 a at a lepton collider in order to determine phases in U ν has been investigated in [94]. We note that for leptogenesis, not only the relation between the U 2 a and the amount of CP violation but also the washout of the different active lepton flavours is crucial. In the present analysis, we calculate the BAU and the expected number of events for each point in the model parameters space to assess the question whether both requirements can be fulfilled simultaneously in a consistent manner.
• We estimate how precisely the experiments can measure the magnitude of the individual U 2 a . If any heavy neutral leptons are discovered in future experiments, the relative size of the U 2 a provides a powerful test whether these particles can generate the light neutrino masses and the BAU in the minimal seesaw model with two RH neutrinos (cf. [64,86] for previous studies). The sensitivity of an experiment to individual U 2 a is therefore important to assess the experiment's potential to not only discover heavy neutral leptons, but understand their role in particle physics and cosmology.
• In addition, we also investigate which values of the heavy neutrino mass splitting are consistent with leptogenesis and discuss strategies to measure them at future colliders.
This article is organised as follows: In section 2 we recapitulate the basics of the seesaw mechanism and the symmetry protected seesaw scenario. In section 3 we provide details of our scan of the leptogenesis parameter space. In section 4 we specify our approach to assess the reach of future lepton colliders in terms of the M i and U 2 a . In section 5 we present and discuss the results of our numerical analysis. We conclude in section 6. In appendix A we give details on the expected number of events at the given lepton collider, while the statistics behind the precision of the measurements of the mixings is explained in appendix B. A derivation of the evolution equations for the heavy neutrinos including lepton number violating processes is given in appendix C. The Lagrangian of the type-I seesaw model, is obtained by extending the SM Lagrangian L SM by n s right handed (RH) spinors ν Ri , which represent the RH neutrinos. The interaction with the SM only happens via the Yukawa interactions Y ia with the SM lepton doublets a for a = e, µ, τ and the Higgs field φ. The superscript c appearing on the RH neutrinos denotes charge conjugation and ε is the antisymmetric SU(2)-invariant tensor with ε 12 = 1. Further M ij are the entries of the Majorana mass matrix M of the RH neutrinos with eigenvalues M i . At tree level the connection between the Lagrangian (2.1) and the low energy oscillation data can be provided by the Casas-Ibarra parametrisation [95] Here (m diag ν ) ab = δ ab m a denotes the light neutrino mass matrix in the mass basis, i.e., m a are the light neutrino masses. The mixing matrix of the light neutrinos is given by with v the temperature dependent Higgs expectation value evaluated at zero temperature: v = 174 GeV. The unitary matrix U ν diagonalises the light neutrino mass matrix with U ±δ = diag(1, e ∓iδ/2 , e ±iδ/2 ). The non-vanishing entries of V (ab) for a = e, µ, τ are given by where θ ab are the neutrino mixing angles. In order to fix these light neutrino oscillation parameters in U ν , we neglect the non-unitarity in eq. (2.3) and use the parameters given in table 1. Further, δ is the so-called Dirac phase, while α 1,2 are referred to as the Majorana phases. In the case of two generations of sterile neutrinos only one combination of the Majarona phases is physical, i.e. for normal hierarchy (NH) it is α 2 , while for inverted hierarchy (IH) it is α 2 − α 1 . Without loss of generality we can set α 2 ≡ α and α 1 = 0. The mixing of the doublet state ν L with the singlet state ν R implies a deviation of V ν from unitary. The complex matrix R in eq. (2.2) fulfils the orthogonality condition R T R = 1. The number n s of RH neutrinos is unknown and not restricted by anomaly requirements on SM interactions because they carry no SM charges. The fact that two non-zero mass splittings m 2 sol and m 2 atm have been observed implies that n s ≥ 2 if the seesaw mech- anism is the sole origin of light neutrino masses because the number of non-vanishing eigenvalues m a in equation (2.5) has to be smaller than or equal to the number of generations of heavy neutrinos. In the following we focus on the minimal scenario with n s = 2, in which the lightest neutrino is massless (m lightest = 0). This effectively also describes leptogenesis and neutrino mass generation in the νMSM. 5 For n s = 2 one can parametrise the matrix R by via a complex mixing angle ω = Reω + iImω where ξ = ±1. After electroweak symmetry breaking there are two distinct sets of mass eigenstates, which we represent by Majorana spinors ν i and N i . The three lightest can be identified with the familiar light neutrinos, In addition, there are n s heavy neutrinos (2.10) These expressions are valid up to second order in |θ ai | 1. The heavy neutrino mass matrix gets diagonalised by the unitary matrix U N , and we can define V N = (1 − 1 2 θ T θ * )U N in analogy to V ν . Even though the heavy neutrinos are gauge singlets they are able to interact with the ordinary matters due to their quantum mechanical mixing. Therefore, any process that involves ordinary neutrinos can also occur with heavy neutrinos if it is kinematically allowed, but the amplitude is suppressed by the mixing angle Thus, it is convenient to express the branching ratios in terms of It is well known that low scale leptogenesis with only two heavy neutrinos requires a mass degeneracy |∆M | M , wherē (2.14) The physical masses observed at colliders are given by the eigenvalues of M N in eq. (2.11) and includes a contribution O[θ 2 ] from the Higgs mechanism. Their splitting ∆M phys can be expressed in terms of the model parameters as with ∆M θθ = m 2 − m 3 for normal ordering and ∆M θθ = m 1 − m 2 for inverted ordering. If the mass splitting ∆M phys is too tiny to be resolved experimentally, experiments are only sensitive to the quantities The overall coupling strength of the heavy neutrinos can be expressed in terms of the Casas-Ibarra parameters as in case of normal hierarchy and for inverted hierarchy.

Symmetry protected scenario
Many models that motivate a low scale seesaw exhibit additional "lepton number"-like symmetries that make small M i with comparatively large neutrino Yukawa couplings technically natural. Such symmetry protected scenarios, see e.g. [15] for a benchmark scenario, are phenomenologically very interesting because they allow to make mixings U 2 ai much larger than suggested by the "naive estimate" In the symmetry protected scenario it is convenient to use the parameters instead of ∆M and Imω. ForM near the electroweak scale, experimental constraints allow individual Yukawa couplings |Y ia | that are larger than the electron Yukawa coupling [64,98], and the smallness of the m i is primarily a result of the smallness of and µ. Specific models that invoke µ, 1 typically predict a relation between these parameters, i.e., specify a path in the -µ plane along which the limit µ, → 0 should be taken. 6 In the present work we take an agnostic approach and treat and µ as independent parameters.
An approximate B − L symmetry emerges in the limit where these parameters are small. This can be made explicit by applying the rotation matrix to the fields ν Ri in eq. (2.1), which brings M and Y into the form where ε a Y a < 1 7 . When setting ε a → 0, one can assign lepton number +1 and −1 to the states where ν Ri refer to the flavour eigenstates of M . In the limit µ, → 0 the state ν Rw decouples. This scenario is often called pseudo-Dirac scenario because the Lagrangian can be expressed in terms of a single Dirac spinor ψ N = (ν Rs + ν c Rw ), The second line in eq. (2.26), which summarises the LNV terms, vanishes in the limit 6 While there is nothing that forbids setting µ = 0, εa must remain finite to ensure that the neutrino masses mi > 0. 7 Note that the entries εa are actually of order O[ √ ] in the parameter defined in eq. (2.21). We use these conventions to be consistent with the notation commonly used in the literature. µ, → 0. In the mass basis we find we find Y 1a = iY 2a = Y a / √ 2 in this limit, and therefore The ratios U 2 a /U 2 are mostly determined by the parameters in U ν [33,64,86,[99][100][101][102]. To leading order in m sol /m atm , and θ 13 [33,64,86,103,104] they can be approximated as for NO, (2.28) This is the limit in which we perform the collider analysis in section 4. Strictly speaking we should distinguish between the lepton number L carried by the SM fields and a generalised lepton numberL that includes the lepton charges associated with some "lepton number"-like symmetry.L is conserved in the limit µ, → 0, while L is still violated by active-sterile oscillations in this limit. L is only conserved if one in addition setsM = 0, in which case the seesaw approximations that we use cannot be applied. To simplify the notation and stick to the commonly used terminology, we in the following refer to the B −L conserving limit µ, → 0 as "approximate B − L conservation".

Leptogenesis
In leptogenesis, the matter-antimatter asymmetry of the universe is generated in the lepton sector and transferred to the baryonic sector via weak sphalerons [34]. These processes are only efficient for temperatures above T sph = 130 GeV, below which the baryon number B is conserved. In the standard Leptogenesis scenario [35], often referred to as thermal leptogenesis, the baryon asymmetry is generated due to the decay of the heavy neutrinos with M i above the electroweak scale such that the N i are not light enough to be produced efficiently at lepton colliders in the decay of real weak gauge bosons. Instead, the baryon asymmetry can be generated via CP-violating flavour oscillations of the N i (baryogenesis from neutrino oscillations) [43] or the decay of Higgs bosons [17]. The main difference between standard leptogenesis and these mechanisms lies in the way how the deviation from thermal equilibrium, which is a necessary condition for baryogenesis [105], is realised. Superheavy N i come into equilibrium at temperatures T T sph . In this case the deviation from equilibrium that is responsible for the asymmetry we observe today is caused by their freezeout and decay, which must occur at T > T sph in order to be transferred from a lepton into a baryon asymmetry ("freeze out scenario"). For M i below the electroweak scale, on the other hand, the relation (2.5) implies that the Y ia can be small enough that the N i do not reach thermal equilibrium at T < T sph ("freeze in scenario").

Full system of differential equations
Characterisation of the asymmetries. It is convenient to describe the asymmetries by the quantities that are conserved by all SM processes including the weak sphaleron transitions. Here B is the baryon number and L a = (g w q a + q Ra ) are the individual lepton flavour charges stored in the right handed leptons q Ra and the doublet leptons q a . The total SM lepton number is given by The L a and L are violated by the Yukawa couplings Y ia . The smallness of the Y ia implies that the ∆ a evolve very slowly compared to the rate of gauge interactions that keep the SM fields in kinetic equilibrium, so that they can effectively be described by chemical potentials, cf. eq. (C.37). Due to their Majorana nature the N i do not carry any charge in the strict sense, but in the limit T M i the helicity states of the heavy neutrinos effectively act as particle and antiparticles. This allows to define the sterile neutrino charges in terms of the difference between the occupation numbers of N i with positive and negative helicity, i.e., the diagonal elements of the matrices n hij defined in eq. (C.28) in the heavy neutrino mass basis, It is useful to introduce the generalised lepton number This quantity is approximately conserved in the temperature regime T M i because helicity conserving processes are suppressed when the N i are relativistic. It is in general not equivalent to the lepton numberL that is conserved in the limit , µ → 0. Though the violation ofL is parametrically small, it can be significant in the early universe because heavy neutrino oscillations generally convert ν Rs and ν Rw into each other, andL is therefore not a useful quantum number to describe leptogenesis. Only in the highly overdamped regime these conversions are suppressed, and one can identify ν Rw with the slowly evolving state.
Quantum kinetic equations. The evolution equations for the deviations of the heavy neutrino number densities from equilibrium δn h form a coupled system of equations together with the evolution equations for the asymmetries. A consistent set of equations to describe low scale leptogenesis has first been formulated in ref. [17] as a generalisation of the density matrix equations that are commonly used in neutrino physics [106]. In this work we employ a modified version of the quantum kinetic equations used in ref. [64], which have been derived in ref. [88]. The main improvement in the present work is the inclusion of LNV processes, which were neglected in most previous studies. Here we only present the quantum kinetic equations, the modifications to the derivation in ref. [88] that are necessary to derive them are presented in appendix C. We use the time variable z = T ref /T , where T ref is an arbitrarily chosen reference scale. It is convenient to set this to the sphaleron freeze out temperature T sph such that the freeze out happens at z = 1. In terms of z, the evolution equations for the deviations of the N i occupation numbers from equilibrium read where a R = m Pl 45/(4π 3 g * ) = T 2 /H can be referred to as the comoving temperature in a radiation dominated universe, H is Hubble parameter, g * is the number of relativistic degrees of freedom in the primordial plasma and m Pl the Planck mass. The collision term for the N i is conventionally decomposed into two contributions: The damping term Γ N and the backreaction termΓ N . The latter pursues chemical equilibration of the heavy neutrinos with the Higgs field and the SM lepton doublets. The source term for the lepton charges S a ≡ S aa can be obtained from the more general expression and is responsible for the creation of the lepton doublet charges due to the off-diagonal correlations δn ij . The momentum independent rates γ ± with γ + = γ av + and γ − = γ |kav| − that appear in S a are discussed in appendix C. The term in brackets in eq. (3.6) is the washout term. The matrices A, D as well as the vector C account for the fact that spectator effects redistribute the charges [107]. They are given by where the following notations are used: (3.14) Further, the momentum independent ratesγ ± withγ + = γ av + andγ − = γ av − that appear in the termΓ a N are discussed in appendix C. These expressions for H th N , Γ N andΓ a N include LNV effects only at leading order inM /T and are valid in the case of two relativistic heavy neutrinos with kinematically negligible mass splitting. This improves the accuracy of our treatment compared to the previous analysis in ref. [64], as explained in section 3.2. There are, however, several effects that we still neglect. This includes the kinematic effect of particle masses (gauge boson, top quark and N i ), the temperature dependent continuous freeze out of the weak sphalerons (we assume an instantaneous freeze out at T = T sph 8 ) and the error that occurs due to the momentum averaging. The latter is briefly discussed after eq. (C.28). We expect that these effects are subdominant in most of the parameter region we study here. However, they may become important if either the baryon asymmetry is generated shortly before the sphaleron freeze out or if the heavy neutrinos have masses comparable to the W boson.
Collision terms. The lepton number conserving and violating contributions to the collision terms exhibit a different dependence on the Yukawa couplings and on T . The different dependency on the Y ia can be expressed in terms of the quantities Υ +h and Υ a +h (lepton number conserving) or Υ −h and Υ a −h (lepton number violating). The exact T dependence is in principle rather complicated because various different processes contribute to the equilibration rates. Here we employ the commonly used linear approximation for the T dependence of the momentum averaged rates, which hides all the complications in the numerical coefficients γ av + , γ av − and γ |kav| − . For the lepton number conserving rates we use the values γ av + = 0.012 given in ref. [109], based on the results obtained in refs. [110,111]. Lepton number conserving processes are highly dominated by hard momenta |k| ∼ T of the heavy neutrinos, such that there is no significant difference between evaluating the rates at the average momentum |k av | ≈ 3.15T and momentum averaging the rates. Note that the production rate γ + and the rate for the backreaction termγ + are in general different since they come with different powers of the equilibrium distribution function before momentum averaging. Equating them induces an error of roughly 30 %, an approximation that is still sufficient for our discussion.
The effect of momentum averaging is different for the lepton number violating processes because the corresponding rates come with an additional factor of M 2 /|k| 2 . This results in an infrared enhancement of these rates. Consequently, we have to distinguish between the lepton number violating rate that is momentum averaged γ av − and the one that is evaluated at the averaged momentum |k av |. Approximate numerical results are derived in appendix C. We use for the backreaction termΓ a N , as it depends on quantities that are in equilibrium. In contrast to that, we use for the terms that depend on the deviations δn, such as the source term S a and the production term Γ N .
Thermal corrections to the effective Hamiltonian. In the absence of helicity flips, the thermal correction to the heavy neutrino masses is simply given by the term involving h th + = 0.23, as mentioned in ref. [88], plus a contribution arising from the expectation value of the Higgs field We evaluate the latter term with the approximation (C.17) that has already been used in ref. [88]. Lepton number violating forward scatterings generate an additional correction The derivation of this expression is sketched in appendix C.

The role of lepton number violating processes
The lepton numberL is violated by the Majorana mass term M . For M i below the W mass, the BAU is generated when the N i are relativistic. The rates ofL-violating processes in this regime are suppressed by M 2 i /T 2 and have therefore been neglected in most previous studies. This is, however, not justified in general. There are two important effects arising from theL-violating processes, the first is a lepton number violating source term, the second is the enhanced equilibration of the right handed neutrino states in the symmetry protected regime.
To understand the importance of theL violating source let us first consider the case ∼ 1, i.e., the "naive seesaw limit" (2.20). There are two competing sources of lepton asymmetries, the CP-violating oscillations of the N i [17,43] and the decay of Higgs bosons with large thermal masses [44] (cf. also ref. [112] for an earlier discussion). The heavy neutrino oscillations do not directly change L. This can be understood intuitively in terms of the suppression of LNV by M 2 i /T 2 and is shown in detail in appendix D.3 of ref. [88]. A total L = 0 (and hence B = 0) is generated with the help of the washout, which erases the individual L a at different rates. This generates a total L = 0 (and hence B − L = 0) even ifL is conserved because it stores part of the asymmetry in the sterile flavours, where it is hidden from the sphalerons. Since the washout is also mediated by the Yukawa couplings Y ia , the final baryon asymmetry in the regime ∼ 1 is O[Y 6 ], cf. also ref. [113] for a pedagogical discussion. In contrast to that, the Higgs decays directly violate L and L, which leads to a contribution O[Y 4 M 2 i /T 2 ] to the baryon asymmetry. 9 In the regime ∼ 1 the seesaw relation (2.5) predicts a |Y ia | 2 ∼ m a M i /v 2 . Hence, the Higgs decays can dominate the baryon asymmetry if it is primarily generated at temperature scales lower than O(v M /m a ). Furthermore, as the asymmetry generated this way is a totalL asymmetry it cannot be erased by the usualL-conserving rate and can therefore dominate the baryon asymmetry even for 1. The second effect is the enhanced equilibration of the heavy neutrino eigenstates, which is most important in the symmetry protected limit, where the matrix Y Y † has two vastly different eigenvalues, the magnitudes of which scale as a ε 2 a ∼ and a Y 2 a ∼ 1/ . The damping of deviations in the heavy neutrino occupation numbers from equilibrium is governed by the eigenvalues of Γ N . Let us for a moment assume that there are noLviolating processes. Then we can define the heavy neutrino interaction eigenstates as the eigenvectors of the helicity dependent flavour matrices Υ +h . In the limit µ, → 0, they can be identified with the states ν Rs and ν Rw that carry the generalised lepton chargeL, cf. eq. (2.25). One pair of interaction eigenstates (approximately ν Rs ) has comparably strong couplings Y 2 a ∝ 1/ , while the couplings of the other pair (approximately ν Rw ) are suppressed by ε a ∝ √ . This leads to an overdamped behaviour of the N i oscillations because the more strongly coupled states come into equilibrium before the heavy neutrinos have performed a single oscillation. The deviations from equilibrium which drive baryogenesis are then given by the slow evolution of the feebly coupled states. The feebly coupled state instead reaches equilibrium through the mixing with the strongly coupled state. In the presence ofL-violating processes, both eigenvalues of Γ N receive corrections from the terms involving Υ −h . For the larger eigenvalue, these can be neglected. However, the correction ∼ Y 2 aM 2 /T 2 ∝ 1/ ×M 2 /T 2 to the smaller eigenvalue, which governs the damping of the feebly coupled interaction eigenstate, is not necessarily negligible compared to the terms ∼ ε 2 a ∝ in Υ +h . For sufficiently small they dominate over the lepton number conserving terms involving Υ +h , which are suppressed by . Since the deviations from equilibrium in the overdamped regime are mostly determined by the occupation numbers of the feebly coupled state, this modification of the damping rates has a strong effect on the behaviour of the entire system of equations and affects the BAU. One can roughly estimate that this can affect the asymmetry if this occurs before sphaleron freezeout ( <M /T sph ) if the BAU is generated in the overdamped regime.

How to perform the parameter scan
We perform a parameter scan by numerically solving the evolution equations (3.5) and (3.6) in order to identify the range of heavy neutrino parameters for which the observed BAU can be explained. We limit ourselves to masses in the range between 5 GeV <M < 50 GeV. At smaller M i , fixed target experiments offer much better perspectives to search for heavy neutrinos than high energy lepton colliders. At larger M i , the approximations in the derivation of the expressions (3.9)-(3.12) are not justified. The scan is performed as follows.
We first randomly choose a value forM between 5 GeV and 50 GeV with a logarithmic prior and a value for Im ω with a flat prior between 0 < Imω < 7. This corresponds to a logarithmic prior in U 2 ∼ e 2Imω . The upper limit Imω < 7 does not limit the scan; in practice we find that the observed BAU cannot be produced for larger values of Imω due to a stronger washout.
After fixingM and Imω, we perform a simple Markov-chain-Monte-Carlo (MCMC) search using the Metropolis-Hastings algorithm over the remaining parameters α, δ, ∆M and Reω. The proposal distribution is a multivariate gaussian distribution in α , δ , Re ω and log ∆M/M , while the acceptance distribution is given by where |Y B | is obtained by numerically solving the evolution equations (3.5, 3.6). In the final analysis we accept all parameter choices that give a BAU within the 5σ obs region of the observed BAU. As the largest mixing angles require a hierarchical flavour pattern U 2 a U 2 , we in addition perform targeted scans in which the initial values of parameters α and δ for the MCMC scan are chosen to minimise the ratio U 2 a /U 2 . These points yield the largest numbers of events one can expect at an experiment.
The upper bound on the mixing angle U 2 in figure 3 is determined by binning the data points consistent with the BAU according to the logarithm of the mass logM into 60 bins, and in each bin choosing the point with the largest mixing angle U 2 . If the N i are produced in the decay of Z bosons in the s channel, then the total number of expected collider events (to be explained in the next section) in a given experiment and for fixedM in good approximation only depends on U 2 , cf. figures 9 and 10. For the Z pole runs, we can therefore uniquely define the area where one can expect more than 4 events in figure 3. If the N i are produced via t-channel exchange of W bosons, then the mixing in at least one of the vertices must be U 2 e because the experiment collides electrons and positrons. The total event rate therefore depends on U 2 e in a different way than on U 2 µ and U 2 τ , and the number of events cannot be determined by fixingM and U 2 alone, cf. figure 11. In figure 3 we therefore indicate two regions: The one where the expected number of events exceeds four under the most pessimistic assumptions about the relative size of the U 2 a for fixed U 2 ("guaranteed discovery"), and the one where it exceeds four under the most optimistic assumptions ("potential discovery"). The lines corresponding to a guaranteed discovery are obtained by picking the smallest mixing angle U 2 in each bin. To obtain the lines with a potential discovery, we instead select the points where the number of expected events is N < 4, and from the bins select those with the largest mixing angle U 2 . For the plots showing the maximal/minimal number of expected events in appendix A, we divide the points into 60 bins according to the logarithm of the massM , as well as 60 bins according to the logarithm of their mixing angle U 2 . From each bin we then select the points with the maximal and minimal numbers of events expected at the future collider considered.

Measurement of the leptogenesis parameters at colliders
In this section we discuss the possibility of measuring the neutrino parameters at possible future high-precision lepton colliders. The precise knowledge of the mixing quantities U 2 a and the flavour mixing ratios U 2 a /U 2 is crucial to test whether the properties of a hypothetical heavy neutral lepton that is discovered at a future collider are compatible with leptogenesis and the generation of light neutrino masses.
The upper bound on the active-sterile mixing angles U 2 a that is consistent with the observed baryon asymmetry of the universe in the minimal seesaw model results in comparatively long lifetimes in the range of picoseconds to nanoseconds for the heavy neutrinos with masses between a few GeV and the W boson mass [114]. Therefore, the heavy neutrinos produced in the particle collisions have a long enough lifetime to travel a finite distance before they decay inside the detector, giving rise to a displaced vertex. The origin of the displaced vertex signature by heavy neutrinos 10 at particle colliders is shown schematically in figure 1. Such an exotic signature allows for extremely sensitive tests of the active-sterile mixing angles for heavy neutrino masses below ∼ 80 GeV, especially for future lepton colliders with their high integrated luminosities, see e.g. ref. [51,57].
The part of the viable leptogenesis parameter region that can be accessed by colliders corresponds to the symmetry protected scenario described in section 2.2 [15,64,88]. 11 Mixing angles much larger than the estimate (2.20) require 1, and explaining the BAU with n s = 2 requires µ 1. For the collider phenomenology it is sufficient to consider the model (2.26) with ε a = µ = 0. Small deviations from this limit have a negligible impact on the production and decay rates and we use the results from ref. [57], wherein displaced vertex searches are discussed in the symmetric limit. As will be discussed in more details in section 5.3, due to the small mass splittings the heavy neutrinos can oscillate into their antiparticles and vice versa. This has no effect on the sensitivity of the considered displaced vertex searches, since we do not distinguish events from heavy neutrinos and antineutrinos anyways.
We focus on the following future lepton colliders with these specific physics programs for our study: • FCC-ee: The Future Circular Collider in the electron positron mode with its envisaged high integrated luminosity of L = 110 ab −1 for the Z pole run 12 .
• CEPC: The Circular Electron Positron Collider running at the Z pole and 240 GeV center-of-mass energy with an integrated luminosity L = 0.1 ab −1 and 5 ab −1 , respectively.
• ILC: The International Linear Collider running at the Z pole and 500 GeV center-ofmass energy with an integrated luminosity of L = 0.1 ab −1 and L = 5 ab −1 , respectively.
Lepton colliders produce heavy neutrinos primarily by the process e + e − → νN . At the Z pole, the production process is dominated by the s-channel exchange of a Z boson while at both center-of-mass energies of 240 and 500 GeV the production process is dominated by the t-channel exchange of a W boson. For its cross section σ νN we can thus take the following mixing angle dependency for the two cases: σ νN (U 2 ) at the Z pole and σ νN (U 2 e ) for the center-of-mass energies above the Z pole.
The produced heavy neutrinos decay into four different classes of final states: semileptonic (N → jj), leptonic (N → ν), hadronic (N → jjν), and invisible (N → ννν). We display the branching ratios for the four classes with varying heavy neutrino mass in figure 2. We note that for the branching ratios in the figure we summed over all the lepton flavours. We are mainly interested in semileptonic decays of the heavy neutrino, which provide a charged lepton of flavour a from which one can probe the squared mixing angle U 2 a and thus test the flavour mixing pattern. The resulting branching ratios of the semileptonic decays with a charged lepton a are given by Br(N → a jj) 0.5 × U 2 a /U 2 . In the narrow width approximation, the expected number of the displaced decay events of the heavy neutrino that decay semileptonic with a charged lepton of flavour a can be expressed as Here L denotes the integrated luminosity of the experiment, and P (x 1 , x 2 , τ) denotes the fraction of the displaced decays of the heavy neutrino with proper lifetime τ to take place between the detector-defined boundaries x 1 and x 2 . The lifetime is given by the inverse of the total decay width, which is proportional to U 2M 5 if the masses of final state particles are neglected. The probability of particle decays follows an exponential distribution such that P can be written for x 2 ≥ x 1 as with the relativistic β = v/c and Lorentz factor γ. We assume that the displaced vertex signature is free from SM background (see ref. [57]) for the boundaries as given by an SiD-like detector [116], given by the inner region (x 1 = 10 µm) and the outer radius of the tracker (x 2 = 1.22 m). The numerical calculation of the cross section for the different discussed performance parameters of the considered colliders is done in WHIZARD [117,118] by including initial state radiation and only for the ILC by including also a (L,R) initial state polarisation of (80%,20%) and beamstrahlung effects.
For the expected number of semileptonic events, N sl = a N a , we demand at least four events over the zero background hypothesis to establish a signal above the 2σ level. In the case of the Z pole run, the active-sterile mixing dependence of N sl is given by U 2 , which allows to uniquely infer U 2 for N sl exceeding four events. In the case of the center-of-mass energies above the Z pole run, however, U 2 cannot be determined uniquely since N sl depends differently on U 2 e than on U 2 µ and U 2 τ due to the production cross section being dependent on U 2 e . This ambiguity leads to the "guaranteed discovery" and "potential discovery" regions discussed at the end of section 3.3.
Moreover, the heavy neutrino massM could be measured from the invariant mass of the semileptonic final states M jj . Its precision can be assumed to be of the same order as the jet-mass reconstruction, which is ∼ 4 % for jet energies of 45 GeV with the Pandora Particle Flow Algorithm [119]. If a sizeable number of events is present, a more precise mass resolution may come from the νµ − µ + final states. For a displaced vertex the neutrino momentum can be inferred from the requirement of pointing back to the primary vertex, which yields the invariant mass M νµµ [120].
In order to determine the achievable precision on measuring the flavour pattern realised by leptogenesis we consider only the statistical uncertainties on the flavour mixing ratios U 2 a /U 2 . The observable random variables areN sl which is Poisson distributed with mean N sl , and theN a which follow a multinomial distribution with probability p a = U 2 a /U 2 . The expected number of semileptonic decays with a lepton of flavour a is given as above by N a = N sl U 2 a /U 2 . The precision on measuring U 2 a /U 2 , expressed as δ(U 2 a /U 2 ) U 2 a /U 2 with δ being the standard deviation for U 2 a /U 2 , comes from the statistical uncertainty of the ratio N a /N sl . Since N a is not independent on N sl the following precision for the flavour mixing unlike for the usual propagation of error where the uncertainties add. This point is discussed in more detail in appendix B. We confront the leptogenesis parameter space region with the achievable statistical precision of the flavour-dependent mixing U 2 a /U 2 for the different lepton colliders in section 5.2.

Sensitivity in theM − U 2 plane
In figure 3 we show the region in theM − U 2 plane consistent with leptogenesis and the potential of future collider experiments from the displaced vertex search to probe sterile (right-handed) neutrinos with this mass and active-sterile mixing.
The grey region at the top of the plots corresponds to the experimental constraint on U 2 from DELPHI [121,122]. We have not included the constraint from displaced vertex searches from LHCb (cf. ref. [62]), which are slightly more sensitive in the range between 5 and 10 GeV but do not directly probe only U 2 (it probes mainly U 2 µ ). The grey area at the bottom is excluded since in this region, assuming two right-handed neutrinos, the observed two light neutrino mass squared differences cannot be generated. The region below the blue line indicates the parameter space for which the baryon asymmetry can successfully be generated by leptogenesis.
The left column of the figure shows the case of normal ordering (NO) for the light neutrino masses, and the right column the case of inverse ordering (IO). We show the results forM < 50 GeV since above 50 GeV the uncertainty for the leptogenesis calculation increases significantly. Regarding the collider sensitivities, we show the lines for four expected events. The FCC-ee with √ s = 90 GeV (solid, green) is displayed in the top row, the ILC with √ s = 90 GeV (solid, red) and with √ s = 500 GeV (solid, brown) in the middle row, while the discovery line of the CEPC with √ s = 90 GeV (solid, purple) and with √ s = 240 GeV (solid, orange) 13 is shown in the bottom row. With the performance parameters and run times considered, the FCC-ee clearly shows the best prospects of probing the right-handed neutrinos involved in the leptogenesis mechanism. It can cover a large part of the parameter region consistent with leptogenesis.
The ILC and the CEPC both have significantly better sensitivity for the IO case, especially regarding the runs with higher center-of-mass energy. In fact, there are no expected events consistent with leptogenesis for ILC with √ s = 500 GeV and CEPC with √ s = 240 GeV in case of normal ordering. The reason is that for the NO case the electron mixing, which mediates the dominant production channel for the higher energy runs, is suppressed compared to the mixing to the other flavours through the requirement of reproducing low energy neutrino parameters.
For the ILC with √ s = 500 GeV and the CEPC with √ s = 240 GeV there are two lines shown, a dashed line and a solid line. This takes into account that the sensitivity for these 13 Since the FCC-ee also features a physics run at √ s = 240 GeV with the same integrated luminosity this result is also valid for the FCC-ee. runs depends not only on U 2 , but also on U 2 e . The solid line means that in the parameter space for consistent leptogenesis there are parameter points which can be probed by the experiment, whereas the dashed line means that for all the leptogenesis parameter points with this U 2 the right-handed neutrinos can be discovered. These two cases correspond to the "potential discovery" and "guaranteed discovery" regions discussed in section 3.3.
We note that compared to the FCC-ee, CEPC plans a much shorter run time for the 90 GeV run, since the current plans focus on Higgs measurements (and therefore on 240 GeV). A longer CEPC run time at 90 GeV could strongly improve the discovery potential for right-handed neutrinos, up to sensitivities close to the ones of the FCC-ee, a comparison plot is provided in figure 4.
Finally we remark that further above the four-event lines shown in figure 3, large numbers of displaced vertex events from the long-lived heavy neutrinos could be observed, especially at the FCC-ee, as shown in figures 9 and 10. As we will discuss in the next subsection, this large number of events can even allow for precise measurement of the flavour composition of the right-handed neutrinos.

Precision for U 2
a /U 2 in the U 2 a /U 2 − U 2 plane for different flavours In figures 5 and 6 we show the precision for measuring U 2 a /U 2 with a = e, µ, τ using the method described in section 4. Due to the potentially large number of events, the future experiments can not only discover the right-handed neutrinos but also measure their flavour-dependent mixing, i.e. U 2 a /U 2 . Together with a measurement of M and U 2 , this can be a first step towards checking the hypothesis that the observed right-handed neutrinos are indeed responsible for the generation of the baryon asymmetry of the universe (as well as for the observed light neutrino masses), as we discuss below.
The coloured regions in figures 5 and 6 correspond to the parameter space where leptogenesis is possible, takingM = 30 GeV as an example. The lines in the different colours correspond to the precision that can be achieved for measuring the ratios U 2 a /U 2 (with a = e, µ, τ ). The sensitivity depends also on the other flavour ratios not shown explicitly, and we display the most conservative precision estimate here, i.e. for the choice of the other parameters where the precision is lowest. For NO, the precision is best for measuring U 2 µ /U 2 and U 2 τ /U 2 since the active-sterile mixing with the e flavour is suppressed in the NO case. For the IO case, on the contrary, the best precision can be achieved for U 2 e /U 2 . The possible large number of events at the FCC-ee allows for precision to the percent level (cf. figure 5). At the ILC and CEPC (cf. figure 6) a precision up to about 10% -5% could be reached for part of the parameter space forM = 30 GeV. We remark that for smaller masses, a larger number of displaced vertex events could be measured, which would improve the relative precision of the flavour mixing ratios.
The plots for inverted ordering feature prominent spikes. They are a result of the fact that leptogenesis with the largest U 2 requires a flavour asymmetric washout, i.e., a strong hierarchy amongst the U 2 a . These can be understood from the fact that a large U 2 implies large Yukawa couplings. Larger Y ia increase both, the source and washout terms in the kinetic equations (3.5) and (3.6). For U 2 10 −8.5 , a complete washout of the lepton asymmetries prior to the sphaleron freezeout can only be avoided if one of the U 2 a is much smaller than U 2 , protecting the asymmetry stored in that flavour from the washout. However, the requirement to explain the observed light neutrino mixing pattern imposes constraints on the relative sizes of the U 2 a , cf. figure 7. For IO it turns out that the electron has to couple maximally U 2 e /U 2 ≈ 0.94. This explains the peak in the bottom left panel of figure 5. Having a large U 2 e /U 2 requires the other ratios to be small: U 2 µ /U 2 +U 2 τ /U 2 0.06. But still there is the freedom to choose which of the ratios is small. It turns out that the largest possible U 2 is given if the muon couples minimally, explaining both the peak in the bottom middle and the bottom right panel of figure 5. Note that consequently the maximal height of the peaks should be equal for all three flavours. For NO the electron has to couple minimally, U 2 e /U 2 ≈ 0.006, in order to allow for the largest U 2 . In this case the range in which the other ratios are allowed is rather large: U 2 µ /U 2 + U 2 τ /U 2 0.994. Consequently, the peaks are not visible in case of NO.
How to read these plots. If heavy neutral leptons are discovered at a future collider, the relative size of their mixings U 2 a can be used to test of the hypothesis that these are the common origin of light neutrino masses and baryonic matter in the universe [64,86,94,123]. With figures 5-7 we can estimate how strong this test can be at CEPC, ILC and FCC-ee. We use the caseM = 30 GeV as an example to illustrate this. If all three U 2 a could be measured exactly by experiments, then one could simply check whether the point with the observed ratio U 2 e /U 2 µ lies within the region in figure 7 that is allowed for the observed U 2 . This can either support or rule out the hypothesis that the discovered particle is involved in neutrino mass generation and leptogenesis. In reality the U 2 a can only be measured with a finite precision. This smears out the corresponding point in figure. Figures 5 and 6 can be used to estimate the expected uncertainty for any given value of the U 2 a . Further input which will affect such consistency checks will of course come from measurements of the parameters in U ν and of the light neutrino mass ordering by neutrino oscillation experiments. Note that we have completely neglected the experimental uncertainties on these parameters in the plots.

Measuring ∆M at colliders
It is important to note that, if the parameters of the right-handed neutrinos pass the necessary condition imposed by the previous test (i.e. the observed flavour-dependent ratios are in agreement with neutrino mass generation and leptogenesis), then this does not yet mean that we can confirm them as the source of the baryon asymmetry. For this, it is in particular crucial to obtain information on ∆M . In figure 8 we show the regions of ∆M which are consistent with leptogenesis and the neutrino oscillation data.
The most straightforward way of obtaining information on the mass splitting is by a direct (kinematic) measurement of ∆M phys . ∆M phys is related to the mass splitting ∆M in the electroweak unbroken phase by eq. (2.15). Realistically, however, this may be possible only for ∆M phys in the GeV range (cf. section 4). In this regime ∆M ∆M phys , so this corresponds to the largest ∆M we found to be consistent with leptogenesis in our scan.
We note that the indirect measurement in neutrinoless double β decay that was proposed in refs. [86,87,89] cannot be used in the range of M i considered here because the contribution from N i exchange to this process is strongly suppressed by their virtuality.
Another possible way to probe ∆M phys is via non-trivial total ratios between the rates of LNC and LNV involving the N i . This is possible for instance at proton-proton or electron-proton colliders, where there are unambiguous LNV signatures. Some work in this direction has been done for the LHC [69,72,[124][125][126][127]. Non-trivial ratios require a decay rate Γ N which satisfies Γ N ∼ ∆M phys . With Γ N ≈ 6.0 × 10 −6 U 2 GeV for our benchmark value ofM = 30 GeV, we obtain that e.g. U 2 = 10 −9 requires ∆M phys = O(10 −14 ), which is much smaller than ∆M θθ . Such small ∆M phys is in principle possible, but requires a large cancellation between the contributions from ∆M and ∆M θθ in eq. (2.15). From eq. (2.15) we can see that the cancellation happens for a correlation between ∆M and Re ω. For larger ∆M phys , the ratios between the rates of LNC and LNV processes gets close to one and it is not possible to infer ∆M phys .
Furthermore, a direct measurement of the heavy neutrino mass splitting ∆M phys could be possible via resolved heavy neutrino-antineutrino oscillations at colliders, as recently discussed in [126]. This also allows to measure ∆M phys when the total (integrated) ratios between the rates of LNC and LNV processes is close to one, since the heavy neutrinoantineutrino oscillations give rise to an oscillating pattern between the rates of LNC and LNV processes as a function of the vertex displacement. The oscillation time is directly related to the mass splitting ∆M phys . Again, from eq. (2.15) we get that a measurement of ∆M phys implies a non-trivial relation between ∆M and Re ω, which could be used to test leptogenesis. An interesting limit is given by ∆M ∆M θθ , or even ∆M = 0, which corresponds to the pure (or approximate) linear seesaw scenario. Then, ∆M phys ≈ ∆M θθ , with good prospects to resolve the oscillation patterns and measure ∆M phys (cf. [126]). But also for ∆M well above ∆M θθ , such that ∆M ∆M phys , the heavy neutrino-antineutrino oscillations could be resolved, for instance at FCC-hh where a large boost factor is possible which enhances the oscillation length in the laboratory frame. CEPC @ s =240 GeV CE PC @ s =9 0 Ge V Figure 3: The blue "BAU" line shows the largest possible U 2 for which the BAU can be generated for givenM , as found in the parameter scan described in section 3.3. The other coloured lines mark the parameter regions in which future lepton colliders can observe at least four expected displaced vertex events from Ni with properties that are consistent with successful leptogenesis. The solid and dashed lines correspond to the "guaranteed discovery area" and "potential discovery area" discussed in section 3.3. The grey area is disfavoured by DELPHI (on the top) and the neutrino oscillation data (at the bottom). We show no lower bound on U 2 from leptogenesis because it is lower than the constraint from neutrino oscillation data in this mass range. More details are given in the main text, cf. section 5.1.

Normal Ordering
Inverted Ordering The grey area is disfavoured by DELPHI (on the top) and the neutrino oscillation data (at the bottom). We show no lower bound on U 2 from leptogenesis because it is lower than the constraint from neutrino oscillation data in this mass range. More details are given in the main text, cf. section 5.1.

Discussion and conclusions
In this paper, we have investigated the question whether leptogenesis, as a mechanism for explaining the baryon asymmetry of the universe, can be tested at future lepton colliders.
Focusing on the minimal scenario of two right-handed neutrinos, we have estimated the allowed parameter space for successful leptogenesis in the heavy neutrino mass range between 5 and 50 GeV. We have improved previous calculations in various ways. The main improvement lies in the consistent inclusion of the lepton flavour violating source from heavy neutrino oscillations and the lepton number violating source from Higgs decays. In addition, we have included the effect of the temperature dependent Higgs expectation value and used updated values for the light neutrino oscillation parameters.
Regarding future colliders we have focused on the FCC-ee, i.e. the Future Circular Collider in the electron positron mode with its envisaged high integrated luminosity of 110 ab −1 for the Z pole run, the CEPC (Circular Electron Positron Collider) running at the Z pole and 240 GeV center-of-mass energy with an integrated luminosity of 0.1 ab −1 and 5 ab −1 , respectively, as well as the ILC (International Linear Collider) running at the Z pole and 500 GeV center-of-mass energy with an integrated luminosity of 0.1 ab −1 and 5 ab −1 , respectively.
We have confronted the parameter region where the heavy neutrinos can simultaneously explain the observed light neutrino oscillation data and the baryon asymmetry of the universe with the discovery potential for heavy neutrinos at the above-mentioned future lepton collider options, with results shown in figure 3. Future lepton colliders can be very sensitive in this mass range via displaced vertex searches. 14 We found that especially the FCC-ee can cover a substantial part of the heavy neutrino parameter space consistent with leptogenesis. A similar sensitivity could be achieved with the CEPC if more time for the Z pole run is devoted than assumed here, cf. figure 4.
Also the ILC and the CEPC (with the current plan for run times) can discover heavy neutrinos for a significant part of the parameters consistent with leptogenesis (cf. figure 3). For an inverse neutrino mass hierarchy we find that the runs with 500 GeV and 240 GeV can be competitive with the Z pole run, while for normal neutrino mass hierarchy only the Z pole runs at CEPC and ILC feature a heavy neutrino discovery potential within the leptogenesis parameter region.
Beyond the discovery of heavy neutrinos, towards testing whether they can indeed generate the baryon asymmetry, we have studied the precision at which the flavour-dependent active-sterile mixing angles can be measured. Due to the possible large number of events at the FCC-ee, measurements with a relative precision at the percent level would be possible (cf. figure 5). At the ILC and CEPC (cf. figure 6) a precision up to about 10 % -5 % could be reached for parts of the parameter space. We also provide a way to check whether the measured flavour-dependent ratios can be the cause of light neutrino masses and the BAU (cf. figure 7). Furthermore, we have studied which values of the heavy neutrino mass splitting are consistent with leptogenesis (cf. figure 8) and discussed how they  GeV (bottom, right) for the case of inverse ordering (IO) of the light neutrino masses. It turns out that besides the FCC-ee, cf. figure 5, only the four channels displayed here can be tested with precision better than 50 %. All other channels, e.g. the ones that test U 2 µ /U 2 and U 2 τ /U 2 at the ILC and CEPC are not plotted here because the precision that can be achieved is below 50 % for all values of U 2 consistent with leptogenesis. The two heavy neutrinos are assumed to be almost degenerate in mass at M = 30 GeV. Details are given in the main text, cf. section 5.2. could be measured at colliders. A detailed study of the prospects for measuring ∆M at future colliders is highly desirable.
Confronting the ratios of the flavour-dependent active-sterile mixing angles and the heavy neutrino mass splittings measured at future lepton colliders with the parameter space for successful leptogenesis can be a first step towards probing this mechanism of generating the cosmological matter-antimatter asymmetry.

A Total number of events
In this appendix we present detailed information about the number of events expected at the FCC-ee at the Z pole (figure 9), CEPC and ILC at the Z pole (figure 10) as well as CEPC and ILC at their maximal energy ( figure 11). The plots have been generated as explained in section 3.3.

B Precision of measuring flavour mixing ratios
Distribution for N sl semileptonic events and N a in a flavour a. The probability distribution function (PDF) for N sl semileptonic events with N a of them being in the flavour a is a product of a Poisson and a binomial distributions, where Br(a) = U 2 a /U 2 is the branching ratio of semileptonic states with a in the final state. The expected numbers of events are N sl = λ sl and N a = λ sl Br(a) ≡ λ a .
However, if one does not keep track of the total number of semileptonic events, and  The allowed mixings U 2 in comparison to the Lagrangian mass splittings ∆M (upper panels) and the physical mass splittings ∆M phys (lower panels) with an average mass M = 30 GeV are shown in blue for normal ordering (left panels) and inverted ordering (right panels), respectively. The red line represents the seesaw limit, below which the parameter region is excluded by neutrino oscillation data. The vertical, dashed, green lines correspond to the difference ∆M θθ between the eigenvalues of MN , cf. eq. (2.11), solely from the coupling to the Higgs field. The physical mass splitting ∆M phys acessible in experiments is related to the mass splitting ∆M in the electroweak unbroken phase and ∆M θθ by eq. (2.15). Note that leptogenesis is possible even for ∆M = 0 due to ∆M θθ = 0 during the electroweak crossover. marginalizes over N sl , the PDF reduces to a pure Poisson distribution: The variance of N a is then equal to its expectation value: Var(N a ) = N a = λ a . Precision of measuring U 2 a /U 2 . The main quantity to calculate is the expected value of Var(U 2 a /U 2 ). Here one may use the usual propagation of error, however, with the caveat that N sl is not independent of N a , which results in: Since the propagation of error assumes large event numbers, for small N sl one has to calculate the expected value and variance of Var(N a /N sl ) from the full PDF, and similarly, with the Euler constant γ E ≈ 0.58. The expected sensitivity is then given as: NO, ILC at √ s = 90 GeV IO, ILC at √ s = 90 GeV  In the large N sl limit this equation agrees with the result obtained through error propagation.
The vanishing uncertainty in the N a → N sl limit might seem concerning. However, note that this is the a-priori uncertainty. In other words, provided that the parameters really are such that N a = N sl , we do not expect any events in the other channels, i.e. the uncertainty of N a /N sl vanishes.

C Derivation of the evolution equations
The closed-time-path (CTP) formalism of non-equilibrium quantum field theory [128][129][130] offers a convenient way to derive quantum kineitc equations for the evolution of charge densities in the early universe. A detailed derivation of the rate equations (3.5, 3.6) in this formalism is given in the appendix of ref. [88]. However, the authors neglected the contribution from LNV processes to the kinetic equations, which we include in the present work. In the following we summarise the differences with the derivation given in ref. [88] if LNV processes are taken into account. An alternative derivation can e.g. be found in ref. [92].

C.1 Quantum kinetic equations for the heavy neutrinos
Nonequilibrium correlation functions -In the CTP formalism, all observables can be expressed in terms of correlation functions. For the heavy neutrinos, the relevant correlation functions are the spectral and statistical propagators Here N i are the Majorana spinors of the heavy neutrinos in the symmetric phase of the SM, i.e., the eigenvectors of M M . In the following, we suppress the spinor indices α, β, and i, j are heavy neutrino flavour indices. We can separate S + N into an equilibrium partS + N and a deviation δS N , i.e., S + N =S + N + δS N and decompose with the helicity projectors The equilibrium propagatorS N can be decomposed into functionsḡ ih in exactly the same way. In order to relate the functions g ih to on-shell quasiparticle occupation numbers, a number of approximations are necessary. We follow the same steps as in appendix B of ref. [88], with the exception that we keep terms up to first order inM 2 /T 2 . In the symmetry protected regimes with µ 1 that we consider in this work, we can still neglect terms involving ∆M in the constraint equation that determines the quasiparticle dispersion relations Ω i . In addition, we also neglect "thermal masses" and set Ω 2 i = k 2 +M 2 . Physically this corresponds to the assumption that these corrections are kinematically negligible. We emphasise that we do not neglect ∆M and the thermal masses in the kinetic equations, where they are absolutely crucial for the flavour oscillations. 15 This leads to the approximate relations which can be compared to eqs. (B.30) in ref. [88]. They allow to express all Lorentz components in terms of the flavour space matrices g 0h , (C.8) Using the on-shell approximations we find in the mass degenerate case In the comoving frame, the equilibrium distribution functions of fermions are given by On-shell kinetic equations -The evolution equation for the flavour matrices δf h can be obtained by the procedure presented in ref. [88], using eqs. (C.9, C.10). At leading order in the chemical potentials for leptons (µ ) and the Higgs (µ φ ) it reads where the indices i, j are the heavy neutrino flavours. The effective Hamiltonian H N can be decomposed into a vacuum mass term H vac N and a thermal correction H th N from forward scatterings and coupling to the temperature dependent Higgs field expectation value, such that (H N ) ij = (H vac N ) ij + (H th N ) ij . H N is responsible for the heavy neutrino flavour oscillations that occur due to the misalignment between their vacuum mass matrix M and Γ N . The thermal part H th N as well as the thermal damping rates Γ N andΓ a N can be expressed in terms of the Hermitian self-energy Σ H N and the anti-Hermitian (or "spectral") self-energy Σ A N , respectively. Up to numerical prefactors, these can be identified with the real and imaginary part of the usual retarded self-energy. We split both self energies into an equilibrium part and a deviation, Σ N =Σ N + δΣ N . Assuming that all SM degrees of freedom are in kinetic equilibrium, we can express δΣ N in terms of chemical potentials µ a and µ φ for the SM leptons and Higgs field. It is convenient to introduce the quantities for the equilibrium parts. This allows to express the effective Hamiltonian as Here we have introduced the vectork = (|k|, k 0k ) that is orthogonal to k in Minkwoski space,k · k = k ·k = 0, −k 2 = k 2 = M 2 i with M i the mass of the heavy neutrinos. We approximate the temperature dependent expectation value of the Higgs field as in ref. [88] by where z v ≈ 0.8 is the characteristic time where the Higgs expectation value starts to differ from zero. Further, Γ N is the damping rate of the deviations δf h towards equilibrium and reads The term describes the backreaction of the matter-antimatter asymmetries in the plasma of SM particles on the evolution of the N i . We use the notation of helicity-even and helicity-odd parts of the deviations from equilibrium The relativistic approximation |k|/k 0 = sign(k 0 ) that was adopted in ref. [88] allowed to express all rates in terms of the quantity γ(k) = 2gw a R k·Σ A N k 0 , which conserves lepton number. When allowing for the next non-vanishing order O(M 2 i /|k| 2 ), one finds that there are both, lepton number conserving and lepton number violating contributions to the damping rates. We refer to the lepton number conserving coefficient (which corresponds to γ(k) in ref. [88]) as γ + (k) in the following, and to the lepton number violating coefficient as γ − (k). At zeroth order in O(M 2 i /|k| 2 ), γ − (k) vanishes. These rates are given by Analogously we find for the Hermitian part Hi , (C. 25) and the term accounting for the effect of the Higgs field expectation value Evolution equations for number densities. To allow for a fast numerical exploration of the parameter space, we use momentum averaged rate equations in this work. The equilibrium number density n eq of the N i and the deviation δn hij from it which appear in the coupled system of differential equations (3.5, 3.6) are given by In order to obtain an equation in terms of n eq and δn hij , we face the usual problem of approximating an integral over a product by a product of integrals, because of the fact that the distributions f eq and δf hij appear together with another quantity that depends on the momentum k on the RHS of the kinetic equation (C.13). To describe the two momentum averages, for a rate X(k) we introduce We may use two different approximation strategies depending on the dependence on k and whether we are integrating over f eq , or δf hij to obtain the system of kinetic equations (3.5, 3.6): (i) Averaging with equilibrium weights. The expressions for the backreaction term (C. 19), as well as the washout term of the active charges appear multiplied with the equilibrium distribution functions of the right handed neutrino and charged leptons.
The backreaction and washout rates will therefore be governed by the lepton number conserving rateγ where we used the observation that one may neglect the contribution of order [f eq (k)] 2 if the rates are not infrared enhanced by a power smaller than k −2 without introducing an error of more than O(40%), which allows us to use the production rate from [109] based on refs. [110,111]. The lepton number violating rate is infrared enhanced, which means that we have to use the full expressioñ Similarly, if one assumes kinetic equilibrium for the right-handed neutrinos, i.e. that the non-equilibrium distribution of the right-handed neutrinos is proportional to the Fermi-Dirac distribution δf ≈ δn n eq f eq , we can calculate the rates , where h is either h th ± or h EV .
(ii) Evaluating the rate at the average momentum. The rate γ − requires a separate treatment in the overdamped regime. When one allows for lepton number violating processes, the equilibration of the weakly coupled state can proceed not only via the mixing with the strongly coupled state, but also through Higgs decays, which are suppressed by a factor z 2 M 2 /k 2 . This causes an enhancement in production for small momenta k, due to which the deviation from equilibrium δf may significantly deviate from the kinetic equilibrium assumption used in the momentum averaging procedure required to calculate (C.33). Therefore, as for a momentum mode which can realistically represent the production of the right-handed neutrinos we instead choose the average momentum of the heavy neutrinos, such that With these approximations we obtain the momentum independent equations (3.5) and (3.6) for the heavy neutrinos and the SM charges.

C.2 Evolution equations for the SM lepton charges
The gauge interactions among the SM fields are fast and keep them in kinetic equilibrium. This allows to describe the slowly evolving deviations of the SM fields from thermal equilibrium by chemical potentials, which can be expressed in terms of the charges introduced in section 3 through the approximate linear relations q X = a 2 R 3 µ X for massless bosons a 2 R 6 µ X for (massless) chiral fermions , (C. 37) cf. e.g. appendix A of ref. [133]. The evolution of the lepton doublets is described by the following differential equation where the first term corresponds to the washout term and the second term is the source term S a ≡ S aa , defined as It is fed by the off-diagonal correlations δn ij . When accounting for lepton number violating effects the source term is slightly modified when compared to the purely flavoured source in ref. [88] S ab (k) As mentioned previously, we evaluate γ − (k) at the average momentum |k av | ≈ 3.15T , while for γ + (k) we use the approximation from ref. [109]. In each case δf is replaced by δn. The momentum independent expression is given by S ab = 2 a R g w i =j Y * ia Y jb s=± γ s iIm(δn even ij ) + sRe(δn odd ij ) .

(C.41)
It is the combination (k −k) ·Σ A N that gives rise to a term that violates lepton number. Note that although the source term, γ − = γ |kav| − seems negligible compared to γ + = γ av + , it needs to be included as it violates the generalised lepton numberL, and the lepton charge generated this way may not be deleted by the lepton number conserving washout as pointed out in refs. [44,45].

C.3 Effects from spectator fields
It is known that spectator fields contribute to the chemical equilibration and redistribute charges during leptogenesis, such that some asymmetries are hidden from the washout. The charge densities of the leptons q and the Higgs field q φ , as well as the baryon charge density can be expressed in terms of the asymmetries ∆ a : q = A∆ , q φ = C∆ , B = D∆ , (C.42)

C.4 Determination of the transport coefficients
In this subsection we provide some details about the computation of the heavy neutrino dispersion relation and damping rates.
Thermal correction to the heavy neutrino mass -The momentum dependent thermal correction is defined by h th ± (k) = g w k 0 (k ±k) ·Σ H N , (C. 44) withΣ H N the Hermitian reduced self-energy of the heavy neutrinos. At leading order it can be computed from the Feynman diagram shown in figure 12. We neglect the thermal masses of the leptons and the Higgs running in the loop. The analytic expression for the Hermitian self-energyΣ H N within the hard-thermal-loop approximation has been derived in ref. [111]Σ where the heavy neutrino mass M i acts as a regulator to the IR enhancement.
Damping rate -In analogy to the thermal correction to the heavy neutrino masses, their equilibration rates can be defined via γ ± (k) = g w k 0 (k ±k) ·Σ A N (C. 49) and can also be computed from the Feynman diagram 12. Here γ − (k) is the momentum dependent lepton-number violating rate, while γ + (k) corresponds to the lepton number conserving rate. The lepton number conserving rate in the regime T M i is discussed in detail in refs. [109][110][111][134][135][136] In contrast to processes that contribute to γ + (k), where decays and inverse decays are subdominant with respect to 2 ↔ 2 scatterings, it has been assumed that for γ − (k) the 1 ↔ 2 decays and inverse decays dominate over the scatterings [44,45,92,134]. In preparation of this manuscript we have verified that 2 ↔ 2 scatterings give a subdominant contribution. In the Schwinger-Keldysh CTP formalism, rates for these 1 ↔ 2 processes can be derived from the self-energŷ Σ < N (k) = d 4 p (2π) 4 d 4 q (2π) 4 (2π) 4 δ 4 (q − k + p)i S < (p)i∆ < φ (q) , (C. 50) when using propagators in the zero-width limit i S < (p) → −2πδ(p 2 − m 2 )f ψ (p 0 )sign(p 0 ) p , (C.51) i∆ < φ (q) → 2πδ(q 2 − m 2 φ )f φ (q 0 )sign(q 0 ) , (C. 52) and whereˆ Σ = γ µΣ µ . Here it is crucial to use the modified dispersion relation of the Higgs field and the leptons due the their thermal masses, which give rise to different kinematic regimes in spite of the fact that the intrinsic particle masses vanish in the symmetric phase of the SM. The thermal masses from the gauge interactions are equal for the leptons and Higgs, but due to the additional contribution that m φ receives from the Higgs selfinteraction and the couplings to fermions, one finds m 2 φ > m 2 . For the following discussion it is useful to mention that the non-thermal masses of the Higgs boson and the leptons are zero since electroweak symmetry is still unbroken. Further, the thermal masses for the heavy neutrinos are suppressed by the Yukawa couplings compared to their Majorana masses M i . For that reason, the relevant masses are given by M i and the thermal masses m 2 φ , m 2 . Since low scale leptogenesis occurs at T M i , only the hierarchy m φ > M i + m can be realized, which allows Higgs quasiparticles to decay into the heavy neutrino and a lepton. 16 When making use of the delta function, the full loop integral can be expressed where I 0 (ω ± ) ≡ log(1 + e βω ± ) − log(−e βk 0 + e βω ± ) , (C.56) I 1 (ω ± ) ≡ x(log(1 + e βω ± ) − log(1 − e βω ± −βk 0 )) + Li 2 (−e βω ± ) − Li 2 (e βω ± −βk 0 ) , (C. 57) which has already been derived in ref. [111]. The integration limits of I 0,1 are determined by the kinematics of the on-shell particles in the loop, such that (cf. [137]) Further, we put the heavy neutrino on its mass shell k 2 = M 2 i . Note that we are interested in a momentum-averaged description of the coupled system of kinetic equations. The momentum averaged rate for heavy neutrinos with average massM can by approximated by γ av − ≈ 1.9 × 10 −2 × z 2M