Three-Flavoured Non-Resonant Leptogenesis at Intermediate Scales

Leptogenesis can successfully explain the matter-antimatter asymmetry via out-of-equilibrium decays of heavy Majorana neutrinos in the early Universe. In this article, we focus on non-resonant thermal leptogenesis and the possibility of lowering its scale. In order to do so, we calculate the lepton asymmetry produced from the decays of one and two heavy Majorana neutrinos using three-flavoured density matrix equations in an exhaustive exploration of the model parameter space. We find regions of the parameter space where thermal leptogenesis is viable at intermediate scales, $T\sim 10^{6}$ GeV. However, the viability of thermal leptogenesis at such scales requires a certain degree of cancellation between the tree and one-loop level contribution to the light neutrino mass matrix and we quantify such fine-tuning.


I. INTRODUCTION
There is overwhelming experimental evidence for an excess of matter over antimatter in the Universe. This asymmetry remains a fundamental and unresolved mystery whose explanation demands new physics beyond the Standard Model (SM). The baryon asymmetry may be parametrised by the baryon-to-photon ratio, η B , which is defined to be where n B , n B and n γ are the number densities of baryons, anti-baryons and photons, respectively. This quantity can be measured using two independent methods that probe the Universe at different stages of its evolution. Big-Bang nucleosynthesis, BBN, [1] and Cosmic Microwave Background radiation, CMB, data [2] give η B BBN = (5.80 − 6.60) × 10 −10 , η B CMB = (6.02 − 6.18) × 10 −10 , at 95% CL, respectively. As the uncertainties of the CMB measurement are smaller than those from BBN, we shall apply the CMB value throughout this work. In order to dynamically produce the observed baryon asymmetry in the early Universe, the mechanism of interest must satisfy the Sakharov conditions [3] 1 : B (or L) violation; C/CP violation and a departure from thermal equilibrium. Leptogenesis [4] satisfies these conditions and produces a lepton asymmetry which is subsequently partially converted to a baryon asymmetry via B + L violating sphaleron processes [5].
Leptogenesis is particularly appealing as it typically takes place in models of neutrino masses, simultaneously explaining the baryon asymmetry and the smallness of the neutrino masses. In its simplest realisation, the lepton asymmetry is generated via out-of-equilibrium decays of heavy Majorana neutrinos [6][7][8][9]. This process occurs approximately when the temperature, T , of the Universe equals the mass scale of the decaying heavy Majorana neutrino.
In general, the scale of thermal leptogenesis is not explored below the Davidson-Ibarra (DI) bound, M 1 ≈ 10 9 GeV [10]. Davidson and Ibarra found an upper bound, proportional to M 1 , on the absolute value of the CP-asymmetry of the decays of the lightest heavy Majorana neutrino. This constrains the regions of parameter space in which successful leptogenesis may occur as a function of M 1 . This translates into the DI bound on M 1 itself as the minimum value required for successful leptogenesis. There have been a number of in-depth numerical studies which support this bound and require M 1 ≥ 10 9 GeV [11,12] in conjunction with a bound on the lightest neutrino mass, m 1 ≤ 0.1 eV [11,13,14].
The original derivation of this bound makes some simplifying analytical assumptions and hence is subject to three caveats: only the lightest heavy Majorana neutrino decays; the heavy Majorana neutrino mass spectrum is hierarchical; and flavour effects, which account for the differing interaction rates of the charged-lepton decay products of the heavy Majorana neutrinos, are ignored. In this work, we shall investigate scenarios of three-flavoured thermal leptogenesis in a more general setting than these conditions allow. We shall then consider lower heavy Majorana neutrino masses at scales M 1 ≈ 10 6 GeV. Given the existence of low-scale lep-2 togenesis models at the TeV scale, we shall refer to this as "intermediate" scale leptogenesis.
There are several reasons to explore leptogenesis at intermediate scales. Firstly, the introduction of heavy neutrinos to the SM leads to a correction to the Higgs mass which may potentially be unnaturally large. This is because the correction to the electroweak parameter µ 2 (the negative of the coefficient in the quadratic term of the Higgs potential), is proportional to the light neutrino masses and to M 3 , with M the heavy Majorana neutrino mass scale [15]. In order to avoid corrections to δµ 2 larger than say 1 TeV 2 one requires the lightest pair of Majorana neutrino masses to have M 1 < 4 × 10 7 GeV and M 2 < 7 × 10 7 GeV [16]. Secondly, there is a tendency for baryogenesis models to reside at the GeVor GUT-scales which leaves intermediate scales relatively unexplored. Finally, thermal leptogenesis at intermediate scales may resolve a problem that arises in the context of supersymmetric models which include gravitinos in their particle spectrum. Gravitinos have interaction strengths that are suppressed by the Planck scale and consequently are long-lived and persist into the nucleosynthesis era. The decay products of the gravitinos can destroy 4 He and D nuclei [17,18] and ruin the successful predictions of nucleosynthesis. Thus, in order reduce the number of gravitinos present at this stage, one requires a reheating temperature less than a few times 10 9 GeV depending on the gravitino mass [19].
The scale of leptogenesis may be lowered through the introduction of a symmetry to the SM. In [20], nonresonant thermal leptogenesis is explored at intermediate scales in the context of small B − L violation. It is shown that the DI bound may be evaded because, in the context of this near-symmetry, the lepton number conserving part of the CP asymmetries can be enhanced as they are not connected to light neutrino masses. It is shown that the scale may be lowered to 10 6 GeV. An alternative symmetry-based approach is to introduce supersymmetry in which one may also reduce the scale of leptogenesis to intermediate scales. In this context, the bound on the absolute value of the CP-asymmetry found by Davidson and Ibarra is greatly enhanced. Consequently, the DI bound is lowered thus allowing for the possibility of intermediate scale leptogenesis [21].
Beyond the application of supersymmetry and heavy pseudo-Dirac neutrinos, there are other means of lowering the scale of leptogenesis; if the decaying heavy Majorana neutrinos are near-degenerate in mass, the indirect CP-violation may be resonantly enhanced [6][7][8][9] and subsequently the mechanism may be lowered to the TeV scale. This has been explored in the context of type-I [22][23][24][25], II [26][27][28][29] and III [30,31] seesaw mechanisms. Another mechanism, proposed by [32], is one in which leptogenesis is achieved via CP-violating heavy Majorana neutrino oscillations [33][34][35][36]. The generation of the lepton asymmetry takes place close to the electroweak scale and the associated GeV-scale heavy Majorana neutrinos may be searched for at a variety of experiments such as LHCb [37,38], BELLE II [39] and the proposed facility, SHiP [40][41][42]. Although, leptogenesis via oscillations is a testable and plausible explanation of the baryon asymmetry, it has been shown its simplest formulations require a certain amount (∼ 10 5 ) of fine-tuning [43].
In this article, we revisit the question: how low can the scale of thermal leptogenesis go? We focus solely on the possibility that the heavy neutrinos are Majorana in nature and find thermal leptogenesis is possible at intermediate scales without resonant effects. In addition, we present an in-depth numerical study of the dependence of the baryon asymmetry produced from non-supersymmetric thermal leptogenesis on the low and high-scale model parameters.
The work presented in this paper is structured as follows: in Section II we review the origins of light neutrino masses in the type-I seesaw framework, further we review the Casas-Ibarra parametrisation of the Yukawa matrix and then introduce a modification of this parametrisation in the presence of large radiative corrections. We end this section by introducing a measure of fine-tuning in the context of the neutrino masses. In Section III we discuss the motivations for and some theoretical aspects of thermal leptogenesis. We follow in Section III A with the density matrix equations we shall solve to calculate the lepton asymmetry. We demonstrate in Section III B, that the fully flavoured Boltzmann equations, which do not incorporate flavour oscillations, may significantly qualitatively differ from the lepton asymmetry calculated from the density matrix equations and justify the use of semi-classical density matrix equations rather than kinetic equations derived from first principles non-equilibrium quantum field theory (NE-QFT). Our numerical methods are described in Section IV. The results of our numerical study for one and two decaying heavy Majorana neutrinos are presented in Section V A and Section V B respectively. In Section VI we explore the analytical consequences of the numerical results and provide an explanation for the fine-tuning. Finally, we summarise and make some concluding remarks in Section VII.

