Low-scale leptogenesis with three heavy neutrinos

Leptogenesis induced by the oscillations of GeV-scale neutrinos provides a minimal and testable explanation of the baryon asymmetry of the Universe. In this work we extend previous studies invoking only two heavy neutrinos to the case of three heavy neutrinos. We find qualitatively new behaviour as a result of lepton number violating oscillations and decays, strong flavour effects in the washout and a resonant enhancement due to matter effects. An approximate global $B - \bar L$ symmetry (representing the difference of baryon and a generalised lepton number) can protect the light neutrino masses from large radiative corrections, while simultaneously providing the ingredients for the resonant enhancement of the lepton asymmetry due to thermal contributions to the heavy neutrino dispersion relations. This mechanism is particularly efficient for large heavy neutrino mixing angles near the current experimental limits, a regime in which leptogenesis is not feasible in the minimal scenario with two heavy neutrinos. In this new parameter regime, low-scale leptogenesis is testable by the LHC and other existing experiments.


JHEP01(2019)164
All elementary fermions with the exception of neutrinos are known to exist with both chiralities, left-handed and right-handed, in the Standard Model (SM) of particle physics. Right-handed neutrinos could, if they exist, explain a number of open puzzles in particle physics as well as in cosmology, cf. e.g. [1] for an overview. Most importantly, they can generate non-zero neutrino masses m i that explain the light neutrino flavour oscillations via the type-I seesaw mechanism [2][3][4][5][6][7]. A key prediction of the seesaw mechanism is the existence of heavy neutrino mass states N i with masses M i m i and weak interactions with the SM leptons a (with a = e, µ, τ ) which are suppressed by small mixing angles θ ai . For M i below the TeV scale, the N i can be searched for experimentally. The experiments ATLAS [8][9][10], CMS [11][12][13] and LHCb [14,15] at the LHC currently perform such searches in the mass range M i > 5 GeV. For the M i below the W gauge boson mass considered in this work, the sensitivity is expected to improve significantly by using a wider range of signatures and improved triggers [16][17][18][19][20][21][22][23]. Further improvement could be achieved with additional detectors [24][25][26]. In the future, a lepton collider could offer an ideal tool to search for heavy neutrinos with masses below the W mass [27][28][29][30][31][32][33][34][35][36]. Searches at smaller masses M i < 5 GeV are preformed at the NA62 experiment [37][38][39] as well as at T2K [40], and in the future at SHiP [31,41].
Further motivation for the existence of heavy neutrinos comes from cosmology. Their Yukawa interactions F ai with the SM flavours a = e, µ, τ generally violate CP and can potentially generate a matter-antimatter asymmetry in the primordial plasma that filled the early Universe, which can be converted into a baryon asymmetry by weak sphaleron processes [42]. This mechanism is known as leptogenesis [43] and provides an attractive explanation for the baryon asymmetry of the Universe (BAU), which is believed to be the origin of baryonic matter in the Universe at present time (cf. e.g. [44] for a discussion). Leptogenesis can either be realised during the freeze-out and decay of the heavy neutrinos [43] ("freeze-out scenario") or during their production [45,46] ("freeze-in scenario"). The freeze-in scenario is particularly interesting from a phenomenological viewpoint because it is feasible for masses M i below the electroweak scale [47], which are within reach of experiments [48].
The number n of right-handed neutrinos is not constrained by theoretical arguments within the SM. However, in the context of many gauge extensions of the SM it should equal the number of SM generations (n = 3) to ensure the anomaly freedom of the theory. From an experimental viewpoint n ≥ 2 is needed to explain the two observed light neutrino mass splittings if the type-I seesaw is the sole origin of the light neutrino masses.
Most phenomenological studies of low-scale leptogenesis in the past have focused on the minimal model with n = 2. This effectively also describes neutrino mass generation and leptogenesis in the Neutrino Minimal Standard Model (νMSM) [46,49], where the third right-handed neutrino is a Dark Matter candidate, and the observational constraints on its properties [50,51] imply that it practically decouples. First estimates of the N i properties in the νMSM were made [52] shortly after the viability of the freeze-in mechanism in the minimal setup had been shown [46]. Following a number of conceptual treatments [53][54][55][56][57], JHEP01(2019)164 the parameter space was first systematically studied in refs. [47,58,59]. Following this, several authors have investigated details of the problem, such as the momentum averaging in the kinetic equations [60][61][62], the thermal production rates [61][62][63][64][65][66][67][68][69], the gradual sphaleron freeze-out [70], lepton number violating (LNV) effects in the decay and scattering rates [36,61,62,71] and from mixing [72], the dependence on the initial conditions [73] and the connection to neutrinoless double β decay experiments [74][75][76]. Recent parameter scans of the minimal n = 2 model that have implemented some of this progress can be found in refs. [36,74,[77][78][79][80] for the minimal seesaw and for its embeddings in inverse and linear seesaw models [81,82]. While this minimal model is extremely predictable and in principle fully testable [74,79], a key disadvantage is that the requirement to protect the generated asymmetries in the early Universe from washout limits the feasibility of the freeze-in leptogenesis mechanism to values of the mixing angles θ ai that are so small that it will be very challenging to produce the particles in sizeable numbers at the LHC.
The scenario with n = 3 has a much larger parameter space, which makes a phenomenological exploration more difficult. In ref. [83] it has been pointed out that this additional freedom can make leptogenesis with much larger mixing angles possible because it allows to make strong hierarchies |F ai | |F bi | amongst the Yukawa couplings consistent with light neutrino oscillation data [84]. This allows to protect the asymmetry in the flavour a from washout while the coupling |F bi | can be large enough to yield observable event rates at the LHC. The numerical analysis in ref. [83] is by now known to be incorrect because it neglected the early equilibration of one of the interaction eigenstates, see e.g. [75]. However, the physical argument can still be expected to be true. Further studies of the model with n = 3 [77,78,83,[85][86][87] have not systematically explored the parameter space, so that the range of heavy neutrino couplings that can be made consistent with leptogenesis and with light neutrino oscillation data in this scenario is not yet known. With the present work we want to address this issue and systematically scan the parameter space of the low-scale seesaw model with three right-handed neutrinos. While we perform an agnostic scan of the entire parameter space, we pay special attention on the region where the seesaw model approximately respects a generalised B − L symmetry [88]. In this parameter region the symmetry protects the light neutrino masses in a way that observable N i production rates at colliders can be made consistent with the observed neutrino oscillation data in a technically natural way [89].
Our systematic analysis demonstrates significant quantitative and qualitative differences with respect to the n = 2 scenario. On the one hand, we confirm that the parameter space which simultaneously accounts for the neutrino oscillation data and the observed baryon asymmetry of the Universe, projected onto experimentally accessible quantities such as the active-sterile mixings and neutrinoless double beta decay effective mass, is significantly enlarged, implying significant discovery space for experiments such as NA62, T2K, Belle II and the LHC. On the other hand, we find qualitatively new dynamical processes in the kinetic equations describing leptogenesis, such as a dynamically generated resonant enhancement, providing new channels to generate the baryon asymmetry of the Universe.

JHEP01(2019)164
The remainder of the paper is organised as follows. In section 2 we introduce our setup, with a particular focus on the role of (approximate) global symmetries. The kinetic equations governing leptogenesis are introduced in section 3, emphasising the subtleties associated with a temperature dependent mass eigenbasis. We come back to this point in section 4, where we discuss the different physical processes involved in the generation of the lepton asymmetry both for n = 2 and n = 3. The details of our parameter scan are given in section 5, with the resulting experimental prospects discussed in section 6. We illustrate the different dynamical processes contributing to leptogenesis by means of some representative benchmark points in section 7 before concluding in section 8. Further technical details can be found in the three appendices.

Review of the model and notation
The most general renormalizable Lagrangian that contains only SM fields and n flavours of right-handed neutrinos ν Ri reads Here we have suppressed SU(2) indices; ε is the totally antisymmetric SU(2) tensor. The F ai are Yukawa couplings between the ν Ri and the SM leptons a , M M is a Majorana mass matrix for the singlet fields ν Ri . 1 In the following we work in the flavour basis where M M is diagonal unless a different basis is explicitly specified. The breaking of electroweak symmetry by the Higgs expectation value φ = (0, v) T (with v = 174 GeV at T = 0) generates a Dirac mass term ν L m D ν R with m D = vF from the Yukawa interaction term F L εφ * ν R .
After electroweak symmetry breaking (EWSB), the complete neutrino mass term reads Here we have added the 1-loop correction δm 1loop ν [90] since we aim to perform an analysis that is consistent at second order in the Yukawa couplings F . The mass matrix (2.2) can be diagonalised as In the parameterisation (2.3) M is first block-diagonalised by a complex 3 × n matrix θ that mediates the mixing between the active neutrinos ν L and the sterile neutrinos ν R . The unitary matrices U ν and U * N then diagonalise the 3 × 3 and n × n blocks m ν and M N in the upper left and lower right corners, respectively, as we can recover the well-known tree-level result The spectrum of neutrino mass states is clearly separated into three light and n heavy mass eigenstates which can be expressed in terms of the Majorana spinors

JHEP01(2019)164
The mixing matrix Θ quantifies the misalignment of the mass eigenstates ν i and N i with the original "active" and "sterile" neutrinos ν L and ν R . It leads to a θ-suppressed weak interaction of the heavy mass eigenstates N i , The first two terms represent the θ-suppressed interactions of the N i via the weak currents. Through these interactions the heavy neutrinos N i replace ordinary neutrinos in all processes if this is kinematically allowed, but with amplitudes suppressed by the angles Θ ai . The third term represents the Yukawa coupling to the physical Higgs field h in the unitary gauge. Here we have employed the relation m W = 1 2 gv involving the weak gauge coupling constant g. It is convenient to introduce the quantities which practically govern the event rates for processes involving the N i .

Approximate lepton number conservation
In absence of any special structure in the matrices F and M M , the seesaw relation (2.13) suggests that which would imply unobservably tiny branching ratios in collider experiments. However, the seesaw relation is a matrix valued equation, and the light neutrino mass squares m 2 i are the eigenvalues of the matrix m † ν m ν . If there are cancellations in m † ν m ν , then small m 2 i can be made consistent with arbitrarily large U 2 ai . Constraints from experiments other than neutrino oscillation ones are comparably weak in most of the mass range between the kaon and W boson masses, cf. e.g. [84,97] and section 5.2, so that U 2 ai ∼ 10 −5 are phenomenologically viable. The cancellations could either occur accidentally (which would require a tuning of at least five orders of magnitude to achieve U 2 ai near the current LHC bounds [13]) or be owed to a symmetry. Indeed, if the Lagrangian (2.1) approximately respects a generalised B −L symmetry [88,89,98], then the cancellations in m † ν m ν occur in a technically natural way because the light neutrino masses must be proportional to small parameters that quantify the amount of B −L violation. Here B denotes the usual SM baryon number andL is a generalised lepton number, that is composed of the SM lepton number L and some charge associated with the ν R (see below). Specific models that realise an approximate B−L symmetry include models with Rparity violation [99][100][101][102][103][104], "inverse seesaw" type scenarios [105][106][107][108][109] (cf. also [110,111]), the "linear seesaw" [112,113] (cf. also [114][115][116]), scale invariant models [86], some technicolorinspired models [117,118], the νMSM [88] and other low-scale seesaw realisations [119][120][121].
Low-scale leptogenesis in connection to an approximate B −L symmetry has previously JHEP01(2019)164 been studied in the framework of linear and inverse seesaw scenarios in refs. [81,82], while the νMSM parameter space has been studied in refs. [47,52,58,59] and numerous follow up publications (cf. references given in section 1). The B −L symmetry enforces that the ν Ri must either i) decouple, ii) have vanishing Majorana masses or iii) arrange themselves in pairs that form (pseudo-)Dirac spinors. For the n = 3 case this can be made explicit with the parameterisation In the limit a , a , µ, µ → 0 the quantity B −L is conserved. In terms of the original ν Ri , one can make the following assignment of charges, spinorL-charge where the subscripts s, w indicate that the corresponding states ν Rs,w are strongly/weakly coupled to the SM states in the high temperature limit, T M . We can now write the Lagrangian in the form where we have introduced the (pseudo-)Dirac spinor ψ N = 1 √ 2 (ν Rs + ν c Rw ) and we sum over multiple occurrences of the flavour index "a".
The generalised lepton numberL is significantly violated by the oscillations amongst the heavy neutrinos even if a , a , µ, µ 1. This means that it is in general not a suitable quantity to describe the time evolution in the early Universe, at least not if the rate of the oscillations is faster than the expansion of the Universe (T < (M 0 |M 2 i −M 2 j |) 1/3 ) and faster than the equilibration time scale of the heavy neutrinos (max(|( , where γ av ∼ 10 −2 is a generic numerical coefficient appearing in the rate (3.6) and M 0 ≡ m P (45/(4π 3 g * )) 1/2 = T 2 /H can be interpreted as the comoving temperature in a radiation dominated Universe with Hubble parameter H, m P = 1.22 × 10 19 GeV and g * counting the relativistic degrees of freedom in the thermal bath. In the regime T M i the helicity states of the heavy Majorana neutrinos N i practically act as JHEP01(2019)164 "particle" and "antiparticle" states. One can use this fact to define another generalised lepton number L = L + L N , (2.22) that is approximately conserved even if the parameters a , a , µ, µ are not small. Here L N is a quantum number that can be assigned to the helicity states P h N i , where P h = 1/2 × [1 + hγ 0 γ i γ 5 (p i /|p|)] is the helicity projector with momentum p, as 2 spinors L-charge It is indeed L (not L orL, both of which are violated) that is usually being referred to when distinguishing between "lepton number conserving" and "lepton number violating" terms in the low-scale leptogenesis literature.

The LHC testable scenario
The νMSM realises the B −L conservation by choosing the seesaw scale Λ below the electroweak scale and keeping all parameters a , a , µ, µ tiny. In that model, a , µ 1 is required to make sizeable F a consistent with light neutrino oscillation data and for successful leptogenesis (which requires µ 1 as only two heavy neutrinos participate in the process). Regarding the third heavy neutrino, which serves as a Dark Matter candidate, a 1 is required to ensure its longevity for any mass allowed by structure formation considerations, while µ 1 is in addition needed for consistency with indirect Dark Matter searches if one assumes a to be sizeable enough that the Dark Matter is produced thermally via weak interactions, cf. e.g. [50,51]. The parameter space of the νMSM is the subject of many past and ongoing studies and will not be further addressed here because it is, from the viewpoint of neutrino mass generation and leptogenesis, practically a scenario with n = 2 due to the extreme smallness of the a , µ required for the stability of the Dark Matter candidate. Instead we focus on scenarios where all three heavy neutrinos have masses M i of roughly the same magnitude. In this context it is worthwhile noting that it is sufficient for the B −L conservation that either µ = 0 or a = 0 (as well as a = 0 = µ), since in this case the third right-handed neutrino either decouples or obtains only a Dirac mass term. Hence, scenarios with moderately small a and µ of order unity are technically natural. It turns out that the choice a , a , µ 1 with µ 1 allows for successful leptogenesis with |F ai | larger than the electron Yukawa coupling. This leads to mixings U 2 a1 and U 2 a2 that are well within reach of current experiments.

Quantum kinetic equations
The quantum kinetic equations for freeze-in leptogenesis [45] in the density matrix formalism [122] read (see e.g. refs. [46,47,60]): where the n×n matrix R N encodes the density matrix of the three heavy neutrinos in kinetic equilibrium normalised to the entropy density, The SM sector is taken to be in thermal equilibrium, and is thus fully characterised by the chemical potentials where µ La are the flavoured left-chiral lepton chemical potentials and µ φ is the Higgs chemical potential, which appear in the corresponding distribution functions. These are connected to the chemical potential associated with B − L, µ ∆ = a µ ∆a by the relation while the terms involving coefficients γ (i) and γ (i) represent L-conserving and L-violating dissipation rates, respectively. Note that we neglect the L-violating part of H . N D = 2 is an SU(2) factor. The thermally averaged rates
The numerical values of c are determined following ref. [36]. 3 Both c The running of the SM gauge couplings (included in the L-conserving rates) is given by where g 0 = 0.652 and g 0 = 0.357 denote the values of the corresponding gauge couplings at Λ = πT = m Z . For the purpose of our numerical scan, we find it convenient to switch the time variable from cosmic time t to x = T EW /T , leading to the system of differential equations (A.9)-(A.11). See ref. [67] for further details.
3 What was computed in ref. [36] is in fact the quantity γ (1b) − γ (1a) that appears in front of the term µF F † in eq. (3.2) after rewriting We extracted γ (1b) and γ (1a) from this by assuming that these coefficients have equal values. Practically this means that we guessed the coefficient in front of the term µ(F (RN − 1)F † + F * (RN − 1)F T ). We do not expect this to have any phenomenological consequences because that term is small at all times: at early times µ is small, and at late times the heavy neutrinos are close to equilibrium.

JHEP01(2019)164
The final B−L asymmetry is obtained by evaluating the chemical potentials µ ∆ at the scale of electroweak symmetry breaking, where g s = 106.75 denotes the effective number of degrees of freedom in the SM above the EW phase transition. SM sphaleron processes only pick up the asymmetry in the active sector, converting it to the baryon asymmetry we observe in the Universe today, Here n B,B counts the number density of (anti-) baryons today, s denotes the entropy density of the Universe and the baryon asymmetry of the Universe is measured to be Y B = (8.6 ± 0.01) × 10 −11 [123]. For later reference we introduce also the asymmetry in the heavy neutrino sector, which in the absence of L-violating processes is identical to Y B−L . Note that the definition (3.16) refers to quasiparticle occupation numbers and should therefore be applied in the basis where the effective Hamiltonian H is diagonal, which does not necessarily coincide with the flavour basis where M M or M N are diagonal, as discussed in the following.

Mass and interaction bases
It is worthwhile to emphasise that some caution should be taken with respect to the flavour basis in which the above equations are defined. In general, neither the basis where M M nor the one where M N is diagonal correspond to the physical (quasi)particle mass basis. On the one hand there is the O[θ 2 ] correction in eq. (2.9) from the Higgs expectation value which affects the physical masses at temperatures below ∼ 160 GeV [124]. In the B−L symmetry protected regime the physical mass splitting at T = 0 (given by the eigenvalues of M N ) can considerably deviate from the splitting between the eigenvalues of M M . This effect, which was already pointed out in ref. [52], can be crucial for the generation of late time asymmetries (and hence DM production [125]) in the νMSM [47]. Recent discussions of possible phenomenological implications can be found in e.g. refs. [36,79], where it is also described how the v(T ) dependent term should be added to the effective Hamiltonian (3.5). Moreover, the mixing induced by the temperature dependent Higgs field value v(T ) can have a significant impact on the generation of lepton asymmetries shortly before sphaleron freeze-out [72]. We ignore both of these effects in the following because we expect that they only lead to O [1] corrections in a limited part of the parameter region that we consider.
On the other hand there are corrections to the dispersion relations in a medium from forward scatterings [126][127][128], known as thermal masses or matter potentials, that affect the properties of SM particles [129] and heavy neutrinos [130]. These are responsible for the term ∝ F † F T /8 in the effective Hamiltonian (3.5). The physical heavy neutrino quasiparticles in the primordial plasma correspond to the eigenstates of the full Hamiltonian (3.5),

JHEP01(2019)164
including the thermal correction. Since the relative size of the thermal and vacuum masses changes with temperature, the physical quasiparticle mass basis rotates throughout the evolution of the Universe. At high temperatures, when the thermal masses dominate, it should be identified with the interaction basis where F † F is diagonal. At low temperatures it coincides with the vacuum mass basis where M N is diagonal. It is important to note that the interpretation of the diagonal elements of the density matrix R N as measuring the corresponding number densities only holds if the density matrix is expressed in the quasiparticle mass basis.
To which degree the rotation of the effective mass basis affects the generation of asymmetries depends on the model parameters. It is only relevant if the heavy neutrinos have reached sizeable occupation numbers while the temperature dependent contribution from the "thermal masses" to the splitting of the eigenvalues in the effective Hamiltonian (3.5) still dominates over the vacuum splittings M 2 i − M 2 j . This happens quite generically for experimentally accessible heavy neutrinos because the approximate B −L symmetry that is required to make sizeable F ai consistent with light neutrino oscillation data implies that at least two M i are quasi-degenerate. We focus on this case in the following. For small mixing angles finite temperature effects usually only lead to sub-dominant corrections because the thermal masses are smaller and the mass splittings are in general larger.
At high temperatures the L-conserving rates γ (i) are parametrically larger than the L-violating rates γ (i) , and one can understand the behaviour in terms of the eigenvalues and eigenvectors of the matrices F † F and M † M M M . In the B−L conserving limit the matrix F † F has three vastly different eigenvalues; one is roughly given by a |F a | 2 while the other two are suppressed by 2 a and 2 a . Hence, one interaction eigenstate (ν Rs , which is always part of the pseudo-Dirac spinor ψ N ) feels the full strengths ∼ F a of the Yukawa interactions, while the other two (ν Rw and ν R3 or combinations thereof) practically decouple.
Unless µ is very close to unity (i.e. M 3 M ), the mixing and rotation mainly occur between the components of the pseudo-Dirac spinor ψ N because the thermal corrections to . We shall consider this case first and neglect ν R3 for the moment. In this case most of the discussion for the case n = 2 in ref. [75] can directly be applied. At high temperatures the states ν Rs and ν Rw defined in eq. (2.20) are approximately both, effective mass and interaction eigenstates. ν Rs picks up a large thermal mass ∼ a |F a | 2 T 2 and is produced at a large rate ∼ a |F a | 2 γ (i) , while the corrections to the mass of ν Rw are only ∼ a | a F a | 2 T 2 and its production rate is only ∼ a | a F a | 2 γ (i) . Despite the large thermal mass splitting, there are no rapid oscillations between the two states because the effective heavy neutrino mass and interaction bases are almost aligned. If ν Rs comes into equilibrium before the oscillations commence, then the BAU is generated in a single overdamped oscillation in the strong washout regime ("overdamped regime"). This is in contrast to the case of small mixing angles, where a large number of oscillations occur before the sphaleron freeze-out in the weak washout regime ("oscillatory regime"), see ref. [75] for details. At late times the feebly coupled state is driven to equilibrium by the overdamped oscillation and by the L-violating damping rates, which are ∝ F * F T and not a -suppressed.

JHEP01(2019)164
If one considers all three heavy neutrinos, then the situation can be much more complicated. For n = 2 there are practically only two relevant time scales in the heavy neutrino sector: the time when the first heavy neutrino state reaches equilibrium (given by its thermal damping rate) and the frequency of the oscillations (given by the mass splitting). On the contrary, there are generally three equilibration and three oscillation time scales in the system with n = 3. In the simplest scenarios all frequencies and damping rates are well-separated, and no heavy neutrinos reach thermal equilibrium before the sphaleron freeze-out. In that case the separation of scales allows to treat each of the oscillations separately in a simplified n = 2 model, similar to the treatment of light neutrino oscillations in the Sun or in the atmosphere. However, for µ 1 and µ ∼ 1 all three states can mix with each other and complicated behaviour can arise. We illustrate this for a few example points in section 7.
Finally we note that the charge L in eq. (2.22) should be defined in the rotating quasiparticle mass basis. The relation between this basis and the vacuum mass basis, however, depends not only on temperature, but also on the model parameters. For simplicity we use in the following the vacuum mass basis and the interaction basis as approximations for the actual quasiparticle mass basis at very low and very high temperatures, respectively.

Most important new physical effects
The generation of a baryon asymmetry via the freeze-in mechanism is a complex nonequilibrium process that involves an interplay between coherent oscillations and decoherent scatterings, both of which can occur in a L-conserving or L-violating manner. For n = 3 there is a large number of (possibly vastly different) times scales involved, and the phenomenology of the leptogenesis parameter space is very rich. While our scan systematically explores this parameter space numerically, a qualitative analytic understanding of the results is highly desirable.

Minimal n = 2 scenario without L-violation
Let us first briefly review the case of only two heavy neutrinos, in order to highlight the qualitative differences in the full three neutrino case studied in this paper. Moreover, we neglect for the moment L-violating processes. In this case the system can be studied by analytic methods [45,46,75,81]. No CP-violation can arise in the heavy neutrino oscillations [75], so that the CP-violation necessary to generate a lepton asymmetry must arise in the active sector and/or in the mixing between the active and sterile sectors. In particular, washout processes play a crucial role in the generation of a net L = 0. In the weak washout regime 4 an analytical expression for the lepton asymmetry was derived in 4 For n = 2, one can parametrically distinguish a weak washout regime (or "oscillatory regime") in which the equilibration of both heavy neutrinos occurs after the freeze-out of weak sphaleron processes (i.e. |(RN )ii|, |(RN )ii| 1 for all i and T > T sph ) and a strong washout regime (or "overdamped regime") in which the occupation numbers of one heavy neutrino interaction eigenstate reach equilibrium before sphaleron freeze-out. For weak washout the BAU scales as ∝ (m 2 [46], while the dependence in the strong washout regime is rather complicated [75].

JHEP01(2019)164
refs. [45,46] (see also [81]). It was found to be proportional to i,a (F ai ) † δ a F ai with which in particular vanishes in the flavour symmetric limit In the strong washout regime an asymmetric coupling f 1 to the active flavours is typically necessary to protect the asymmetry from the strong washout processes in the sterile sector, see e.g. the discussion in refs. [75,82], where The max and min are to be understood with respect to the index a, and the term between parenthesis averages the sum with a weight proportional to the relative size of the Yukawa couplings for the right-handed neutrino i. The approximate second equality holds only in the B −L conserving limit, where F a are the large entries in the parameterisation (2.19).
If one active flavour is coupled significantly weaker than the other generations, the asymmetry in this flavour can be preserved for a considerable time even when the right-handed neutrinos approach thermal equilibrium. For n = 2 the heavy neutrino mixing pattern is strongly constrained by light neutrino oscillation data [39,74,79,131,132], and in particular f > 5 × 10 −3 [39]. This imposes an upper limit on the maximal U 2 i for which leptogenesis is feasible: leaving aside highly fine-tuned parameter choices, the two heavy neutrinos necessarily form a pseudo-Dirac spinor ψ N if their mixings are much larger than the estimate (2.17). The larger Yukawa couplings F a a F a in (2.21) then drive the entire heavy neutrino sector towards equilibrium in an overdamped manner [75], and the only way to protect the BAU from washout is a strong hierarchy f 1. How strong this hierarchy must be to ensure the survival of some asymmetry until sphaleron freeze-out depends on the magnitude of the Yukawa couplings and masses (and hence U 2 i ). As a result, the experimental constraint on f from neutrino oscillation data imposes an upper limit on U 2 i for a given M i .

n = 3 scenario without L-violation
The situation is qualitatively different in the case of three right-handed neutrinos. One can distinguish three different new physical effects: 1) Larger flavour hierarchies are allowed. A hierarchy in the couplings of the heavy neutrinos to individual SM flavours (f 1) can protect a part of the lepton asymmetry from the washout even if the heavy neutrinos reach equilibrium (R N RN ∼ 1) if f is small enough to keep one of the washout rates below the Hubble rate. In the scenario with n = 2 the requirement to reproduce the light neutrino oscillation data practically requires f > 5 × 10 −3 [39] (cf. also [74,79]). This hierarchy is not strong enough to protect JHEP01(2019)164 the asymmetry from washout [36] for mixings U 2 i near the current LHC bounds [13]. As already pointed out in ref. [83], the relaxed lower experimental bound on f in the scenario with n = 3 [84,133] allows to protect the BAU from washout for much larger U 2 i if one SM flavour couples only very feebly to the pseudo-Dirac pair.
2) Asymmetry in the heavy neutrino oscillations. Contrary to the n = 2 case discussed above, the n = 3 case allows for a generation of a net asymmetry L = 0 during the heavy neutrino oscillations (even if L-violating effects are neglected), without requiring any flavour asymmetric Yukawa couplings. In appendix B we explicitly derive the corresponding source terms entering the quantum kinetic equations by means of a perturbative expansion in the lepton asymmetries. We emphasise the presence of a new source term for the asymmetry in the heavy neutrino sector, which arises (for n ≥ 3 only) from the first term in eq. (3.1). This allows for the generation of an asymmetry even in the absence of (flavour asymmetric) washout processes, contrary to the situation for n = 2 [75].
3) Resonantly enhanced asymmetry. The produced asymmetry strongly depends on the heavy neutrino mass splitting and is resonantly enhanced if the splitting between two of the heavy neutrino masses is very small [46]. In the primordial plasma the effective quasiparticle masses are given by the eigenvalues of the effective Hamiltonian (3.5). Due to the interplay between temperature dependent and independent terms in the effective Hamiltonian, the effective mass splittings are time dependent. As a result, a maximal resonant enhancement can be generated dynamically, even if the mass spectrum in vacuum is only moderately degenerate. 5 This is similar to the Mikheyev-Smirnov-Wolfenstein (MSW) effect that affects light neutrino oscillations in matter. In contrast to the MSW effect for light neutrinos it does not require the presence of lepton chemical potentials because the Yukawa couplings can give different thermal masses to the neutrinos (while the light neutrinos' gauge interactions are flavour blind, so that different effective masses can only be realised through chemical potentials). In particular, an (avoided) level crossing necessarily occurs for µ > 1, i.e., if the state ν R3 with couplings F a3 ∼ a F a has a vacuum mass M 3 >M larger than the pseudo-Dirac spinor ψ N with couplings ∼ F a . This is because the component ν Rs of ψ N defined in (2.20) receives a comparably large thermal mass ∼ |F a | 2 T 2 , which necessarily exceeds the effective mass of ν R3 at sufficiently high temperature. If this crossing occurs during the time when the asymmetry is generated, the resonant effect can maximally enhance it, even if the vacuum masses are only moderately degenerate. In contrast, in the B −L protected regime of the n = 2 case, the interaction and Majorana mass bases have to be maximally misaligned to reproduce the small active neutrino masses, and hence any avoided level crossing comes with a mass gap which is typically too large to resonantly enhance the asymmetry. For n = 3 with µ > 1 the level-crossing temperature T crossing can be JHEP01(2019)164 estimated in the limit of approximate B −L symmetry (| a |, | a |, µ 1 in eq. (2.19)), yielding

Effects of L-violation
In the weak washout or oscillatory regime, the BAU is generated at temperatures T M i . In this case L is in good approximation conserved during the heavy neutrino oscillations. However, for U 2 i that are large enough to be probed with the LHC, the couplings F a are generally large enough to drive part of the heavy neutrinos (in particular ν Rs ) to equilibrium before sphaleron freeze-out. In this strong washout or overdamped regime, the asymmetry is often generated near the sphaleron freeze-out, and L violating effects can in general not be neglected. This leads to several new effects. . This is in particular important in the case n = 2, where the oscillations themselves cannot generate a L = 0 (cf. point 2) ) and leptogenesis always relies on the flavour asymmetric washout, so that the BAU is necessarily of order O[F 6 ai ]. It has been pointed out in refs. [71,135] ai ] for certain parameter choices. For n = 3 we expect this effect to be less relevant because a L = 0 can already be produced in the oscillation in absence of L-violation, cf. appendix B. This latter effect is ∼ O[F 6 ai ] [75], but not suffering from any M i /T -suppression, it can be active over a much longer period of time. We should stress that the above power counting only holds when the washout is weak enough that there is a clear separation between the times when the heavy neutrinos start to oscillate and when they come into equilibrium ("oscillatory regime"). If some heavy neutrino degrees of freedom reach equilibrium at early times ("overdamped regime"), the parametric dependence is different [75].