II. NEUTRINO MASSES AND THE TYPE-I SEESAW MECHANISM
In addition to the excess of matter over antimatter, the SM cannot account for non-zero neutrino masses which were discovered through the observation of neutrino oscillations by Super-Kamiokande twenty years ago [45] and were subsequently confirmed by a large number of other oscillation experiments. These oscillations occur because the neutrino flavour and mass eigenstates do not coincide. Such a misalignment between bases may be described by the Pontecorvo-Maki Nakagawa-Sakata, PMNS, matrix U , which relates the flavour and mass eigenstates of the One-loop level diagrams showing the physical particle contributions to neutrino mass at one-loop. The W -boson contribution proportional to the momentum is neglected. When, the external momentum is taken to be zero, the Z-and Higgs-boson contributions (Z 0 , H 0 ≡ φ 0 − v) together provide a correction to the tree-level mass that is finite and independent of the choice of renormalisation scale.  light neutrinos through and is conventionally parametrised as 2 where c ij ≡ cos θ ij , s ij ≡ sin θ ij , δ is the Dirac phase and α 21 , α 31 are the Majorana phases [46] which are physical if and only if neutrinos are Majorana in nature.
The current best-fit and 1σ range of the neutrino parameters [44] are provided in Table I. The Dirac CP-violating phase, δ, enters the neutrino oscillation probabilities sub-dominantly and remains mostly unconstrained by experimental data. However, as has been anticipated [47], the complementarity of long-baseline experiments such as T2K [48] and NOνA [49] with reactor 2 We have adopted the PDG parametrisation of the PMNS matrix [1].
In addition cosmology provides complementary bounds on the sum of the neutrino masses thanks to the imprint that neutrinos leave on the CMB and large-scale structure (LSS) in the early Universe. These bounds are significantly more stringent than the bounds from tritium-decay experiments with the 95% CL constraint m ν ≤ 0.2 eV [2]. In order to account for different analyses and underlying cosmological models we shall impose a constraint 3 m ν ≤ 1.0 eV, throughout this work. As discussed previously, the SM cannot explain neutrino masses in its minimal form. Arguably the simplest extension of the SM that incorporates small neutrino masses is the type-I seesaw mechanism [23][24][25]. This mechanism introduces a set of heavy Majorana neutrino fields N i and augments the SM La-3 For our best-fit points we find none that exceed mν = 0.63 eV (see Appendix D) and thus all are in within the more relaxed cosmological bound mν < 0.72 eV provided by Planck TT + lowP [2]. 4 grangian to include the following terms in which Y is the Yukawa coupling and Φ the Higgs SU (2) doublet, Φ T = φ + , φ 0 andΦ = iσ 2 Φ * , and L T = ν T L , l T L is the leptonic SU (2) doublet. For convenience we have chosen, without loss of generality, the basis in which the Majorana mass term is diagonal. In our work, we shall assume that there are three heavy Majorana neutrinos N i , with a mildly hierarchical mass spectrum in which M 1 < M 2 < M 3 .
After electroweak symmetry breaking, at the tree-level, the light neutrino mass matrix (at first order in the seesaw expansion) is 4 in which m D = vY is the Dirac mass matrix that develops once the Higgs acquires the vacuum expectation value v.
The tree-level mass matrix is not generically an accurate approximation of the light neutrino mass matrix over all regions of the parameter space. This is because there is no guarantee that the radiative corrections to the neutrino self-energy are negligible. Indeed there exist regions of parameter space in which radiative corrections are comparable to, or larger than, the tree-level contribution to the mass (see Table VI). For this reason, we find it necessary to incorporate the effects of the oneloop contribution to the masses given by [76] This parametrisation expresses the Yukawas in terms of both low energy measurable parameters (in m ν and U ) and high energy, currently untestable parameters (in the complex orthogonal matrix R and the Majorana mass matrix M ). An advantage to using this parametrisation is that one can automatically achieve the correct structure of mass-squared differences.
Naturally, when the loop-corrections are negligible we may replace f (M ) −1 with M and Eq. (3) reduces to the usual CI parametrisation. Throughout the remainder of this work, we shall apply the parametrisation of Eq. (3) to ensure radiative corrections are accounted for. We parametrise the R-matrix in the following way where c ωi ≡ cos ω i , s ωi ≡ sin ω i and the complex angles are given by In general the structure of the R-matrix cannot be constrained; however in [79], the authors demonstrated if the heavy Majorana neutrino mass matrix is invariant under a residual CP-symmetry, the R-matrix is constrained to be real or purely imaginary [79]. It was shown in [80,81] that the PMNS phases were a sufficient source of CP-violation to generate the observed baryon asymmetry if thermal leptogenesis occurred during an era in which flavour effects were non-negligible.
Both the light and heavy Majorana neutrino mass matrices determine the structure of the Yukawa matrix. Once the value of the lightest neutrino mass is fixed, 5 we shall assume the best-fit value for the solar and atmospheric mass squared splittings and hence the light neutrino mass matrix is determined. In addition, we constrain the sum of the neutrino masses such that it is within experimental bounds, specifically m ν ≤ 1.0 eV. In order to ensure the lepton asymmetry does not become resonantly enhanced [82], we choose the righthanded neutrino mass spectrum to be mildly hierarchical: [83]. In summary, the model parameter space of leptogenesis, as given by the Casas-Ibarra parametrisation of Eq. (3), is 18-dimensional and we shall henceforth denote this quantity as p.
In anticipation of our results, we shall define a parameter that quantifies the degree of fine-tuning for a given solution. First we define the fine-tuning measure to be where SVD[m 1-loop ] i and SVD[m ν ] i denote the ith singular values of the m 1-loop and m ν neutrino mass matrices respectively. As the neutrino mass matrix is the sum of the tree-and one-loop contributions, a cancellation between the two leads to the fine-tuning measure being larger than unity. In the limit that the tree-level contribution dominates, the fine-tuning measure tends to zero. We declare a technical limitation that we shall accept in this work. In lowering the value of M 1 , we find fine-tuned solutions in which the tree-level and one-loop contributions cancel to produce a neutrino mass matrix smaller than either alone. However, the higher-order radiative corrections cannot be assumed to perform a similar cancellation and thus we should take care that the two-loop contribution is not too large in comparison with the one-loop correct light neutrino mass matrix (see Appendix D).

III. THERMAL LEPTOGENESIS
Minimal thermal leptogenesis [4] proceeds via the outof-equilibrium decays of the heavy Majorana neutrinos. CP violation, arising from the interference between treeand loop-level diagrams, causes the CP-asymmetric decays of the heavy Majorana neutrinos which induce a lepton asymmetry. The production of the asymmetry from decays competes with a washout from inverse decays of the heavy Majorana neutrinos. The final lepton asymmetry is partially reprocessed to a baryon asymmetry via electroweak sphaleron processes which occur at unsuppressed rates at temperatures above the electroweak scale [5].
The time evolution of the lepton asymmetry may be calculated using semi-classical or NE-QFT methods. In both approaches, these kinetic equations account for the decay of the heavy Majorana neutrino and washout processes. In the simplest formulation, these kinetic equations are in the one-flavoured regime, in which only a single flavour of charged lepton is accounted for. This regime is only realised at sufficiently high temperatures (T 10 12 GeV) when the rates of processes mediated by the charged lepton Yukawa couplings are out of thermal equilibrium and therefore there is a single charged lepton flavour state which is a coherent superposition of the three flavour eigenstates. However, if leptogenesis occurs at lower temperatures (10 9 T 10 12 GeV), scattering induced by the tau Yukawa couplings can cause the single charged lepton flavour to decohere and the dynamics of leptogenesis must be described in terms of two flavour eigenstates. In such a regime, a density matrix formalism [84][85][86][87][88] allows for a more general description than semi-classical Boltzmann equations, since it is possible to calculate the asymmetry in intermediate regimes where the one and two-flavoured treatments are inadequate.

A. Density Matrix Equations
As previously mentioned, the most basic leptogenesis calculations were performed in the single lepton flavour regime. The one-flavoured regime is realised at high temperatures (T 10 12 GeV) where the leptons and antileptons that couple to the right-handed neutrinos, N i , maintain their coherence throughout the era of lepton asymmetry production. This implies there is a single lepton (anti-lepton) flavour, 1 ( 1 ), which may be described as a coherent superposition of charged lepton flavourstates, (e, µ, τ ), where the amplitudes c iα are functions of the Yukawa matrix, Y . In such a regime, the interaction rate mediated by the SM lepton Yukawas are out of thermal equilibrium (Γ α < H) and this implies there are no means of distinguishing between the three leptonic flavours. However, if leptogenesis occurs at lower scales (10 9 T (GeV) 10 12 ), the interactions mediated by the tau charged lepton Yukawa come into thermal equilibrium (Γ τ > H) and the Universe may distinguish between τ and τ , where τ is a linear combination of the electron and muon flavoured leptons orthogonal to τ . The one-and fully two-flavoured description of leptogenesis is appropriate at T 10 12 GeV and 10 9 T 10 12 GeV, respectively. There exists the possibility that thermal leptogenesis occurs at even lower temperatures, T < 10 9 GeV, during which the interactions mediated by the muon have equilibrated. In such a regime the kinetic equations should be given in terms of all three lepton flavours. In this work, we shall focus on this particular scenario. Our discussion of the density matrix equations and notation will closely follow the prescription of [88] where the theoretical background and derivation of the density matrix equations are 6 fully presented. We refrain from rederiving the details of this formalism and instead refer the interested reader to the aforementioned reference. The most general form of the density matrix equations, assuming three decaying heavy Majorana neutrinos, is given by where Greek letters denote flavour indices, n Ni (i = 1, 2, 3) is the abundance of the ith heavy Majorana neutrino 5 , n eq Ni the equilibrium distribution of the ith heavy Majorana neutrino, D i (W i ) denotes the decay (washout) corresponding to the ith heavy Majorana neutrino, and are given by [14] and with K 1 and K 2 the modified Bessel functions of the second kind with H is the Hubble expansion rate and Λ α is the self-energy of α-flavoured leptons. The thermal widths Λ τ , Λ µ of the charged leptons is given by the imaginary part of the self-energy correction to the lepton propagator in the plasma (see Appendix A). Finally, the P 0(i) αβ ≡ c iα c * iβ , denotes the projection matrices which describe how a given flavour of lepton is washed out and the CP-asymmetry matrix describing the decay asymmetry generated by N i is denoted by αβ . These CP-asymmetry parameters may be written as [6,86,88,89] (i) and Eq. (6) may be used to calculate the lepton asymmetry in all flavour regimes and even accurately describes the transitions between them [84][85][86][87][88]. The off-diagonal entries of n αβ , which in general may be complex, allow for a quantitative description of these transitions. If leptogenesis proceeds at temperatures 10 9 T (GeV) 10 12 , (for example), the terms (Λ τ ) /Hz damp the evolution of the off-diagonal elements of n αβ . This reflects the loss of coherence of the tau states when the SM tau Yukawa couplings come in to equilibrium. The remaining equations provide a description of leptogenesis in terms of Boltzmann equations for the diagonal entries of n αβ and for n Ni . Although a more accurate treatment of leptogenesis is provided by the NE-QFT approach, the density matrix equations that we choose to solve are accurate so long as we are in the strong washout regime. In Appendix B, we demonstrate that this is the case and thus justify our use of the density matrix formalism. In general, n αβ is a 3 × 3 matrix whose trace gives the total lepton asymmetry: The latter is then multiplied by a factor a/f ≈ 0.01, where a = 28/79 describes the partial conversion of the B − L asymmetry into a baryon asymmetry by sphaleron processes, and f ≡ n rec γ /n * γ = 2387/86 accounts for the dilution of the asymmetry due the change of photon densities (n γ ) between leptogenesis (n γ = n * γ ) and recombination (n γ = n rec γ ) : η B 10 −2 n B−L [14].