5)
Damping of ν Rw . In absence of L-violating processes the heavy neutrino damping rates are approximately proportional to F † F . In the B −L symmetric limit, one eigenvalue of this matrix is much larger (∼ F 2 a ) than the other two (∼ 2 a F 2 a , 2 a F 2 a ). This means that the states ν Rw and ν R3 approach thermal equilibrium very slowly. In particular, ν Rw is primarily driven to equilibrium via the overdamped oscillation with ν Rs [75]. The heavy neutrino damping rate due to L-violating processes is proportional to ∼ F T F * . This in particular means that ν Rw is driven to equilibrium at a rate ∝ F 2 a (M /T ) 2 , which can be much larger than the L-conserving rate ∝ 2 a F 2 a , since the weakly coupled eigenstate of F T F * in general does not coincide with ν Rw . In contrast to that, the state ν R3 , which is not part of the pseudo-Dirac system ψ N , remains feebly coupled with equilibration rate ∝ 2 a F 2 a (unless M 3 M , in which case there can be significant mixing between ν R3 and the other two states).

JHEP01(2019)164
6) Washout efficiency. If the L-violating processes are neglected, the washout of the helicity-based L charges in the heavy neutrino sector enforces a simultaneous decrease of the L charges in the active sector, suppressing the SM lepton number asymmetry L, i.e. Y B−L . In the presence of L-violating processes this is not necessarily the case, allowing for a sizeable asymmetry in the active sector even if the washout in the sterile sector is effective. This results in a larger final BAU [36].