B. Thermal Leptogenesis with Three Flavours
In this section we demonstrate, in the case of one decaying heavy Majorana neutrino, the need to solve the density matrix equations, rather than the more approximate Boltzmann equations (in which the off-diagonal entries of the density matrix are set to zero). Although we shall give explicit expressions for the density matrix equations in this section, we emphasise that our main results are always found by solving Eq. (6) in which all three heavy Majorana neutrinos decay. The Boltzmann equations, with one decaying heavy Majorana neutrino, are written as where we have used the abbreviation 11 . This set of equations is appropriate for M 1 10 9 GeV, when the flavour-components of the charged leptons each experience strong and distinct interactions with the early Universe plasma. The density matrix equations, with specified flavour indices, follow straightforwardly from an explicit expansion of the commutators in Eq. (6). The resulting equations, for a single decaying heavy Majorana neutrino are The Boltzmann equations, Eq. (12), are recovered in the limit (Λ µ ) /Hz, (Λ τ ) /Hz → ∞ as the off-diagonal density matrix elements become fully damped. This limit is only valid for M 1 10 9 GeV. However, for a given point in the model parameter space, p, it is not a priori obvious if these Boltzmann equations well approximate the density matrix equations.
We illustrate the quantitative difference between the density matrix equations (Eq. (13)) and the Boltzmann equations (Eq. (12)) by solving both for four benchmark points: BP1, BP2, BP3 and BP4 with a vanishing initial abundance of N 1 , (see Table II). In these scenarios, the light mass spectrum is chosen to be normally ordered, m 1 = 10 −2 eV, M 1 is allowed to vary with M 2 = 3.5 M 1 M 3 = 3.5 M 2 which satisfies a mildly hierarchical mass spectrum.
As can be seen from Fig. 2, the deviation between the two is generally small (< 5%) for M 1 ∼ 10 6 GeV. In the case of BP1, the Boltzmann equations do not deviate from the density matrix solutions until M 1 ≈ 10 9 GeV and for M 1 ≈ 10 6 GeV, the discrepancy between the two is negligible. A more pronounced deviation is exhibited in BP2, BP3 and BP4, with an underestimation from the Boltzmann equations particularly evident in BP4 in which the R-matrix has relatively large elements. In this example, the deviation between the solutions even for low masses M 1 = 10 7 GeV and M 1 = 10 6 GeV is ∼ 20% and ∼ 5% respectively. The discrepancy grows as a function of M 1 . As can be seen from these benchmark points, the fully three-flavoured equations may not well approximate the density matrix equations well even for M 1 10 9 GeV. As we are interested in exploring the parameter space over a range of values of M 1 , we shall use the more accurate density matrix equations.
Here we summarise the approximations and physical effects that we shall exclude from our calculation but whose inclusion would increase the accuracy of our calculations. Such effects include lepton number-changing scattering processes, spectator effects [90][91][92], thermal corrections [93,94] and the inclusion of quantum statistical factors [95][96][97][98]. |∆L| = 1 scattering and related washout processes occur as a result of Higgs and lepton mediated scattering involving the top quark and gauge bosons. It has been demonstrated that scatterings involving the top quark are most important at relatively low temperatures, T < M [99]. Therefore, the effects of |∆L| = 1 scattering on the strong washout regime (where the bulk of lepton asymmetry is produced at T > M ) are small and have been estimated to affect the final lepton asymmetry to a level less than ∼ O(10)% [100,101]. However, for weak washout these corrections are necessary for a correct calculation of the final lepton asymmetry [102]. Spectator processes cause the redistribution of the asymmetry generated in the leptonic doublets amongst other particle species in the thermal bath. These processes typically protect the lepton asymmetry from washout and therefore increase the efficiency of leptogenesis [103]. Although the inclusion of spectator effects could further lower the scale of successful thermal leptogenesis, we relegate the inclusion of these effects for further studies. Besides, the neglect of these effects leads to an overly-conservative estimate of the amount of baryon asymmetry produced. Quantum kinetic equations are derived from the first principles of NE-QFT based on the Closed-Time Path (CTP) formalism. This approach resolves unitarity issues and properly accounts for the effect of quantum statistics on the lepton asymmetry. However, it has been shown there is little qualitative difference between the density matrix and CTP approach in the strong washout regime [98]. This is because in the strong washout regime, where the decays and inverse decays of the heavy Majorana neutrinos occur much faster than the expansion rate of the Universe, the majority of the lepton asymmetry is produced at temperatures greater than the mass of the decaying heavy Majorana neutrinos. As a consequence, the contributions of the particle distribution functions are heavily Boltzmann-suppressed.
We demonstrate in Appendix B that the regions of parameter space we explore in this work correspond to strong washout. This allows us to make two justifiable simplifications to our kinetic equations which are more easily implementable for a phenomenological study. Firstly, we ignore the impact of lepton number-changing scatterings and secondly we solve kinetic equations using the density matrix formalism rather than equations derived from NE-QFT.

IV. COMPUTATIONAL METHODS
The computational core of this work is solving a set of coupled differential equations as shown in Eq. (6). We use the Python interface for complex differential equations [104] to the LSODA algorithm [105] that is available in Scientific Python [106].
Our aim is to find regions of the model parameter space, p, that yield values of η B (p) that are consistent with the measurement η BCMB = (6.10 ± 0.04) × 10 −10 . In order to do so, we have to use an efficient sampling method. This is mainly for three reasons. Firstly, the parameter space has a relatively high dimension. Secondly, the function η B (p) itself does not vary smoothly with changes of p. In fact, tiny variations of the input parameters yield function values differing in many orders of magnitude and sign. Thirdly, the computation of η B (p) for a single point is relatively expensive and can take up to the order of seconds. Thus any attempt of a brute-force parameter scan is doomed to fail. Finally, we are not only interested in a single best-fit point but also a region of confidence that resembles the measurement uncertainty.
We found the use of Multinest [107][108][109] (more precisely, pyMultiNest [110], a wrapper around Multinest written in Python) to be particularly well suited to address all the aforementioned complications associated to this task. The Multinest algorithm has seen wide and very successful application in astronomy and cosmology. It provides a nested sampling algorithm that calculates Bayesian posterior distributions which we will utilise in order to define regions of confidence.
In all our scenarios, Multinest uses a flat prior and the following log-likelihood as objective function Once a Multinest run is finished, we use Super-Plot [111] to visualise the posterior projected onto a two-dimensional plane.

V. RESULTS
We present the solutions to the density matrix equations of Eq. (6) for the case of one and two decaying heavy Majorana neutrinos in Section V A and Section V B, respectively. In principle, it is necessary to consider the decay of all three heavy Majorana neutrinos, however we first consider the decay of the lightest heavy Majorana neutrino as computationally this scenario is less expensive than the two and three decaying heavy Majorana neutrinos case. In Section V B, we demonstrate the scale of thermal leptogenesis involving the decay of two heavy Majorana neutrinos does not change significantly from the scenario of one decaying case. These two scenarios are qualitatively and quantitatively similar and so we do not proceed to the case where the third heavy Majorana neutrino contributes to the baryon asymmetry through their decays.

A. Results from N1 Decays
As detailed in Section IV, solving the density matrix equations of Eq. (13) for an 18-dimensional model parameter space, p, is a challenging numerical task. In order to reduce the volume, p, we shall fix certain parameters. Firstly, as the solar and reactor mixing angles are relatively precisely measured, we shall use the values for these angles from global fit data [44]. Although we best-fit points of S1, S2 and S3 (S1, S2 and S3) respectively.  allow the lightest neutrino mass (m 1 for NO and m 3 for IO) to vary within the experimentally allowable region, given by the sum of neutrino masses, the other two light masses are determined from the best-fit values of the atmospheric and solar mass squared splittings from global fit data [44]. Finally, we fix the heavy Majorana mass spectrum leaving only 11 of the 18 parameters of p to be varied.
In all scenarios we choose a set of initial values for M 1 , M 2 and M 3 , in which, as mentioned, we ensure M 3 > 3M 2 and M 2 > 3M 1 . We explore the parameter space and find the regions consistent with η B CM B to a 1 and 2σ level. Through the inspection of the fine-tuning of the solutions in the regions of 1σ agreement, we decide either to lower the scale of M 1 or not (while keeping the ratios M 2 /M 1 and M 3 /M 2 fixed). The lower the scale, the higher the fine-tuning and thus the greater the impact of higher-order corrections. Thus we do not further lower the scale when either the two-loop contributions becomes greater than a few percent or when the finetuning exceeds O(1000) (see Appendix D) . If one were to incorporate the effects of higher radiative orders, the parameter space could be explored at even lower scales where the fine-tuning is greater.
In the case of one decaying heavy Majorana neutrino contributing to the lepton asymmetry, we shall pick out six scenarios in total, all of which satisfy a mild hierarchical spectrum:  for normally ordered light neutrino masses. In the case of inverted ordering, we shall denote these scenarios as S 1 , S 2 and S 3 . Scenarios S 1 (S 1 ) and S 2 (S 2 ) have the same mass ratios, with S 1 (S 1 ) corresponding to the lowest value of the scale M 1 with acceptable fine-tuning values and S 2 (S 2 ) presented for comparison. S 3 (S 3 ) corresponds to the lowest scale for its given set of mass ratios. In Fig. 3, we provide the temperature evolution of the absolute magnitude of the lepton asymmetry number densities, |n αα |, α = e, µ, τ for the best-fit points of each scenario.
The parameters of the PMNS matrix are varied within their allowable or measured 3σ range: δ ∈ (0, 360) • ,  For the two-dimensional posterior plots of scenario S 1 , as shown in Fig. 4, the region of the model parameter space consistent to a 1σ level with the observed baryon asymmetry favours larger values of the CP-violating Dirac phase, 120 ≤ δ( • ) ≤ 360. The likelihood function appears to be more sensitive to α 21 than α 31 : from Fig. 4, we observe 80 ≤ α 21 ( • ) ≤ 270 while 65 ≤ α 31 ( • ) ≤ 720 is consistent with the measured baryon asymmetry to a 1σ level. Although the atmospheric mixing angle may take most value within its 3σ range, the likelihood function favours values close to 45 • and in the upper octant. The values of the lightest neutrino mass which are consistent with the observed η B tend to be close to the upper limit, which for normal ordering is m 1 3.32 × 10 −1 eV. This strong dependence of η B on the lightest neutrino mass agrees with work which investigated (two) flavoured thermal leptogenesis [112].
In general, the likelihood function is more sensitive to the imaginary than the real components of the R-matrix. For example, we find that η B is relatively insensitive to x 1 and x 3 : x 1 , x 3 ∈ (−180 • , 180 • ) is consistent with the measured η B to 2σ level. On the contrary, the likelihood function is highly sensitive to x 2 with preferred values of approximately 90 • . We note that the two-dimensional projections onto parameters x 1 and x 3 are not included in the triangle plots as η B (p) exhibits flat directions in both these parameters and the two-dimensional projection plots show little interesting structure 7 . The complex components of the R-matrix are likely to be within a small range: y 1 180 • , y 2 3 • and y 3 180 • where 7 However, these plots are included in the aforementioned link.
the explanation for this structure has been detailed in Section VI. Given the mass of the decaying heavy Majorana neutrino is relatively light, it would be expected that large phases of the PMNS and R-matrix are favoured as these ensure the Yukawa couplings are sufficiently large.
The triangle plots for larger masses of M 1 and more hierarchical heavy Majorana neutrino spectra of S 2 and S 3 are shown in Fig. 10 and Fig. 11 respectively. Unsurprisingly, on comparison of scenario S 1 and S 2 (which share the same mass splitting but different values of M 1 ) we observe the scenario with the larger heavy Majorana neutrino masses has a larger region of the model parameter space consistent with the measured η B CM B . Moreover, as expected, the constraints on the R-and PMNS-matrix parameters in scenario S 2 are weaker yet qualitatively similar to S 1 . In particular, the m 1 -dependence in S 2 is less severe than in the scenario of S 1 ; for example in Fig. 4 the 2σ allowed region for the lightest neutrino mass is 1.25×10 −1 ≤ m 1 (eV) ≤ 3.32×10 −1 while in the case of Fig. 10 For smaller values of m 1 , successful leptogenesis is possible for larger values of the heavy Majorana neutrino mass M 1 . For larger heavy Majorana neutrino mass splitting, we anticipate the model parameter volume consistent with data will be reduced. This is because the CP-asymmetry becomes increasingly suppressed for larger mass splittings. This effect is confirmed upon comparison of Fig. 10 and Fig. 11 where the former has milder mass splitting. In contrast to S 1 , in the case of both S 2 and S 3 , the likelihood function favours values of θ 23 close to 45 • and in the lower octant.
The triangle plot showing the two-dimensional posterior distributions of the 11-dimensional model parameter space for S 1 is shown in Fig. 5. The dark (light) red contours correspond to the regions of parameter space consistent with 68% (95%) confidence levels. As anticipated, the points of the model space consistent with the measurement are different from the normal ordering case and the volume of parameter space p consistent with the measured η B CM B is less constrained. In particular we observe that the likelihood function is relatively insensitive to changes of δ, α 31 and θ 23 . However, this scenario displays a similar feature to S 1 , where the likelihood function favours values of α 21 ≤ 360 • .
Additionally, the likelihood has a flat direction in the x 1 and x 3 parameters of the R-matrix (as discussed in Section VI). We observe that all values of x 1 and x 3 are consistent to a 2σ level with the measured η B ; however, the likelihood is very sensitive to x 2 with x 2 90 • . Similarly, to the normal ordering scenario the imaginary phases of R are constrained with y 1 180 • , y 2 2 • and y 3 180 • . The triangle plots for larger masses of M 1 and more hierarchical spectra of S 2 and S 3 are shown in Fig. 12 and Fig. 13 respectively. As seen in the case of normal ordering, the scenario with the slightly more hierarchical mass spectrum (M 2 = 5 M 1 , M 3 = 5 M 2 ) has a slightly smaller volume of parameter space consistent with the data than the case of the milder hierarchy.
As mentioned previously, we remain agnostic about the level of fine-tuning Nature finds tolerable and simply present the fine-tuning measure defined in Eq. (5) for the regions of the model parameter space within 1σ of the measured η B . To be explicit, the top (bottom) three plots of Fig. 6 shows the distribution of the fine-tuning measure within the 1σ region of S 1 , S 2 and S 3 (S 1 , S 2 and S 3 ) shown in Fig. 4, Fig. 10 and Fig. 11 (Fig. 5,  Fig. 12 and Fig. 13) respectively. It should be noted we did not veto points in p which exceeded a given value of fine-tuning, moreover increasing the spread from 1σ to 5σ would allow for a broader spread of fine-tuning values, both smaller and larger.
In general, for normal ordering, the fine-tuning measure for points within 1σ is O (100). We observe that the minimal fine-tuning value for S 1 ≈ 330. Somewhat unsurprisingly, the scenario with the large mass of decaying heavy Majorana neutrino, S 2 , has smaller fine-tuning due to the fact the complex phases of the R-matrix may attain a broader range of values. We observe the minimum fine-tuning measure in the case of S 2 to be ≈ 180. However, in the case of S 3 (where the decaying heavy Majorana neutrino mass is the same as S 2 the mass splitting between the heavy Majorana neutrinos is larger) the finetuning values are in general larger due to the increased mass of N 3 .
The fine-tuning present in the case of inverted ordering is, in general, less than in the case of normal ordering. The minimum value of fine-tuning present in S 1 100. Again, the same pattern emerges as in the case of normal ordering where the fine-tuning in S 2 (S 3 ) is less (greater) than S 1 . In fact, for S 2 the minimum fine-tuning ≈ 40. Again, we emphasise the fine-tuning we present here is for points in p within 1σ of the best fit value of η B CM B and allowing for an increase in the spread around the best fit value would allow for smaller (and larger) values of fine-tuning.
At such scales, T 10 9 GeV, it is impossible to have successful leptogenesis without some degree of cancellation between the tree and one-loop level contributions. However, we did investigate if there existed regions of p such that thermal leptogenesis was viable (within 1σ of the central value of η B CM B ) where either the tree or one-loop level contribution dominates. In the latter scenario, where the radiative corrections dominate over the tree-level contributions, the fine-tuning measure should be close to unity as |m 1-loop |/| m tree + m 1-loop | ≈ 1 for m tree m 1-loop . We applied the same numerical procedure to solve the density matrix equations with one decaying heavy Majorana neutrino and vetoed points in p if the fine-tuning measure was not within the boundary 0.9 ≤ F.T ≤ 1.1. After scanning a series of differing heavy Majorana neutrino mass spectra, we found the loop-dominated scenario was possible, assuming normal ordering, for M 1 = 10 9 GeV with M 2 = 3.15M 1 and M 3 = 3.15M 2 . The best-fit point is denoted as F.T loop in Table III and the triangle plot of the twodimensional posterior distributions may be found on the provided webpage. In the former scenario, where the tree-level contributions dominates, the fine-tuning measure will be close to zero. Using Multinest to search for regions of p consist with tree-domination we required the fine-tuning to be within the boundary 0 ≤ F.T ≤ 0.2. We found no solutions compatible with this condition for M 1 < 10 9 GeV. However, we did find a single single point consistent with a fine-tuning ≈ 0.18 for a mass spectrum of M 1 = 10 9 GeV, M 2 ≈ 3.15M 1 and M 3 ≈ 3.15M 2 . Note that a two-dimensional projection of the posterior is not possible and we simply provide the value of this point as F.T tree in Table III. For larger values of M 1 more points will exist that satisfy the condition and so we regard F.T tree as the solution of lowest M 1 in which the tree-level is the dominant contribution. The absolute values of the Yukawa matrix elements are listed for all scenarios for reference in Appendix E. We note that it is possible to reduce the fine-tuning by considering the scenario where M 2 = M 3 . Such a scenario may result from the introduction of a partial symmetry into the type-I seesaw. As, in this section, we only consider the case that N 1 decays this does not lead to resonant leptogenesis. As an example, consider S 1 but with M 2 = M 3 ≈ 5.05 × 10 6 GeV. Such a point in p leads to η B = 6.1 × 10 −10 , which is in good agreement with the experimental value. In this case, N 2 and N 3 act as two Majorana components of a pseudo-Dirac pair. The contribution of N 2 and N 3 to the tree-level mass is cancelled (as together they are lepton number conserving) and a dramatic reduction in our fine-tuning measure occurs, resulting in F.T. ≈ 2.1. This is similar to the scenarios considered in [20] and will not be further discussed in this paper.
In summary, taking advantage of flavour effects, a mildly hierarchical heavy Majorana neutrino mass spectrum and fully exploring the Casas-Ibarra model parameter space p allows for successful thermal leptogenesis at surprisingly low scales with M 1 = 10 6 GeV for both normal and inverted ordering. In the case of normally ordered light neutrinos, larger values of the δ are favoured in conjunction with an atmospheric mixing angle close to θ 23 = 45 • (slightly above or below depending on the scenario, see Table III). We observe that larger masses of m 1 are favoured as this compensates for decreasing M 1 . In the scenario of an inverted ordered mass spectrum, the likelihood function shows little sensitivity to changes in the low-energy neutrino parameters. On the other hand, the R-matrix is comparatively highly constrained. In addition, we present the distribution of the fine-tuning measure within 1σ of the measured η B and found the fine-tuning was in general smaller for inverted than normal ordering and usually took values ∼ O (100). The minimal fine-tuning of the six scenarios we studied was found for S 2 which achieved F.T. ∼ 40. Moreover, we anticipate allowing regions of p consistent with 5σ of the central value of η B would allow for still smaller values of fine-tuning.