7)
Equalising flavoured asymmetries. When some of the heavy neutrinos come into equilibrium, the L-conserving processes wash out the flavoured asymmetries L a , but cannot erase the total SM asymmetry L (and hence the BAU) unless all charges in the heavy sector have been erased. In this situation the washout of the total L is driven by the L-violating processes. Since L is in equal parts composed of e, µ and τ asymmetries, the direction in flavour space in which all L a are equal is only slowly erased, while deviations from L e = L µ = L τ are erased by the much faster L-conserving processes. Therefore the asymmetries in all SM flavours tend to be equal in this situation, cf. e.g. figure 4. Note that the total L would also be erased in absence of the L-violating processes once the sterile charges are driven to equilibrium because the total L vanishes for our initial conditions. 8) Direct L-generation through active-sterile mixing. So far we have mainly considered Lviolating decays in the symmetric phase of the SM (in particular Higgs decays). Similar arguments apply if one includes L-violating scatterings in the symmetric phase of the SM. There is, however, moreover a brief period between the moment when the Higgs field develops a non-zero value v(T ) at T ∼ 160 GeV and the sphaleron freeze-out at T ∼ 130 GeV. During this period the mixing between active and sterile neutrinos directly violates L [72]. We neglect this effect in the present work.
In summary, the dynamics governing leptogenesis in the n = 3 case is much more diverse than in the n = 2 case. After quantifying the parameter space yielding successful leptogenesis in sections 5 and 6, we will illustrate some of the effects described above by means of a few exemplary points in section 7.
5 Strategy for the parameter scan