B. Results from N2 Decays
In this section, we explore the possibility that the decay of two heavy Majorana neutrinos contributes to the baryon asymmetry. In this setup, the density matrix equations follow rather straightforwardly from Eq. (6) and the numerical procedure to find the two-dimensional posterior plots is the same as discussed in Section V A. The qualitative difference between this case and the former as discussed in Section IV is that now N 2 may decay in addition to N 1 . As M 2 > M 1 , N 2 will decay before N 1 with the average time between the two decays determined by the hierarchy of their masses.
In [113] the authors explored thermal leptogenesis using the decay of two heavy Majorana neutrinos in the limit the third is decoupled from the theory. Using analytic estimates, they found the minimal mass of the lightest heavy Majorana neutrino, for successful leptogenesis, to be M 1 ∼ 1.3 × 10 11 GeV assuming a mildly hierarchical mass spectrum. In this scenario, we explored a number of heavy Majorana neutrino mass scenarios and found the lowest mass of N 1 which allowed for successful leptogenesis was M 1 = 10 6.7 GeV with M 2 ≈ 6.3M 1 and M 3 ≈ 4M 2 . We denote these two scenarios as S 4 and S 4 for normal and inverted ordering respectively and the best-fit point and corresponding triangle plots are shown in Fig. 8 and Fig. 9.
Naively, one would think that the decay of two heavy Majorana neutrinos would further lower the scale of leptogenesis; however, this is not the case as there is nontrivial interplay between the decays and washout processes of N 2 and N 1 . We note that contribution of the third heavy Majorana neutrino to the lepton asymmetry in these scenarios is negligible as the CP-asymmetry (3) αβ is several orders of magnitude lower than that of the other two and its washout term W 3 decays far faster.
Unlike in the previous section, we find the twodimensional posterior projections in this case for both orderings do not appear to be too dissimilar. In both cases, the likelihood function is insensitive to δ. In addition, the atmospheric mixing angle can be in the lower or upper octant and there is strong dependence on large values of m 1 (m 3 ) in S 4 (S 4 ). The dependence of the likelihood on the R-matrix parameters is similar to the cases discussed in Section V A; we find x 1 and x 3 may take any values while x 2 90 • . Likewise, two of the imaginary components of the R-matrix are constrained to be large y 1 , y 3 180 • while the other is nearly vanishing y 2 2.5 • . For reference, the corresponding absolute value Yukawa matrices are given in Section E. In a similar fashion to Section V A, we present the fine-tuning measure for the regions of the model parameter space within 1σ of the measured η B . We observe for normal and inverted ordering the fine-tuning ∼ O (100).

VI. DISCUSSION OF FINE-TUNED RESULTS
We may gain an understanding of why fine-tuned solutions were found by the numerical machinery through inspection of the structure of the Yukawa matrix at the best-fit points. Looking at the solutions for one and two decaying heavy Majorana neutrino scenarios, we observe that generically |y 1 | ≈ 180 • , y 2 ≈ 0 • , |y 3 | ≈ 180 • and |x 2 | ≈ 90 • . Consider as a typical example S 1 , for which the orthogonal R-matrix assumes the following form which has the structure The appearance of y 1 and y 3 in the exponentials, and the proximity of In the case of the asymmetries (1) αα , generated in the N 1 decays, and for the best-fit values of the parameters listed in Table III,    Under the approximation m 1 = m 2 , the part of the asymmetry proportional to f 1 (which we call

and
(1) Numerical estimates at the best-fit values of Table III show that this second contribution (the resonance function contribution) is somewhat larger than the first one, although, the baryon asymmetry in the cases studied by us is produced in the non-resonance regime.
In the density matrix equations, the CP-asymmetry parameters enter in the combinations (1) in the three-flavour regime. Thus, although for our best-fit scenarios τ τ (f 2 ) may be O 10 −22 , this does not mean that the (1) αα (f 2 ) give a negligible contribution in the generation of the lepton (baryon) asymmetry.
We note that there is a factor (Y † Y ) −1 11 in the diagonal CP-asymmetries (1) αα (Eq. (10)) for the lightest heavy Majorana neutrino and a factor (Y † Y ) 11 (Eq. (8)) appears in the washout term W 1 . Thus, we naively expect that in order achieve successful leptogenesis, by reducing the washout, (Y † Y ) 11 should be made small. Expanding this quantity, in terms of the R-matrix elements and the remaining CI parameters, we find Thus, with the assumption that this quantity should be small, the relative smallness of the elements R 1i is explained and with it the values of x 2 and y 2 . Similarly, given the dependence on |R 21 | in (1) αα (f 2 ), it may be expected that we should maximise the values of y 1 and y 3 . With these imaginary parts of ω 1 and ω 3 large, the values of the corresponding real parts x 1 and x 3 is immaterial. This is reflected in the relative flatness of their directions in the parameter space plots. The dependence on m 1 in (Y † Y ) 11 may initially lead one to expect m 1 to be minimised. That this is not the case is due to the factors m 2 1 or m 3/2 1 √ m 3 appearing in the expressions for (1) αα . In order to maximise these CP-asymmetries, one would expect m 1 to be found at its largest allowed value (determined by the constraint on the sum of the neutrino masses).
Let us now examine how these choices of parameters affect the expressions for the tree-and one-loop light neutrino masses. We may estimate the light masses using  the largest value of the Yukawa matrix (∼ 10 −2 in the case of S 1 , see Appendix E) and the smallest heavy mass M 1 = 10 6 GeV: This mass is too large from the point of view of the experimental bound and yet the numerical machinery is enforcing neutrino masses which sum to < 1 eV. Let us investigate why this estimate fails. This structure of the R-matrix leads to the following structure for the Dirac  mass matrix: in which |δ 2 | |δ 1 | |u| where each of δ 1 , δ 2 and u are 3-component complex vectors. We may rewrite the tree and one-loop masses in terms of this relatively simple where the commutativity of the diagonal matrices M and f has been exploited. For the one-loop contribution we find

This ensures that their sum is
Due to the relative smallness of the elements of δ i , the light neutrino mass matrix may be considerably smaller than would be expected from a naive estimate based on the size of u. Neglecting terms containing a δ i , we find that This is the mechanism by which the fine-tuned mass matrices are arrived at. Although in this analysis, the results of S 1 were used, the other solutions differ essentially only in the sign used for y i . This introduces a different pattern of minus signs in the matrix of Eq. 15 (and hence also in the expression for m D √ f ) which does not affect the overall argument. Note that this argument is true even for the solutions of the two-decaying heavy Majorana neutrinos equations.

VII. SUMMARY AND CONCLUSIONS
In this work we have explored the viable model parameter space of thermal leptogenesis associated to a type-I seesaw mechanism. To do so, we numerically solved the three-flavoured density matrix equations [88] for one and two-decaying heavy Majorana neutrinos. Of the eighteen dimensional model parameter space, seven parameters were fixed from neutrino oscillation data, cosmological constraints and consideration of a mildly hierarchical heavy Majorana neutrino mass spectrum.
To find the regions of parameter space consistent with the measured baryon-to-photon ratio we used pyMulti-Nest which implements a nested sampling algorithm to calculate Bayesian posterior distributions which are utilised to find regions of confidence. In addition, we ensured the Yukawa matrix entries respected perturbativity and we protected against resonance effects by assuming a mildly hierarchical heavy Majorana neutrino mass spectrum. In the case of one decaying heavy Majorana neutrino, we found the lightest heavy Majorana neutrino mass that could successfully generate the baryon asymmetry, with our choice of upper bound on R-matrix components, to be M 1 10 6 GeV. Naively put, this is possible as the eleven parameters which were allowed to vary compensated for the smaller heavy Majorana neutrino masses. Moreover, with normal ordering, maximally CPviolating values of δ and θ 23 close to 45 • (in most cases slightly larger than 45 • , see Table III) is preferred. In addition, there was strong dependence on the mass of the lightest neutrino. On the other hand, we found in the case of inverted ordering there were no strong constraints on low energy neutrino parameters. For this scenario, the level of fine-tuning was ∼ O (100). In the case of one decaying heavy Majorana neutrino, we found the scenario with the smallest fine-tuning, at intermediate scales, was S 2 , (F.T ∼ 40) with a heavy Majorana neutrino spectrum M 1 = 10 6.5 GeV, M 2 ≈ 3.15M 1 and M 3 ≈ 3.15M 2 . We showed also that fine tuning would not be necessary at all if M 2 = M 3 , when the one loop contribution to the light Majorana neutrino mass matrix is strongly suppressed. We also explored the possibility that either the tree or one-loop radiative corrections dominate the neutrino mass matrix. We found the lowest scale possible for this scenario, assuming a mildly hierarchical spectrum, was M 1 = 10 9 GeV. As discussed, a motivation for exploring leptogenesis at intermediate scales is to avoid large corrections to the Higgs mass. Although, we found regions of the parameter space of three-flavoured thermal leptogenesis consistent with the observed baryon asymmetry, we did not seek to minimise δµ 2 and relegate this to a future study.
Finally, we investigated the case of two decaying heavy Majorana neutrinos. We found the lowest scale for both normal and inverted ordering to be M 1 = 10 6.7 GeV. This scale is higher than in the one decaying heavy Majorana neutrino case because the scale of the washout is larger for N 2 and its CP-asymmetry is small in comparison with N 1 . Although the washout for N 2 decays much more quickly than for N 1 , it still has an appreciable effect on the final lepton asymmetry and so one must raise the scale of the heavy Majorana neutrino masses to achieve successful leptogenesis. In this paper, we did not include spectator effects which could potentially further lower the scale of thermal leptogenesis and may be investigated in future work.

ACKNOWLEDGMENTS
We would like to thank Pasquale di Bari, Bjoern Garbrecht and Andrea de Simone for useful discussions and advice regarding this work. We are grateful to Alexis Plascencia, Carlos Tamarit For the normally ordered mass spectrum, following experimental constraints on the masses, m 2 and m 3 are increasing functions of m 1 . Thus, if the elements of R are fixed, K 1 is smallest when m 1 = 0. A random scan over the angles of R (allowing x i in [0, 360] • and y i in [0, 180] • ) for 10 6 points leads to the conclusion that > 99.9% of points in the parameter space lead to K 1 > 1 and ∼ 99.7% of points lead to K 1 > 10.
In IO, a random scan of 10 6 points found none for which K 1 < 1 only 9 points for which K 1 < 10. Thus, the experimental constraints in both the IO and NO case greatly favour strong washout.
In conclusion, very few points in the parameter space satisfy the experimental constraints on the neutrino mass-squared differences and mass-sum whilst simultaneously achieving weak washout. Thus it is a safe assumption that the washout is strong and our numerical methods are accurate.

Appendix C: The resonance region
The analytical expressions used CP asymmetry parameters have been calculated under the assumption that the the heavy Majorana neutrinos have well-separated masses such that the usual Feynman rules may be used in perturbation theory. The meaning of well-separated here is such that the mass differences are significantly larger than their decay rates. In this appendix we investigate this assumption.
The total CP asymmetry parameter is defined in terms of Γ 1 , the decay rate for N 1 → φ † l andΓ 1 , the rate for 21 CP conjugate process N 1 → φl † , as and the decay terms are As analytical expressions for these are well-known we may put them to use in finding the decay rate. We have The Hubble parameter H in a radiation-dominated Universe is, from the Friedmann equation, which may be expressed in terms of M 1 , z and the Planck mass M P with Thus, the decay rate may be written in terms of the functions that are typically written in the Boltzman equations by In order to avoid the resonance region we require that To test this, the PMNS angles; x 1 , x 2 , x 3 ; and M 1 , M 2 , M 3 were fixed according to the best-fit points for NO (Fig. 4, Fig. 10 and Fig. 11) and also for IO (Fig. 5, Fig. 12 and Fig. 13) and a random scan over the remaining parameters for 10 5 points was performed with the criterion We found there were no points which verified this condition and thus the assumption of non-resonance is justified.
Appendix D: Higher-order radiative corrections We have been careful to include the one-loop radiative corrections to the light neutrino masses. In doing so we  have expanded the region of the parameter space in which we may accurately explore leptogenesis. Of course, there may also be regions in which the higher-order corrections are important. We may ask the question how can we be sure that the neglect of two-loops, three-loops etc. was legitimate?
A pragmatic approach is to perform an order-ofmagnitude estimate of the effects of the higher-order corrections for those points in the parameter space of most significance to our result: the best-fit points for the scenarios S 1 to S 4 and F.T loop , F.T tree . If, in these scenarios, the higher-order corrections appear small then our main conclusions are left untouched.
Our estimate of the two-loop effect (which we shall assume generically dominates three or more loops) will be given by two extra factors of the Yukawa couplings and the conventional loop factor (4π) −2 to the one-loop effect. Let us use with |Y max. | the largest element of the matrix of absolute values of the Yukawas, as a conservative estimate (over-estimate) of the second-order radiative correction to neutrino masses. (This is similar to the estimate used in [78].) From table VI, we see that the two-loop contributions generally provide small corrections and therefore that corrections beyond one-loop order are safely neglected at these points.

22
Appendix E: Yukawa matrices Here we provide a table of the absolute values of the Yukawa matrices (|Y |) for the best-fit points of each scenario considered in Table III and Table IV.   Table III and Table IV