General strategy
Our goal is to perform a systematic parameter scan to identify the range of the heavy neutrino properties that are consistent with all experimental constraints and can reproduce the observed BAU. A major obstacle is the high dimensionality of the parameter space. For n heavy neutrinos, the seesaw model contains 7n−3 free parameters in addition to those of the SM. Only 5 of those (two mass splittings and three mixing angles) are constrained by light neutrino oscillation data. For n = 2 it is possible to perform a complete parameter scan to clearly identify the boundaries of the region where leptogenesis is possible [36,47,59,79] or perform a Bayesian analysis [74]. For n = 3 the dimensionality of the parameter space is so high that even a systematic combination of all experimental constraints in the mass

JHEP01(2019)164
region under consideration here (without leptogenesis) is numerically challenging [84]. If one includes the computation of the BAU, which requires solving the coupled differential equations (3.1) and (3.2) for each parameter choice and is numerically much more demanding than imposing laboratory constraints, then it becomes practically impossible to explore the entire parameter space. However, from a phenomenological viewpoint, one is mostly interested in the projection of the viable parameter region on the M i − U 2 ai planes. In section 6 we present scatter plots in these planes which illustrate that, for masses M i below the electroweak scale, the leptogenesis region covers the entire experimentally allowed range of mixings U 2 ai . Since both, the BAU and experimental constraints, depend smoothly on the model parameters that determine the U 2 ai , it seems physically reasonable that the entire region between the scattered points can be filled if the scan would run for infinitely long. In the remainder of this section we explain how the parameter scan is performed.
It is well-known that leptogenesis is feasible in the n = 2 model with values of U 2 i at most four orders of magnitude above the estimate (2.17), i.e. U 2 i < 10 −6 GeV/M i , for any value of M i in the range considered here [79]. Since the n = 2 parameter space is a subset of the larger parameter space of the n = 3 model under consideration here (in the limit a → 0), the same must apply in the present model. We are therefore primarily interested in studying leptogenesis with U 2 i > 10 −6 GeV/M i . The strongest constraints on the N i properties come from the seesaw relation (2.13) and light neutrino oscillation data. The requirements to reproduce the observed data without fine-tuning for U 2 i > 10 −6 GeV/M i practically enforces the B−L symmetry discussed in section 2.2. The parameterisation (2.19) in principle is ideal to explore this region, but the preproduction of the light neutrino parameters in a randomised scan is a search for the needle in the haystack due to the small error bars of these parameters and the complicated relations between model parameters and observables. To keep the numerical effort at a feasible level, we adopt a three-step strategy that treats neutrino oscillation data different from other constraints: 1. For the generation of parameter points, we employ the Casas-Ibarra parameterisation [136], see below. This parameterisation is not ideal to explore the B−L symmetry protected region, but guarantees consistency with light neutrino oscillation data at the perturbative level.
2. We then remove all points which are not consistent with the experimental constraints described in section 5.2.
3. For each remaining parameter choice we compute the BAU. We consider leptogenesis feasible if the BAU deviates from the observed value by less than a factor five. 6 To improve the numerical performance we rewrite the quantum kinetic equations (3.1) 6 The experimental uncertainty on the BAU is much smaller than this. The larger range for an acceptable YB adopted here reflects theoretical uncertainties as well as the strong sensitivity of the final value of YB on the model parameters -starting from a parameter point whose computed YB deviates from the observed value by an O(1) factor we expect to be able to reproduce the observed value by minimally varying the input parameters.

JHEP01(2019)164
and (3.2) as described in appendix A. We furthermore only track the first ten oscillations of the heavy neutrinos and then set off-diagonal elements of the density matrix to zero. We explicitly check that this does not induce any discontinuity in the evolution of the asymmetries. A similar procedure has been proposed in ref. [58] and has been analytically verified in ref. [57].
Parametrisation. In order to reproduce the observed neutrino masses and lepton mixing we adopt a generalisation of the Casas-Ibarra parameterisation [136], extended to include 1-loop corrections (2.12) to the light neutrino mass matrix [96]. For small θ we may approximate the relation (2.10) by 7 i.e., M can be approximately expressed in terms of the entries of M diag N . In this formalism the Yukawa couplings are determined after having specified the low-energy neutrino oscillation data, the right-handed neutrino masses and a 3-dimensional orthogonal matrix R, The matrix R can be parameterised as a product of three rotations, where the angles ω ij are in general complex numbers and and analogous definitions hold for V 13 and V 23 . We work in the basis where the charged lepton Yukawa couplings are diagonal, so that U ν can be identified with the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix U PMNS .
The full system is thus characterised by 13 real free parameters: 6 real numbers for the three imaginary angles ω ij , 3 heavy neutrino masses, 3 complex phases in the PMNS mixing matrix (one Dirac δ CP and two Majorana α 1,2 ) and the overall light neutrino mass scale. On the other hand, neutrino oscillation data fix (within experimental uncertainties) the mixing angles in the PMNS matrix and the mass differences in the light neutrino spectrum, although current data does not allow to disentangle between two possibilities for the ordering of light neutrino masses (normal ordering (NO) and inverted ordering (IO)). 7 If the splitting between two eigenvalues of MM is smaller than the light neutrino mass differences, then the O[θ 2 ] term in eq. (2.9) can cause large deviations of UN from unity. However, in the B −L symmetric regime one still observes |(F UN )ai| 2 |Fai| 2 and hence |θai| 2 |Θai| 2 = U 2 ai due to the structure of the F and MM in eq. (2.19).

JHEP01(2019)164
We note however that global fits of neutrino oscillation data currently show a preference for NO at 3σ [137,138], and current experiments are starting to constrain the value of the Dirac phase δ CP [139][140][141]. In our scan we randomly generate Yukawa couplings F accordingly to the relation (5.2), using the ranges of input parameters reported in table 1. For the heavy neutrino mass splittings and the complex angles in R, we randomly alternate between drawing our parameters from a linear versus a logarithmic distribution. This enables us to effectively sample the different regions of the parameter space, including the B −L protected regime as well as different flavour and mixing structures. The PMNS mixing angles and light neutrino mass splittings (as well as the PMNS Dirac phase δ CP ) are allowed to vary in the 3σ ranges as determined by the NuFIT collaboration [142].
Targeted scans. In order to efficiently explore the most interesting regions of the parameter space in a reasonable timescale, we run three different scans: • generic scan: all the generated points complying with neutrino constraints and reproducing the observed BAU value are collected; • large mixing targeted scan: the BAU value is only computed for points featuring a mixing U 2 > 10 −4 (GeV/M 2 ) 2 . This region is especially interesting because it can be probed by Belle II, LHCb, ATLAS and CMS; • low mass region targeted scan: we only generate solutions where the lightest of the heavy neutrino masses (M 1 ) is lighter than 0.35 GeV. This region can be probed in the decay of kaons, for instance by T2K.

Further experimental constraints
The realisations of the seesaw mechanism constructed as outlined in section 5.1 (by construction compatible with the neutrino oscillation data) are compared to the following experimental constraints (cf. e.g. ref. [84] for a more detailed discussion): • The sum of light neutrino masses must remain below i m i < 0.12 eV, as determined by the Planck collaboration [143].
• We impose constraints on the deviation from unitarity of the PMNS mixing matrix, as determined in [84,149]. • We require an upper bound on the lifetime, requiring the sterile neutrinos to not be too long-lived in order to not spoil predictions from Big Bang Nucleosynthesis (BBN). We set a conservative upper bound of 0.1 seconds, cf. ref. [172].

Theoretical considerations: parameter volume and tuning
In addition to the experimental constraints listed above we apply a number of theoretical arguments.
• Perturbativity. Although the parameterisation in eq. (5.2) allows for an efficient exploration of the parameter space, the complex angles ω ij cannot acquire arbitrary values: the magnitude of the Yukawa couplings F grows exponentially with the modulus of the imaginary parts Im ω ij , and too large couplings are excluded, either because they break the perturbative expansion leading to eq. (2.8), thus rendering the full parameterisation unreliable, or because they give rise to a strong dynamics. In our scan we allow for sizeable values of the imaginary angles, and explicitly diagonalise the full 6 × 6 mass matrix (including 1-loop corrections) in order to verify the

JHEP01(2019)164
agreement with experimental data, excluding realisations that do not comply with them; moreover we require each entry in |F | to be smaller than 4π.
• Fine-tuning. In the exploration of the parameter space we do not impose any symmetry, but we allow the underlying parameters in the theory to vary as reported in table 1, in order to generate symmetry protected as well as generic solutions. We then quantify a posteriori the level of fine-tuning for each solution, by defining the following quantity

Results
In this section we discuss the results obtained performing the parameter scan described in section 5. Projecting the high-dimensional data set consisting of all parameter points meeting the experimental constraints (including the requirement of successful leptogenesis) on to different physically meaningful two-dimensional planes, we illustrate the qualitative new features arising in the n = 3 case of " freeze-in leptogenesis". Figure 1 depicts the allowed range of active-sterile mixing after imposing all experimental constraints as a function of the heavy neutrino mass. We find that large mixing angles U 2 ai right up to the current experimental bounds are allowed in the n = 3 case across the entire mass range we consider. This is in contrast to the model with n = 2, where a gap of one order of magnitude was reported in [36] for M i ∼ 5 GeV that grows to about three orders of magnitude for M i ∼ 50 GeV. Moreover, we find points with very low fine-tuning (according to the criterion of eq. (5.5)) in the entire viable parameter space projected on to the mass-mixing plane. This provides rich prospects for ongoing and planned experiments searching for sterile neutrinos in the GeV range. In particular, in contrast to the n = 2 scenario, searches for prompt decays of N i at the LHC [9,13,15] and Belle II [162,179] can probe the viable leptogenesis parameter space for n = 3. Moreover, in the region of large mixings and for M i below ∼ 20 GeV, displaced vertex searches at the LHC [9,17,18,22,23,180,181] could see thousands of events, assuming that displacements in the mm range can be resolved. This would allow for a determination of the heavy neutrino flavour mixing pattern [35,36], which is crucial to test the hypothesis that JHEP01(2019)164 these particles are responsible for leptogenesis [74,79]. For n = 2 the sensitivity of such searches could barely touch the viable leptogenesis parameter space [48], and it seems unlikely that the flavour mixing pattern can be measured at a level that allows to draw any conclusions. Hence, the perspectives to test low-scale leptogenesis are much better in the scenario with n = 3.

The range of allowed mass and mixing
A few comments on the distribution of the points in the scatter plots in figure 1 are in place. The main purpose of these plots is to illustrate that leptogenesis is feasible in the entire mass-mixing plane without fine-tuning in the sense of eq. (5.5). The density of points within the allowed area should not be misinterpreted as a measure for any theoretical or experimental preference for particular values. Instead, it is primarily a result of the parameterisation (5.2) and the randomisation procedure described in section 5.1. In particular, some of the most prominent features in the distribution of points appear because we performed a number of targeted scans as described on page 20. In addition to the variation in the density of points within the allowed region, there are also parts of the mass-mixing planes that appear to be empty. This does not necessarily imply that there are no viable parameter choices in these regions, but may also simply indicate that our scan failed to fully exploit these regions. For instance, the distribution of points above M i = 2 GeV suggests that leptogenesis is feasible for mixings all the way up to the experimental upper limit on the individual U 2 ai , but not all the way up to the experimental upper limit on the total U 2 i . We suspect that the reason is that it is difficult to find points where all three mixings U 2 ai are maximal for one of the N i within the parameterisation (5.2), while there is no reason why such points would not exist. Similarly, it is very difficult to explore the region of large U 2 τ i below M i = 2 GeV, while there is evidence that, at least from the point of view of neutrino mass generation, this region is experimentally allowed for n = 3 [84]. This is in contrast to the n = 2 model, where the results presented in ref. [79] indicate that this region is indeed ruled out by the combination of different constraints. Finally, a similar problem arises in the determination of the lower bound on the mixings. While the light neutrino oscillation data and the requirement for the N i to decay before BBN both impose lower bounds on the U 2 i that depend on m lightest [84,133], neither of them can impose a lower bound on the individual U 2 ai for n = 3. The BBN constraint can always be avoided if the N i decays into a SM final state of different flavour, while the neutrino oscillation data can always be explained if another heavy neutrino provides the required mixing with the flavour a.
Relative mass degeneracy. Figure 2 shows the parameter points of figure 1 projected onto a plane spanned by the two mass splittings among the heavy neutrinos. As discussed above, the density of the points carries little physical meaning. It is however remarkable that we find viable leptogenesis points in the entire parameter plane, for all possible hierarchies of the heavy neutrino masses, and covering a wide range of values for the physical mass differences. The regions along the top and right axes correspond to a situation with one pair of very degenerate neutrinos and a third neutrino with at least an O(1) hierarchy. This can be realised in the B−L symmetry protected regime for µ 1 and µ ∼ 1 or physically equivalent configurations in which the labels of the N i are permutated. This   Figure 1. Active-sterile mixing for the viable BAU solutions as a function of the heavy neutrino mass, for a normal (left) and an inverted (right) ordering in the light neutrino mass spectrum. From top to bottom: electron U 2 ei , muon U 2 µi , tau U 2 τ i and summed U 2 i mixings. The grey region is excluded by direct searches of heavy neutral leptons (cf. section 5.2), the lines show the expected sensitivities for the ongoing experiments T2K [182], NA62 [39], Belle II [183], LHCb [180] with an integrated luminosity of 380 fb −1 , and for ATLAS and CMS with an integrated luminosity of 300 fb −1 . The latter include different proposed searches: [22] (continuous line), [17] (dashed line), [21] (dotted line). region contains effective n = 2 models if the third neutrino decouples. In the upper right corner both, µ and µ , are sizeable, and there is no protecting symmetry for the neutrino masses. On the other hand, the central and bottom left area of figure 2 is characterised by three very degenerate neutrinos, with in general all three of them contributing to leptogenesis. The low value of the fine-tuning (according to the criterion of eq. (5.5)) indicates that again this is a B−L protected region. Finally, in the top right corner we find the fully non-degenerate N i spectra, which can accommodate leptogenesis only at the cost of fine-tuning.

JHEP01(2019)164
A further important difference between the n = 2 and the n = 3 case appears in the flavour structure, i.e. in the relative coupling strength of a given N i to the 3 active flavours, For n = 2, the requirements of successful neutrino mass generation and leptogenesis limit the allowed values of these ratios, see refs. [36,74,79]. On the contrary, for n = 3 we find parameter points yielding successful neutrino mass generation and leptogenesis in the entire parameter space. This provides a further interesting possibility to test these leptogenesis mechanisms.

Effect on neutrinoless double β decay
It is well known that the exchange of N i with masses below the electroweak scale can make a significant contribution to the rate of neutrinoless double β decay [111,[145][146][147] in the region where freeze-in leptogenesis is feasible [74,76,78,97,148]. The decay rate is proportional to the quantity where the first term comes from light neutrino exchange and the second one from N i exchange. Here

JHEP01(2019)164
where p is the momentum exchange in the decay and depends on the isotope, cf. e.g. [147]. In our analysis we use the numerical value p = 125 MeV, resulting from an average over different decaying nuclei (see e.g. [146]). Using the estimate (2.17) one would generically expect that the relative size of the two contributions is roughly given by f A (M i ). Since we found viable parameter points for which U 2 i exceeds the estimate (2.17) by several orders of magnitude, one may wonder whether leptogenesis with large mixing angles generally predicts that m ββ greatly exceeds the standard contribution, from light neutrino exchange. However, large U 2 i can only be achieved without fine-tuning if the light neutrino masses m i are protected by the B −L symmetry. This symmetry automatically suppresses m ββ and sets the rate of neutrinoless double β decay (as well as the light neutrino masses) to zero if the symmetry is exact. It is instructive to study how much "tuning" is required to obtain a large decay rate if the symmetry violating parameters are not exactly zero. To see this explicitly we bring eq. (6.1) into the form by using the unitarity relation 3)). Further using eq. (2.7) and the fact that Θ θ in the B −L conserving regime, we can recast this as The first term in this equation is always smaller than the standard prediction, so large contributions can at most come from the second term. 8 The contribution from N 3 is proportional to θ 2 e3 = (v/M ) 2 × F 2 e ( e /µ ) 2 . A priori this term looks potentially large in low-scale seesaw models because of the factor (v/M ) 2 and because of the second power of µ in the denominator, which threatens to cancel the suppression from the e in the numerator. However, there are also at least two powers of µ in the numerator, one from e v 2 /(µ M ) that N 3 makes to the neutrino masses. The current upper limit on the sum of neutrino masses [143] is comparable to the limit on m ββ [184], hence a large contribution to the decay rate could only be achieved at the cost of cancellations in m ν and/or m ββ that are not explained by the B −L symmetry. For µ 1 the contributions 8 It is straightforward to show that all eigenvalues of mν vanish if the parameters a, a and µ in eq. (2.19) are set to zero, which implies that also m ν ββ exactly vanishes in this limit.

JHEP01(2019)164
from N 1 and N 2 are individually large because of the much larger mixing angle. However, due to the B −L symmetry, they interfere destructively, which is manifest in the imaginary unit i in eq. (2.19). To estimate their contribution, we expand to linear order in the B −L violating parameters This may again be compared to the contribution ∝ 2F 2 e (−2 e + µ)v 2 /M that N 1 and N 2 make to neutrino masses through their mixing with ν Le . The term ∝ e F 2 e in m ββ is always smaller than its counterpart in m ν forM < v, while the term ∝ µF 2 e is parametrically of comparable size. Hence, for generic choices of the parameters that are not dictated by the symmetry, the current neutrino oscillation data clearly disfavours large contributions to m ββ from the heavy neutrinos. Using the expression (6.3) for m ν ββ and the numerical values of the mixing (U ν ) ei , the same argument suggests that the contribution from the N i exchange to m ββ is comparable or smaller than that from light neutrinos. Based on this, one can estimate that m ββ should not greatly exceed the standard prediction |m ν ββ | unless the model parameters are either highly tuned to cause accidental cancellations amongst the contributions involving different SM flavours in m ν , or there exist additional flavour structures/symmetries that lead to such cancellations. There is, however, one way to avoid this conclusion that has already been discussed for n = 2 in refs. [74,76,78]: large deviations from m ν ββ can be obtained in a technically natural way if the M i are of the same order as p. This results from the combination of two factors. On the one hand the contribution of heavy neutrinos is maximal if their masses are comparable with the exchanged virtual momentum p, cf. eqs. (6.1), (6.2); on the other hand loop corrections to the light neutrino parameters are proportional to the heavy neutrino masses, cf. eq. (2.7). This is confirmed by the results shown in figure 3. The plot confirms the claim from ref. [78] that leptogenesis in the n = 3 low-scale seesaw model is compatible with both, a rate of neutrinoless double β decay that is much larger or much smaller than the standard prediction |m ν ββ |. However, a much larger rate tends to require a considerable tuning in the sense of eq. (5.5). The lower panels in figure 3 confirm that sizeable contribution from the heavy neutrinos to m ββ can be achieved together with a low fine-tuning (in the sense of eq. (5.5)) for masses of order O(100 MeV). • Both, larger and smaller mixings U 2 i can be made consistent with baryogenesis and neutrino mass generation, cf. figure 1. In the entire mass range studied here, the upper limit on U 2 i is practically given by experimental constraints. That is, for any JHEP01(2019)164 value of U 2 i that is allowed by experiments, one can find a set of model parameters for which baryogenesis is feasible. This considerably improves the perspectives for current and planned experiments to test the mechanism of baryogenesis. At the same time, there is no lower bound on the individual U 2 i . This is in contrast to the case with n = 2, where the estimate (2.17) practically acts as a "floor" for experimental searches.

Benchmark points
• The constraints on the heavy neutrino mass spectrum are relaxed. In particular, no mass degeneracy is needed to generate the BAU.
• The constraints on the flavour mixing pattern, i.e., the relative size of the heavy neutrino couplings to different SM flavours, are relaxed.
Some of these effects have been predicted in the past. For instance, in ref. [83] it was argued that the relaxed constraints on the flavour mixing parameter f should allow for baryogenesis with much larger U 2 i than for n = 2. The fact that baryogenesis with n = 3 is feasible for

JHEP01(2019)164
non-degenerate heavy neutrino spectra was discussed in detail in ref. [85]. Different aspects of the L-violation have been discussed in refs. [36,61,62,71,72,80]. However, it turns out that the behaviour for n = 3 in general is much richer than anticipated in these works. In general, the evolution of charges is governed by a complex interplay of several amongst the effects 1) -8) in section 4, and one cannot uniquely relate the viability of a particular parameter choice to any individual of these mechanisms. It is nevertheless instructive to illustrate some of the most important physical effects for a few selected benchmark points. The parameters of these points are summarised in table 2.
Benchmark point I): resonant enhancement due to level crossing. (Figures 4, 5) The first parameter point we consider is given by the choice This point is similar to the first one in the sense that there is also an approximate B −L symmetry and the masses of all three heavy neutrinos are quite degenerate. The fine-tuning is somewhat higher, but with f.t.(m ν ) = 3.2 × 10 −2 still small. It also leads to overdamped behaviour and exhibits two avoided level crossings between the state that corresponds to ν Rs at high T and the two other states. 9 The first one resonantly enhances the asymmetry production, while the second one is much weaker and occurs when L-violating processes are already relevant and two of the heavy neutrinos have reached equilibrium. The main difference, however, lies in the strongly hierarchical flavour structure, f = 6.2 × 10 −3 . This prevents the equalising of all SM flavours by effect 7) because the muon flavour couples only very feebly to the other charges. This helps to avoid washout in spite of the fact that two heavy neutrino degrees of freedom reach equilibrium around x 0.02 as a result of effect 5) by the L-violating processes. Both level crossings lead to a re-distribution of charges, which is visible in the left panel of figure 6, but only the second one leads to a sign change in the sterile charges. It is worthwhile noting that the zero crossing of the total sterile charge caused by the second level crossing does not enforce a zero crossing of the total active charge L due to effect 6). we checked that this has no significant effect on the asymmetries because they "average out" at later times.
Benchmark point III): large mass splittings. (Figures 8, 9) The final point that we consider is given by   Loop corrections remain comparably small, f.t.(m ν ) = 0.14, in spite of the fact that the parameters µ and µ are not small. There is a moderate flavour hierarchy f = 9.5 × 10 −2 . The evolution of charges corresponds to the standard mild washout scenario. The point serves as an example that leptogenesis can be realised with O(1) mass splitting without resorting to extreme fine-tuning.

Conclusions
The ARS mechanism [45] for "freeze-in leptogenesis" is a remarkable and testable idea to implement leptogenesis within a minimal extension of the Standard Model by adding heavy neutrinos with masses below the electroweak scale. In this paper we perform the first systematic investigation of the ARS mechanism with three right-handed neutrinos (n = 3), extending previous analyses which encompassed only two right-handed neutrinos actively participating in leptogenesis (n = 2). For n = 2 there are only two characterstic time scales associated with the heavy neutrinos -the oscillation period of the two neutrinos, set by their mass difference, and the thermalisation rate, set by their coupling to the SM. On the contrary, for n = 3, a much richer phenomenology arises. As we show in this work, this does not only enlarge the parameter space, enhancing the possibility of a detection in present collider experiments due to a large mixing with the SM neutrinos (as anticipated in [83]), but moreover we find qualitatively new mechanisms to generate the lepton asymmetry, which do not have a counterpart in the n = 2 analysis. The most striking of these qualitatively new effects is a resonant generation of a lepton asymmetry associated with an (avoided) level crossing of the effective mass eigenvalues of the three heavy neutrinos. As is well known from the analysis of the n = 2 case, the generation of a lepton asymmetry is enhanced for a small mass splitting within the neutrino pair. In the case of three neutrinos, a tiny mass splitting can occur dynamically through thermal corrections to the mass eigenstates, which induce a level crossing in the eigenvalues of the effective Hamiltonian. If this occurs when the respective heavy neutrinos have already JHEP01(2019)164 been produced in significant numbers, but have not yet reached full equilibrium, then the lepton asymmetry is resonantly enhanced. This enables successful leptogenesis with only a mild degeneracy in the vacuum masses of the heavy neutrinos and without any fine-tuning in the flavour mixing pattern.
Moreover, compared to the n = 2 case, we find richer flavour structures leading to successful leptogenesis. With only two right-handed neutrinos, successful leptogenesis with mixing angles that may be accessible with the LHC prefers a hierarchical flavour structure, where the generated asymmetry can be protected from washout when stored in a very weakly coupled SM flavour, whereas the requirement to reproduce the observed neutrino oscillation data sets an upper bound on the flavour hierarchy. This tension forbids large mixings between the heavy and the SM neutrinos, making experimental tests challenging. In the case of three heavy neutrinos, these constraints are relaxed in two ways. Firstly, the additional parameter freedom due to the additional state allows to comply with the neutrino oscillation data while simultaneously allowing for a large flavour hierarchy. Secondly, we demonstrate that contrary to the n = 2 case, a lepton asymmetry can be generated even with flavour democratic couplings, due to a new term in the kinetic equations which only arises for n ≥ 3. Consequently, and again contrary to the n = 2 case, we find large mixing between the heavy and the SM neutrinos, right up to the current bounds, to be compatible with both neutrino oscillation data and successful leptogenesis.
We point out that this large mixing, as well as the formation of pseudo-Dirac pairs of right-handed neutrinos, is natural in the context of an approximate global B −L symmetry, where B denotes the SM baryon number andL denotes a generalised lepton number under which also the right-handed neutrinos are charged. With this in mind, we define 'fine-tuned' solutions as parameter points for which the radiative one-loop contributions to the light neutrino masses are large compared to the tree-level contributions. In this sense, we find that experimentally accessible large mixing is possible without any fine-tuning, whereas an enhancement of the neutrinoless double β decay rate is possible only at the cost of fine-tuning unless the heavy neutrino masses are rather close the the momentum exchange in the process. Furthermore, the resonant generation of a lepton asymmetry due to a level crossing of the mass eigenvalues occurs quite generically in the regime protected by the B −L symmetry, since one state of the pseudo-Dirac pair receives large thermal corrections whereas the quasi decoupled third right-handed neutrino does not. Moreover, the participation of the quasi decoupled heavy neutrino automatically protects the generated asymmetry from subsequent washout. All this renders the B −L symmetry protected regime particularly interesting for ARS leptogenesis.
At high temperatures far above the heavy neutrino mass scale, the different helicities of the right-handed neutrinos are conserved quantum numbers. Approaching the EW phase transition, this approximation breaks down, allowing for 'lepton number violating' ( Lviolating) processes. Although active for only a fairly short period of time, these processes can significantly alter the predicted lepton asymmetry. We highlight the different physical processes at work, showing that they can both enhance or reduce the final asymmetry.
In summary, we find that leptogenesis invoking the oscillations of three right-handed neutrinos just before the EW phase transition comes with some qualitative and quantitative differences to the well-studied n = 2 case. New channels of leptogenesis lead to an enhanced JHEP01(2019)164 lepton asymmetry. The viable parameter space, reproducing both the observed neutrino oscillation data and the baryon asymmetry of the Universe, projected onto the mass versus active-sterile mixing plane shows promising opportunities for ongoing experiments, such as NA62, T2K, Belle II and the LHC. This calls for a more detailed study of some of the effects that we have neglected here. These include effects of the electroweak transition (temperature dependent Higgs field value, gradual sphaleron freeze-out, particle masses generated by the Higgs mechanism and the L-violation due to the active-sterile mixing), the full momentum dependence of the equations, a fully systematic perturbative computation of the L-violating rates, as well as, a verification of the validity of the gradient expansion (which justifies the usage of the density matrix equations) during the level crossing.

A Notation for the quantum kinetic equations
We provide in this appendix additional details on the system of differential equations used in our numerical scan to compute the baryon asymmetry, starting from the original set of Boltzmann equations given in eqs. (3.1) and (3.2), see ref. [82] for more details.
Let us consider first the equation for the abundances of the heavy neutrinos, described by R N and RN . Their oscillation processes are described by the term proportional to H , R N,N where the Hamiltonian H can be split as H = H 0 + V N with V N denoting the effective potential and H 0 denoting the vacuum Hamiltonian, where ∆M 2 ij = M 2 j − M 2 i , and in the last expression we have dropped the part of the matrix proportional to the unity matrix, since it drops out in the commutator of eq. (3.1). Thermal averaging yields

JHEP01(2019)164
For the numerical solution of the system of Boltzmann equations it is convenient to adopt x = T EW /T as new time variable. The change of variables is described by the relation It will be convenient to introduce a new parameterisation of the effective potential, and analogously, The analogous terms for the equations for RN are obtained by setting F → F * and µ → −µ.
Performing the change of variables t → x also in the equations for the chemical potentials µ ∆a , we finally obtain the system:

JHEP01(2019)164
where the functions φ (i) are related to the thermally averaged rates γ (i) by: (A.14) The coefficients c (i) X are given in eq. (3.11). The corresponding expressions for φ (i) are obtained by replacing γ (i) by γ (i) .
The two terms [R N , W N ] and [R N , r] in eq. (A.9) (and the corresponding terms in the equation for RN ) represent the terms originating from the potential V N and the vacuum Hamiltonian H 0 , respectively, with where T L,i denotes the typical leptogenesis temperatures associated to the oscillations between the "1st" and "ith" (i = 2, 3) heavy neutrino eigenstate.

B Perturbative expansion
In this appendix we perform a perturbative expansion of the system of Boltzmann equations in terms of the chemical potentials in the active sector µ a , following the procedure outlined in [82] for n = 2. This allows us to gain an analytical understanding of some of the main processes involved in the n = 3 ARS leptogenesis and to identify the qualitative differences with respect to the case of a single pair of quasi mass-degenerate neutrinos. In particular, we will discover an additional (leading order) source term for the lepton asymmetry, enabling successful leptogenesis in the absence of flavour asymmetric Yukawa couplings (see point 2) on page 15). While we do not employ this formalism for the main parameter scan of this paper, we have confirmed for a range of parameter points that it accurately reproduces the results of the full equations. For simplicity and in order to facilitate the comparison with the results of [82], we will omit in this appendix the L-violating terms.

B.1 0th order in the chemical potential
To leading order in µ a eq. (3.1) reads,

JHEP01(2019)164
with H = H 0 + V N introduced in appendix A. Performing unitary rotations of this equation we will in the following identify the different physical effects involved in the generation of a lepton asymmetry. 10 The first step consists in defining an 'oscillation' basis, in which the neutrino oscillations driven by the vacuum Hamiltonian H 0 are removed. This is done by performing a rotation of the form R N = E † R N E with E(t) defined as [58] with the typical oscillation temperatures encoded in the parameters r i , see eq. (A.15). With this, and eq. (B.1) can be written as with W N = E † W N E and φ (0) introduced in eq. (A.12). It will be convenient to introduce a third basis, which we refer to as 'interaction' basis, in which W N is diagonal. This is accomplished by means of the unitary matrix U , Since E and U are unitary, the eigenvalues y i of W N are proportional to those of F † F , and in particular time-independent and real. Motivated by this we construct the time independent part U c of U as To switch between the oscillation and flavour basis we introduce where we note that the matrix D is anti-hermitian and time-independent. Denoting the leading order density matrix of the right-handed neutrinos in the interaction basis by S 0 N , as in the case of two right-handed neutrinos. 10 These 'basis changes' obtained by unitary rotations of the density matrix should not be confused with the different basis discussed in section 2, which are obtained by rotating the spinors (νi, Ni).

JHEP01(2019)164
To obtain the corresponding equation for SN we need to replace the Yukawa coupling by its complex conjugate, F → F * . Denoting the quantities in the SN equation with overbars, this implies Here the first equality follows since Y contains the real eigenvalues of The second is trivial since no powers of F are involved in the definition of E, and the third follows from Finally the fourth equality in eq. (B.10) follows fromD =Ū † (dŪ /dx) withŪ = E † U * c . Note that in the case of two right-handed neutrinos the matrix D is symmetric 11 and hencē D = D T = D. For n > 2, D is anti-hermitian but not symmetric, so this simplification does not apply. With this, the equation for the opposite helicity (N ) neutrinos reads (B.14) Defining An explicit computation shows that the quantity D ∼ iU † diag(0, 1) U has purely imaginary, symmetric off-diagonal elements if and only if ϕ1 = ϕ2. Since there is a free phase in each column of U , this condition can always be met. In the case of 3 right-handed neutrinos, this freedom of choosing the phases of the columns is not sufficient to make D symmetric for a generic 3 × 3 unitary matrix U .

JHEP01(2019)164
right-hand side of eq. (B.18) is absent, and ∆S 0 − = 0 is a solution to eq. (B.18). This reflects that for appropriate initial conditions, the reduced number of CP -violating phases for n = 2 impedes the generation of asymmetries in the sterile sector (see also appendix D of ref. [75]). On the contrary, in the case of three right-handed neutrinos this is no longer the case, leading to ∆S 0 − = 0 already at leading order. Secondly, in the case of two righthanded neutrinos, the last term in eq. (B.17) is absent. One might be tempted to discard this term, since it is proportional to a (small) asymmetry, however at early times when the oscillations are large, the off-diagonal terms of ∆S 0 − can in fact be rather large. We note that in particular in the case of (mildly) hierarchical Yukawa couplings this term can be crucial to obtain the correct thermalisation time scales of the different right-handed neutrino species.
B.2 1st order in the chemical potentials Sterile sector. In the oscillation basis, eq. (3.1) reads (B.19) where as above X = E † XE. Using eq. (A.4), as well as the functions φ (i) and the o µ ,ō µ defined in appendix A, this becomes where we have dropped the subleading term proportional to µ∆S − in the φ (1b) term and As indicated above, in the context of our perturbative expansion, the leading order term driving the asymmetry in the sterile sector is the first term in the second line in eq. (B.24), which is present already at 0th order but is absent for n = 2.
To good approximation, we may set S + 2S 0 in the first term of the second line of eq. (B.24). In this approximation, the equation of motion for ∆S + decouples, and the equations of motion describing the sterile sector are (B.17) and (B.24). In the case of only two sterile neutrinos, this is in fact an exact result to first order in µ a . For completeness, we give here also the equations for ∆S + : Note that contrary to eq. (B.24) (see also discussion below eq. (B.18)), the equation for ∆S + has no source term in the limit ∆S − , µ a → 0, justifying the approximation S + 2S 0 above.
Active sector. The starting point for the equation of the active sector is eq. (3.2). Replacing the γ (i) with φ (i) and using eq. (A.4), this can be rephrased as Both these expressions are hermitian matrices, implying that the diagonal components are real, i.e. we obtain

(B.33)
Note that the asymmetry in the sterile sector sources an asymmetry in the active sector through S aux . More precisely, in the absence of L-violating terms, the final asymmetries in the active and sterile sectors are of equal magnitude but opposite sign.
In summary, all processes relevant for ARS leptogenesis involving three right-handed neutrinos are well described by the system of differential equations (B.17), (B.24) and (B.33). The results obtained from this simplified system agree up to percent-level with the results obtained by solving the original system (3.1) and (3.2) in the absence of LNV processes.

C Approximate analytical solution describing the level crossing
In this appendix we present the approximate solution to the leading order right-handed neutrino number density evolution, eq. (B.1), in the B −L symmetric limit. We introduce the notation for the effective Hamiltonian term and the production term, respectively. After approximately diagonalising the pseudo-Dirac block, the equilibration matrix takes the form: For µ |µ 2 − 1|, the effective Hamiltonian term can be approximated by: Assuming that the off-diagonal correlations are either oscillating quickly, or overdamped, their mean value approaches: with k = l = i = j. In principle, the Yukawa couplings F a can be large enough to cause early equilibration of the sterile neutrinos. In that case, the diagonals of the density matrix R N are approximately given by: where we have neglected the equilibration through mixing of the pseudo-Dirac pair. Note that the equations are given in the interaction basis, where the subscript s corresponds to the strongly coupled state ν Rs and w to the weakly coupled one ν Rw . The number density of the heaviest right-handed neutrino is governed by the equation: dR N 33 dx = − (R N 33 − 1) Γ 33 + Γ 11 |∆H 13 | 2 − |Γ 13 | 2 (∆H 11 − ∆H 33 ) 2 + (Γ 11 /2) 2 + (C.6) + (∆H 11 − ∆H 33 ) Re (∆H 13 Γ * 13 ) (∆H 11 − ∆H 33 ) 2 + (Γ 11 /2) 2 .
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.