Flavoured Resonant Leptogenesis at Sub-TeV Scales

We consider sub-TeV scale flavoured resonant leptogenesis within the minimal type-I seesaw scenario with two right-handed singlet neutrinos $N_{1,2}$ forming a pseudo-Dirac pair, concentrating on the case of masses of the pseudo-Dirac pair having values $M_{1,2} \lesssim 100$ GeV. The case when the CP violating asymmetries in the individual lepton charges $L_l$, $l=e,\mu,\tau$, and in the total lepton charge $L$ of the Universe are generated in $1 \leftrightarrow 2$ decay processes is investigated. We show that successful leptogenesis is possible for $M_{1,2}$ lying in the interval $M_{1,2} = (0.3 - 100)$ GeV. Our results show also, in particular, that for vanishing initial $N_{1,2}$ abundance, flavour effects can play an important role in the generation of the baryon asymmetry, leading to an enhancement of the asymmetry by a factor up to $\sim 300$ with respect to the"unflavoured"leptogenesis scenario.


Introduction
Understanding the origin of the excess of matter over antimatter -the matter-antimatter or baryon asymmetry -in the Universe remains one of the fundamental problems in particle physics and cosmology. The asymmetry can be parametrised by the baryon-to-photon ratio, η B , which is defined as where n B , nB and n γ are the number densities of baryons, anti-baryons and photons, respectively. The value of η B can be determined using the data on the Cosmic Microwave Background (CMB) radiation [1]: η B CMB = (6.02 − 6.18) × 10 −10 , 95% C.L.
A very attractive mechanism of generation of the baryon asymmetry is leptogenesis associated with the type-I seesaw scenario of neutrino mass generation [2][3][4][5][6][7][8]: it links the existence and smallness of neutrino masses to the existence of the baryon asymmetry. Integral to this mechanism are the RH neutrinos ν lR (RH neutrino fields ν lR (x)). They can be added to the Standard Model (SM) as SU (2) L × U (1) Y W singlets without modifying any of the fundamental features of the SM. The minimally extended SM with two RH neutrinos is the minimal scheme in which leptogenesis can be realised. The RH neutrinos are assumed to possess a Majorana mass term as well as Yukawa type coupling L Y (x) with the Standard Model lepton and Higgs doublets, ψ lL (x) and Φ(x), respectively. In the basis in which the Majorana mass matrix of RH neutrinos and the charged lepton mass matrix are diagonal, L Y (x) and the Majorana mass term have the form: where Y li is the matrix of neutrino Yukawa couplings (in the chosen basis) and N i (N i (x)) is the heavy Majorana neutrino 1 (field) possessing a mass M i > 0.
In what follows we will consider the "freeze-out" and "freeze-in" flavoured leptogenesis scenarios in which the Yukawa couplings in Eq. (3) are not CP conserving and the different rates of the decays of the Majorana neutrinos N j , N j ↔ l + + Φ (−) , N j ↔ l − + Φ (+) , and of the Higgs boson, Φ (−) → l − + N j , Φ (+) → l + + N j , generate CP violating (CPV) asymmetries in the individual lepton charges L l , and in the total lepton charge L, of the Universe. These lepton asymmetries are converted into a baryon asymmetry of the Universe (BAU) by (B −L) conserving, but (B + L) violating, sphaleron processes which exist in the SM and are effective at temperatures T ∼ (132 − 10 12 ) GeV.
The scale and spectrum of masses of the Majorana neutrinos N j determine the scale of leptogenesis. In GUT scale leptogenesis N j have masses a few to several orders smaller than the scale of unification of the electroweak and strong interactions, M GU T ∼ = 2 × 10 16 GeV. If the heavy neutrinos N j have hierarchical spectrum, M 1 M 2 M 3 , the observed baryon asymmetry can be reproduced provided the mass of the lightest one satisfies M 1 ∼ > 10 9 GeV [9]. Moreover, quantitative studies in which flavour effects in leptogenesis [10][11][12] (see also [13][14][15]) were taken into account have shown that the CP violation necessary for the generation of the observed baryon asymmetry can be provided exclusively by the Dirac and/or Majorana phases in the Pontecorvo, Maki, Nakagawa, Sakata (PMNS) neutrino (lepton) mixing matrix U PMNS [16][17][18][19][20][21][22][23]. More recent analyses revealed [24,25] that in the case of heavy Majorana neutrino mass spectrum with mild hierarchy, M 2 ∼ 3M 1 , M 3 ∼ 3M 2 , i) successful leptogenesis can take place for M 1 ∼ > 10 6 GeV, and that ii) also in this case the required CP violation can be provided exclusively by the Dirac or Majorana CPV phases of the neutrino mixing matrix. In [26] this was confirmed to be the case as well in the so-called "Neutrino Option" seesaw scenario [27] in which the mass term in the Higgs potential, responsible for the electroweak symmetry breaking in the Standard Theory, is generated at one loop level by the neutrino Yukawa coupling in Eq. (3). In the "Neutrino Option" scenario with two Majorana neutrinos N 1,2 , successful leptogenesis was shown to be possible only in the socalled "resonant regime" [28][29][30], with N 1,2 forming a pseudo-Dirac pair [31,32] with masses M ≡ (M 1 + M 2 )/2 ∼ (1 − 8) × 10 6 GeV and splitting between them, which is of the order of the N 1,2 decay widths Γ 1,2 : ∆M/Γ 1,2 ∼ 1, ∆M/M ≡ (M 2 − M 1 )/M ∼ 10 −8 .
One attractive feature of Resonant Leptogenesis (RL) is that the baryon asymmetry can be produced at relatively low scales, e.g., at the TeV scale. Studies have shown that it is possible to have successful RL at scales exceeding approximately 100 GeV (see, e.g., [33] and references quoted therein) or even at smaller scales if thermal effects are taken into account leading to the possibility of CPV Higgs decays into N 1,2 plus a lepton [34]. Scenarios with low scale RL typically lead to predictions that potentially can be tested at colliders (LHC or future planned) and/or at low-energy experiments (see, e.g., [33,34]).
In the present article we consider sub-TeV scale flavoured RL within the minimal type-I seesaw scenario with two (RH) singlet neutrinos N 1,2 forming a pseudo-Dirac pair. We concentrate on the case when the masses of the pseudo-Dirac pair have values M 1,2 ∼ < 100 GeV and consider scenarios in which the baryon asymmetry is generated in CP violating Higgs and N 1,2 decays. We do not consider in this study the "freeze-in" scenario [35,36] in which the BAU in leptogenesis is generated via N 1 ↔ N 2 oscillations during the epoch when N 1,2 are out of equilibrium, which has been extensively studied (see, e.g., [37][38][39][40][41][42][43][44][45][46] and references quoted therein) 2 .
Our study differs from the study performed in [34], in which the Higgs decay scenario has first been considered, in two aspects: i) we use the improvements in accounting for thermal effects in the processes of interest, which have appeared in the literature since the basic study [47] the results of which were employed in [34], and ii) we take into account the flavour effects, which were not accounted for in [34]. We find, in particular, that flavour effects can play very important role in RL of interest especially in the case of zero initial abundance of the heavy Majorana neutrinos.
Our work differs also from the studies reported in [48,49]. In [48] the authors investigated the generation of the baryon asymmetry by both the Higgs decay and the "freeze-in" oscillation mechanisms without separating the contributions of each of the two mechanisms. The aim of our work is, in particular, to identify just the Higgs decay contribution to the generation of BAU. In [48] flavour effects are taken into account, but the enhancement of the range of the Higgs decay contribution due to thermal effects (the collinear emissions of soft gauge bosons (see further)), is not. In our study the latter effects are accounted for. We have 2 In the "freeze-in" N1 −N2 oscillation scenario [35] the generation of the baryon asymmetry proceeds, as was shown in [36], principally via a lepton number conserving (LNC) terms involving fourth power of the neutrino Yukawa couplings. In the scenario considered by us the baryon asymmetry is generated predominantly by lepton number violating (LNV) terms (involving also fourth power of the neutrino Yukawa couplings). considered also two different types of initial conditions for the two heavy Majorana neutrinos: thermal initial abundance (TIA) and vanishing (zero) initial abundance (VIA). In [48] only the case of VIA has been investigated.
In [49] the authors also studied the generation of the baryon asymmetry by both the Higgs decay and the "freeze-in" oscillation mechanisms including the flavour effects, but without separating the contributions of each of the two mechanisms. Only the VIA case has been considered. In addition the masses of N 1,2 were constrained to lie in the interval M 1,2 = 5−50 GeV. In our study the minimal M 1,2 is determined by the requirement of reproducing the observed value of BAU varying all other parameters in the problem. We find, in particular, that in the VIA case the minimal M 1,2 is by more than a factor of 10 smaller than the value assumed in [49].
In the relatively recent article [50] the authors presented a unified framework of generation of the baryon asymmetry at the sub-TeV scale via the neutrino oscillations ("freeze-in") and the decay RL mechanisms within the scenario with two Majorana neutrinos N 1,2 . They have shown that i) the observed baryon asymmetry can be generated for all experimentally allowed values of the Majorana neutrino masses M 1,2 100 MeV and up to M 1,2 ∼ 1 TeV, and that ii) leptogenesis is effective in a broad range of the relevant parameters, including mass splitting between the two Majorana neutrinos as big as (M 2 − M 1 )/(0.5(M 2 + M 1 )) ∼ 0.1, as well as couplings of N 1,2 in the weak charged lepton current large enough to be accessible to planned intensity experiments or future colliders. The results are presented in [50] without separation of the contributions of the "freeze-in" neutrino oscillation mechanism and the decay mechanism in the VIA case. In the TIA case the separation of interest is made only in Fig. 2 and we find similar lower bound on M 1,2 as reported in [50], with the region of the parameter space of successful RL reported in [50] being somewhat larger than that we found in our work. As we have already emphasised, we have concentrated in our work on the decay mechanism of baryon asymmetry generation. We wanted to identify the parameter space in which the observed value of BAU could be generated via the decay mechanism. Thus, effectively we explore part of the parameter space explored in [50].
We should add finally that in [48][49][50] the authors use the density matrix formalism for the calculation of the baryon asymmetry, while we use a system of Boltzmann equations with a proper source term for the treatment of the resonance effects in leptogenesis based on the decay mechanism and take into account the flavour and thermal effects in the decay mechanism. In the strong wash-out regime, in which the results of our study are obtained, we expect the results based on the two approaches and obtained within the decay mechanism to be largely compatible, with differences that should not exceed a factor of ∼ (2 − 3) in the value of the generated baryon asymmetry. For the reasons explained earlier it is impossible at present to make a comparison between our results and those obtained in [48][49][50] in the VIA case; in the TIA case a comparison is possible to a certain degree with the results reported in [50] (see Section 3.3).
The paper is organised as follows. In Section 2, we summarise the basics of the type I seesaw scenario and the conventions we will employ throughout. In Section 3 we introduce the equations relevant for RL at scales T 10 3 GeV of interest. Then we proceed to show results in two possible scenarios by which the BAU can be produced. They correspond to two different "initial conditions", i.e., N 1,2 initial abundances: i) N 1,2 thermal initial abundance (TIA), and ii) N 1,2 vanishing (zero) initial abundance (VIA). We conclude in Section 4 with a brief summary of our results.

Seesaw, Neutrino Masses and Mixing
In the present Section we set the notations and review some of the elements of the seesaw theory that will be used in our further analysis (see, e.g., [51]).
In the basis in which the charged lepton Yukawa couplings and mass matrix are diagonal but the Majorana mass term of the RH neutrinos ν lR is not, the Lagrangian L Y,M (x) has the form: whereỸ is the matrix of neutrino Yukawa couplings in the considered basis, (ψ lL (x)) T = (ν T lL (x) l T L (x)), l = e, µ, τ , ν lL (x) and l L (x) being the left-handed (LH) flavour neutrino and charged lepton fields, where ν c βR ≡ C(ν βL ) T and α, β = e, µ, τ ; in the case of three right-handed neutrinos we can choose 3 κ, ρ = e, µ, τ . The two matrices M D and M N are complex, in general.
The diagonalisation of the mass term under the condition that M D is much smaller than M N 4 leads to the well-known effective Majorana mass (term) for the LH flavour neutrinos [4][5][6][7][8]: wherem ν = diag(m 1 , m 2 , m 3 ), m 1,2,3 being the masses of the light Majorana neutrinos ν 1,2,3 , m i ∼ < 0.5 eV, and U is a 3 × 3 unitary matrix. The flavour neutrino fields are related to the fields of light and heavy neutrinos ν i (x) and N j (x) with definite mass m i and M j , Here ν jL (x) and N jL (x) are the left-handed components of ν i (x) and  [56]. Notice that the 3σ range for the Dirac phase δ is quite large, so we will treat it as an unmeasured parameter.
on flavour observables [52,53]. For M j ∼ > 500 MeV, depending on the element of η, they are in the range of 10 −3 − 10 −4 at 2σ C.L. For M j larger than the electroweak scale, the constraint on η eµ = η µe is even stronger: The PMNS matrix (in the diagonal charged lepton mass basis employed by us) has the form: Thus, the matrix η parametrises the departure from unitarity of the PMNS matrix. Given the existing limits on the elements of η, we have to a very good approximation: where c ij ≡ cos θ ij , s ij ≡ sin θ ij , δ is the Dirac CP violation (CPV) phase, and α 21 and α 31 are the two Majorana CPV phases [55].
As is well-known, the mass spectrum of neutrinos ν 1,2,3 can be with normal ordering (NO), m 1 < m 2 < m 3 , or with inverted ordering (IO), m 3 < m 1 < m 2 (see, e.g., [54]). In what follows we will concentrate on the case of NO neutrino mass spectrum.
As we have already indicated, we will consider the type-I seesaw scenario with only two "heavy" (singlet) Majorana neutrinos N 1,2 . This is the minimal scenario compatible with the oscillation data [54]. In this case the lightest of the three neutrinos ν 1,2,3 is massless at tree and one-loop level, m 1 ∼ = 0 (NO spectrum) and we have: The neutrino mass spectrum is normal hierarchical (NH): m 1 m 2 m 3 . Of the two Majorana phases, α 21 and α 31 , only the phase difference α 21 −α 31 ≡ α 23 , is physical 5 . Technically, we will use the formalism employed for the presence of three heavy Majorana neutrinos N 1,2,3 in which, however, m 1 = 0 and N 3 is decoupled.
In our numerical analyses we will use the values of the three neutrino mixing angles θ 12 , θ 23 and θ 13 , and the two neutrino mass squared differences obtained in the global neutrino oscillation data analysis performed in [56] and quoted in Table 1.
Equation (6) allows to relate the matrix of the neutrino Yukawa couplings Y and the PMNS matrix U [57]. In the diagonal mass basis of the heavy Majorana neutrinos, which is convenient to use in the leptogenesis analyses, we have: where O is a complex orthogonal matrix, In the case of interest with two active "heavy" Majorana neutrinos N 1,2 and the formalism employed by us, we havê M = diag(M 1 , M 2 , M 3 ), but with m 1 = 0 and the form of the O matrix given for the NH spectrum below, the mass M 3 of the third decoupled heavy Majorana neutrino does not play any role in our analysis: where θ = ω + iξ, ω and ξ being two real parameters 6 . The parameters ω and ξ play important roles in the leptogenesis scenario we are going to consider next. For large values of ξ, such that e ξ e −ξ the first term of the above expression dominates, being enhanced by the exponential. The RV -matrix too is enhanced by large values of ξ and the sum of the square modulus of its entries reads (NH): The O-matrices in Eqs. (11) has det(O) = 1. In the literature, the factor ϕ = ±1 is sometimes included in the definition of O to allow for the both cases det(O) = ±1. We choose instead to work with the matrix in Eqs. (11) but extend the range of the Majorana phases α 21 (31) from [0, 2π] to [0, 4π], which effectively accounts for the case of det(O) = − 1 [21], so that the same full set of O and Yukawa matrices is considered.

Flavoured Resonant Leptogenesis at sub-100 GeV scales
In this study we solve the Boltzmann system of equations for the N 1,2 and BAU abundances taking into account both 1 ↔ 2 decays and inverse decays and 2 ↔ 2 scatterings including the thermal effects. In the case of three-flavoured RL of interest it has the form (see, e.g., [10,58,59] where z ≡ M/T , M 1,2 ∼ = M . The quantities N N i and N ∆α are respectively the number of heavy neutrinos N i and the value of the asymmetry ∆ α ≡ 1 3 B −L α , α = e, µ, τ , in a comoving volume, normalised to contain one photon at z = 0, i.e., N eq N i (0) = 3/4. This normalisation within the Boltzmann statistics is equivalent to using N eq N i (z) = 3 8 z 2 K 2 (z), where K n (z), n = 1, 2, .., are the modified n th Bessel functions of the second kind.
Equations (13) and (14) are an excellent approximation to the density matrix equations of [58], which reduce to Eqs. (13) and (14) for T 10 7 GeV (the case in which all lepton flavours are fully decohered). We note that they do not include relativistic corrections and also they are written under the assumption of kinetic equilibrium.
The term p iα which multiplies the wash-outs, is the projection probability of heavy neutrino mass state i on to flavour state α and is given by The projection probabilities p iα , i = 1, 2, α = e, µ, τ , are strongly flavour dependent. We use the following conversion from the B − L number density to the baryon-to-photon ratio: where 28/79 is the SM sphaleron conversion coefficient and the 1/27 factor comes from the dilution of the baryon asymmetry by photons.

The Decay and Scattering Terms
The terms D i and W D i are due to the 1 ↔ 2 decays and inverse decays, while account for scattering processes. The terms S ) are contributions respectively from t-and s-channel 2 ↔ 2 scattering processes (∆L = 1) involving the SM gauge fields [47]. Similarly, S For the total contributions of the 2 ↔ 2 processes involving the SM gauge fields and the top quark to the production of N i and to the wash-out terms we get from Eqs. (17) - (20): We take into account the thermal effects in the production of N 1,2 in Eqs. (13) and (14) using the results derived in [62] for D i , S (gauge) i and S (quark) i in the relevant case of relativistic N 1,2 . As was shown in [62], the indicated three contributions vary little in the interval of temperatures of interest, T ∼ (100 − 1000) GeV, and we have approximated them as constants equal to their respective average values in this interval. Adapting the results obtained in [62] to the set-up utilised by us we get: where the parameter κ i is the ratio between the total decay rate of N i at zero temperature, and the Hubble expansion rate (H) at z = 1. It proves convenient to cast κ i in the form: with m * = (8π 2 v 2 /3M p ) (g * π)/5 ≈ 10 −3 eV, M P being the Planck mass.
Using the generic relations for the wash-out terms, we get: The sum of the three terms is compatible with the result obtained in [63]. We stress that Eqs. (25) - (27) and (32) - (34) are valid only for z < 1. Moreover, as z 2 K 2 (z) 2 for z 1, all the terms given in Eqs. (25) - (27) and (32) - (34) are basically constant at z 1. The behaviour at z > 1 is not relevant for our analysis as we are considering M ≤ 100 GeV, so that the production of the baryon asymmetry stops at z sph < 1.

The CP Violating Asymmetry
We define the CP violating (CPV) asymmetry with the inclusion of thermal effects as in [34], but taking also into account the flavour effects, namely [29,30,33] 8 : where In Eq. (35) the quantity x (0) ≡ ∆M (0) /Γ 22 , ∆M (0) being the N 2 − N 1 mass splitting at zero temperature 9 . The term that multiplies sgn(M i − M j ) I ij,αα in Eq. (35) is due to heavyneutrino mixing effects 10 . Thermal corrections to the N 2 − N 1 mass splitting, ∆M T , with the total mass splitting given by ∆M = ∆M (0) + ∆M T , are relevant in the denominator of the expression for (i) αα only and are accounted for by the term x T (z) [34]: (35) quantifies the thermal effects to the N 1,2 self-energy cut [34] and is determined by Here p and q are the charged lepton and heavy Majorana neutrino four-momenta respectively, L is defined in [64] and the angular brackets indicate a thermal average. The function γ(z) depends on the masses of the Higgs boson M H (T ), charged leptons M l (T ) and heavy Majorana neutrinos M 1 ∼ = M 2 = M (it does not depend on the mass splitting ∆M ). At T > T EW ∼ = 160 GeV, T EW being the temperature at which the electroweak phase transtion (EWPT) sets in [65], the Universe is in the symmetric phase and the Higgs vacuum expectation value is zero. The Higgs boson and the charged leptons possess only thermal masses 11 . For the charged lepton thermal masses we use the expression given in [66]: where g = 0.65 and g = 0.35 are respectively the Standard Model SU (2) L and U (1) Y W gauge coupling constants. At T < T EW , the Higgs VEV v(T ) grows approximately as (see, e.g., [34,65] GeV is the VEV value at zero temperature. Correspondingly, the charged lepton mass M l receives a non-zero contribution M l (v(T )) in the interval T sph ≤ T < T EW due to v(T ) = 0: M 2 l = M 2 l (v(T )) + M 2 l (T ), l = e, µ, τ . The EWPT contribution under discussion M l (v(T )) is proportional to the zero temperature experimentally determined mass m l of the charged lepton l: M l (v(T )) = m l v(T )/v. It is not difficult to convince oneself that for T in the interval T sph ≤ T < T EW 9 The first term in the numerator of the expression in Eq. (36), as can be shown, is lepton number violating (LNV), while the second term is lepton number conserving (LNC). Our numerical analyses show that under the conditions of the "freeze-out" leptogenesis mechanism we are studying, the dominant contribution in the generation of the baryon asymmetry compatible with the observations is given by the LNV term, with the LNC term giving typically a subdominant contribution; in certain specific cases the LNC contribution is of the order of, but never exceeds, the LNV one. 10 If the CPV asymmetry (i) αα is produced in N1,2 decays one has to include in αα also the effect of heavyneutrino oscillations [30,33] (which is different from the the effect dominating the oscillation leptogenesis scenario [35,36]). Here we consider the case when the CPV asymmetry (i) αα is generated in the Higgs decays, Φ (−) → l − + Nj, Φ (+) → l + + Nj. 11 The thermal corrections to the masses of the heavy Majorana neutrinos N1,2 are negligibly small [47,66].
one has M 2 l (v(T )) M 2 l (T ). Thus, for T ≥ T sph of interest the charged lepton masses are given by their thermal contributions specified in Eq. (39).
For the Higgs mass M H (T ) we employ the results obtained from the thermal effective potential given, e.g., in [67,68], which takes into account the effects of the EWPT in the interval of temperatures T sph T T EW , T sph = 131.7 GeV being the sphaleron decoupling temperature [65]. The discussion of the behaviour of M H (T ) in the interval of temperatures of interest is outside the scope of the present study; it can be inferred from the aforementioned EWPT effective potential. We give here only the expression for the thermal contribution to the Higgs mass: Here It is easy to check using M l (T ) and M H (T ) given in the preceding equation that for the values of z lying approximately in the interval 0.34 z 0.93 the Higgs and heavy Majorana decay processes are kinematically forbidden [34,66]. For z 0.34 the only allowed processes are Higgs decays to heavy Majorana neutrinos and charged leptons, whereas at larger z ∼ > 0.93 only the heavy Majorana neutrino decays are allowed. This is reflected in the behavior of the function γ(z) shown in Fig. 1: for 0.34 z 0.93 one has γ(z) = 0. We find also that for z 1, γ(z) ∼ = 23.5 and for z 1, γ(z) ∼ = 1. For T > T sph of interest, the charged lepton masses satisfy M L (T ) > 0.296T sph ∼ = 39 GeV. This implies that in order for the decays of the heavy Majorana neutrinos N 1,2 to be in principle kinematically possible at T > T sph , the masses of N 1,2 must satisfy M 1,2 > M L (T sph ) ∼ = 39 GeV. Taking into account also the Higgs mass leads obviously to a larger lower bound on M 1,2 (see, e.g., [34]).
We note further that in the interval of temperatures T sph ≤ T < T EW , the thermal contribution to the masses of the SM top quark, Higgs, W ± and Z bosons are all of the order ofgT ,g being one of the SM couplings g, g , h t and λ, and that the contribution to these masses of the non-zero temperature dependent Higgs vacuum expectation value v(T ), determined earlier, is of the same order. At the same time, the momenta of the particles in the thermal bath are of the order of πT and are much larger than the masses, so all the particles relevant for our discussion are ultrarelativistic at the temperatures on interest [69]. This implies also, in particular, that the expressions for the decay and scattering terms introduced in the preceding subsection are valid actually for T ≥ T sph ∼ = 131.7 GeV.
When collinear emissions of soft gauge bosons, present in the thermal bath of the Universe at the epoch of interest, are also included in the decay processes, the disallowed region discussed earlier becomes accessible to the Higgs decays due to the increased range of kinematic possibilities [66]. In the calculations performed by us we estimate the effects of these emissions by adding an interpolation of γ(z) across the "gap" interval 0.34 z 0.93 12 . This is 12 The effects of collinear emissions of soft gauge bosons are included also in the decay terms Di and W D i given in Eqs. (25), (31) and (32), as discussed in the preceding section. illustrated in Fig. 1 showing γ(z) as a function of z. For a given value of M , the behaviour of γ(z) at z ≥ z sph ≡ M/T sph (or at T ≤ T sph ) is not relevant and should be ignored since the sphalerons decouple at T sph .
We note that at high temperatures γ(z) is sensitive to precise values of Higgs and charged lepton thermal masses. As we have indicated, in our numerical analysis we use Higgs and charged lepton thermal masses in agreement with that of [66]. We estimate that a different choice of thermal masses, as well as a different interpolation of γ(z), would only affect our final results by an insignificant factor.
In the absence of thermal effects, namely setting γ(z) = 1 and x T (z) = 0, the second term in Eq. (35) is maximised for x (0) 0.5. This is the "resonant" behaviour typical to RL without thermal corrections. When thermal corrections are taken into account, it is not possible to choose one value of x (0) for which we have resonance at all temperatures. For this reason, the results we will show further are given for different values of x (0) .
For e |ξ| e −|ξ| , the term involving the Yukawa couplings given in Eq. (36) is proportional to sin(2ω)e −2|ξ| . Hence, large values of |ξ| suppress the CPV asymmetry, while ω = (2n+1)π/4 maximises it (in absolute value). For |ξ| 1, a slightly different dependence on ω appears in both the Yukawa coupling term (in the denominator) and in the mixing term of Eq. (35), so that the maximal value of the CPV asymmetry is actually reached for different values of ω (depending on ξ). We find, however, that a more precise choice of ω would not lead to significant differences in the BAU. Therefore, in obtaining our results we set ω = π/4 or 3π/4 (to match the sign of the BAU) in order to maximise the CPV asymmetry at large values of |ξ| 13 . In the analysis which follows we consider only values of ξ ≥ 0 since the results are symmetric for the corresponding ξ < 0.
Throughout we use the ULYSSES [70] Python package for numerical solutions of the Boltzmann equations. We will present results of the numerical analysis for the Majorana phase α 23 set to zero, adding comments in certain cases on how the results change for α 23 having a non zero value in the interval 0 < α 23 ≤ 4π.

Thermal Initial Abundance
Consider the case of thermal initial abundance (TIA) of the heavy Majorana neutrinos, . We can set the ratio N N i /N eq N i = 1 in the r.h.s. of Eq. (14) since under the indicated initial condition the deviations of N N i from N eq N i for any z > z 0 of interest for our analysis are sufficiently small and can be neglected. For the sum of three wash-out factors, W i , in this case we get: The flavoured washout terms in Eq. (14) are given by W iα ≡ p iα W TIA i . Due to the projection probability p iα the wash-out terms exhibit strong flavour dependence.
In what follows we discuss the results of our numerical analysis. Figure 2 shows the evolution of the lepton and baryon asymmetries, the difference |N eq N i − N N i |, the decay, scattering and wash-out rates for δ = 3π/2, M = 10 GeV, x (0) = 100 and ξ = 2.05 -the maximal value of ξ for which we can have successful leptogenesis for M = 10 GeV -and x (0) = 100. The baryon asymmetry η Bl originates from the CPV asymmetry in the lepton charge (flavour) L l (l), l = e, µ, τ . Thus, the total baryon asymmetry is η B = η Be + η Bµ + η Bτ . The figure illustrates the typical scenario of "freeze-out" leptogenesis, namely, the case when the departure from equilibrium of N N i (z) is what drives the generation of the lepton (and baryon) asymmetry. The total baryon asymmetry, to which all flavour CPV asymmetries contribute, freezes at z sph = M/T sph ∼ = 0.076, T sph = 131.7 GeV being the sphaleron decoupling temperature, which is marked by the vertical grey line in the figure.
For the choice of parameters in the figure, the asymmetries |η Bµ | and |η Bτ | exhibit almost identical evolution for z ≤ z sph and are by a factor ∼ 100 larger than |η Be |. This difference reflects the difference between ee . Thus, in this case, η B ∼ = η Bµ + η Bτ . The fact that |η Bµ | ∼ = |η Bτ | |η Be | can have important implications in what concerns the possibility of wash-out of the baryon asymmetry by lepton number non-conserving effective operators of dimension higher than four that might be "active" at the energy scales of interest [71,72]. For z > z sph , and therefore after sphaleron freeze-out, the asymmetry |η Be | converges to the asymmetries |η Bµ | and |η Bτ |. Qualitatively similar behaviour is seen for a range of x (0) and M values, with the main difference being the overall scale of the asymmetry evolution.
To understand the impact of flavour effects, we compare the obtained results with the results in the unflavoured case. The unflavoured approximation is equivalent to taking in Eq. (14) p iα = 1 for every α and then sum over all the flavours. It roughly corresponds to taking the total asymmetry to be the asymmetry in the dominant flavour (either muon or tauon in  Figure 2: The evolution of the lepton flavour and baryon asymmetries, of |N eq N i − N N i | and of the corresponding decay, scattering and washout rates that govern the evolution of the asymmetries in Eqs. (13) and (14) in the case of TIA of N 1,2 . The figure is obtained for δ = 3π/2, M = 10 GeV, x (0) = 100 and ξ = 2.05. The vertical grey line is at z sph = 0.076 and is the endpoint of evolution of the baryon asymmetry η B . The horizontal grey line at 2 is roughly indicating where the different processes get into equilibrium. See text for further details.
in the considered case). As it is then clear from Fig. 2, this approximation would only lead to a O(2) difference in the value of η B . A more detailed analysis shows that flavour effects in the TIA case lead, in general, to a moderate enhancement by a factor of ∼ (2 − 3) of the baryon asymmetry.
We show in Fig. 3 the maximal values of ξ, for which we can have successful RL as a function of the mass scale M . We recall that to the maximal values of ξ there correspond maximal values of |(RV )| 2 (see eq. (12)). The curves of different colours correspond to the contours at different x (0) for which the maximal baryon asymmetry coincides with the present observed value of ≈ 6 × 10 −10 . In the blue region, the baryon asymmetry is too small compared to the observed value. In the white region instead, we can always vary x (0) , ω and/or the CPV phases in the PMNS matrix so to get the correct value for η B .
In the considered scenario, for a given x (0) , the lower the mass, the less is the time for the system to depart from equilibrium before z sph , and so greater must be the CPV asymmetry, and slower the processes keeping N i in equilibrium, so that the baryon asymmetry freezes at the observed value. Correspondingly, by lowering the mass, the maximal value of ξ for which we can have successful RL decreases. This leads to a lower bound on the mass M that depends on x (0) , for which RL can be successful. In the region of small ξ the dependence on ξ of the CPV asymmetry and wash-out terms is less trivial and strongly dependent on the leptonic CPV phases, leading to the growth with the mass, as shown in Fig. 3.
As long as x (0) x T (z), the CPV asymmetry grows with z and x (0) , whilst for x (0) x T (z)  Fig. 3 with those derived for α 23 = 0 we do not find any significant dependence on the value of the Majorana phase α 23 . We can conclude, in particular, that: i) the minimal lower bound on the mass M for having successful flavoured RL is practically the same, namely, M 5 GeV, for all values of α 23 considered; ii) the maximal value of ξ is the same as found for α 23 = 0, namely, ξ ∼ = 3.3, and similarly to the α 23 = 0 case, is reached for M ∼ = 100 GeV; iii) for the maximal value of ξ reached at M = 20 GeV we find ξ = 2.63 at α 23 = 360 • ; it is somewhat larger than the value of ξ = 2.55 found for α 23 = 0 and corresponds to max( l,i |(RV ) li | 2 ) 2.8 × 10 −10 .
Our detailed analysis showed that the wash-out factor for the asymmetry in the e-lepton charge and the related baryon asymmetry η Be exhibit weak dependence of α 23 , while the wash-out factors for the asymmetries in the µ-and τ -lepton charges and, correspondingly, the flavour baryon asymmetries η Bµ and η Bτ , change significantly with α 23 . However these changes are "anti-correlated" in the sense that the sum η Bµ +η Bτ remains practically constant when α 23 is varied. As a consequence, also the total baryon asymmetry η B = η Be + η Bµ + η Bτ does not show any noticeable dependence of α 23 .
We compare next briefly the results obtained by us with those derived in [50] in the TIA case. The results in the decay scenario in the TIA case obtained in [50] are indicated graphically in Fig. 2 in [50]. The minimal value of M 5 GeV for successful leptogenesis we find in our analysis is similar to the value of ∼ 7 GeV found in [50]. In what concerns l,i |(RV ) li | 2 (which is equivalent to the parameter |U | 2 in Fig. 2 in [50]), for the maximal value of this important for the phenomenology parameter at M = 10, 20 and 100 GeV we find, respectively, max( l,i |(RV ) li | 2 ) 2.1 × 10 −10 , 2.8 × 10 −10 and 2.2 × 10 −10 , while, according to their Fig. 2, the authors of [50] get ∼ 5 × 10 −10 , ∼ 1.0 × 10 −9 and ∼ 1.2 × 10 −9 . Thus, our results for the maximal value of the observable ( l,i |(RV ) li | 2 ) at M = 10, 20 and 100 GeV are somewhat smaller that those obtained in [50] and thus can be considered as more conservative.
For completeness we summarise our results for the observables max( l,i |(RV ) li | 2 ) and min( l,i |(RV

Vanishing Initial Abundance
We analyse next the case of a vanishing initial abundance (VIA) of N 1,2 , i.e. N 1,2 (z 0 ) = 0. For the wash-out terms, we assume on the basis of the results reported in [47] that Under these conditions 14 the wash-out term in Eq. (14) in the case of interest has the form:   Figure 4: The evolution of the lepton-flavour and baryon asymmetries, of |N eq N i − N N i | and of the corresponding decay, scattering and washout rates that govern the evolution of the asymmetries in Eqs. (13) and (14) in the case of vanishing initial abundance (VIA) of N 1,2 . The figure is obtained for δ = 0, M = 15 GeV, x (0) = 10 6 and ξ = 1.53. The vertical grey line at z sph = 0.114 is the endpoint of the evolution of the baryon asymmetry η B . The horizontal grey line at 2 is roughly indicating where the different processes get into equilibrium. At z z sph , the asymmetry |η Be | (solid red curve) is smaller in magnitude than the asymmetries |η Bµ | and |η Bτ |. However, the weaker washout of η Be results in it dominating η Bµ and η Bτ by the time of sphaleron decoupling and in η B ∼ = η Be . See text for further details.
The flavoured washout terms in Eq. (14) are given by W iα ≡ p iα W VIA i . In Figs. 4, 5 and 6 we show the evolution of the leptonic asymmetries N ∆α respectively for i) δ = 0, M = 15 GeV (z sph = 0.114), x (0) = 10 6 and maximal ξ = 1.53, ii) δ = 300 • , M = 20 GeV (z sph = 0.15), x (0) = 10 6 and ξ = 1.33, and iii) δ = 3π/2, M = 16 GeV (z sph = 0.121), x (0) = 10 3 and maximal ξ = 2.18, respectively. Also shown is the growth of N N towards the evolving equilibrium distribution N eq N (z), governed by the combination D i + S t i + S s i . The sphaleron transition occurs at z sph marked by the vertical grey line after which the baryon asymmetry η B is "frozen" and remains constant at the value at z sph . Sharp dips in the asymmetries correspond to sign changes as we always plot absolute values.
As is seen in Fig. 4, the asymmetries η Bµ and η Bτ , which are generated by the CPV µand τ -flavour asymmetries, are strongly suppressed in the interval 0.07 z ≤ z sph due to the relatively large wash-out factors. This is reflected in the sudden dips of the corresponding curves as they are driven through zero by the wash-out effects. As a consequence, by the time of sphaleron decoupling most of the baryon asymmetry is due to the lepton CPV asymmetry residing in the electron flavour, η B ∼ = η Be . Since η Be was mostly generated during the production of RH neutrinos, i.e., before N N i reached N eq N i , this case corresponds to the "freezein" scenario of generation of baryon asymmetry.
We highlight the fact that, in contrast to the TIA case, in the VIA scenario illustrated in The unflavoured approximation would then neglect the electron CPV asymmetry and consequently η Be , which actually contributes most in the flavoured scenario. In this particular case, flavour effects lead to a O(300) enhancement with respect to the unflavoured case 15 .
We find that the enhancement of the baryon asymmetry due to flavour effects depends strongly on the CPV phase δ. This is illustrated in Fig. 5, which shows a "freeze-in" scenario of leptogenesis for δ = 300 • , in which, in contrast to case with δ = 0 reported in Fig. 4, the flavour enhancement of the baryon asymmetry is approximately by a factor of 60. The results presented in Fig. 5 show also that, depending on the values of leptogenesis parameters, all three lepton CPV asymmetries residing in the electron, muon and tauon flavours can give significant contributions to the baryon asymmetry so that η B = η Be + η Bµ + η Bτ .
The features reported in the preceding discussion can be obtained for other choices of the parameters and we find, in general, that in the mass range of interest, varying x (0) and δ accordingly, flavour effects can lead to enhancement of the generated baryon asymmetry by a factor ranging from a few to a few hundred.
In Fig. 6 instead, the final baryon asymmetry is generated after all the three lepton flavour CPV asymmetries, initially generated during the production of RH neutrinos, are fully erased by the washout processes. Therefore this corresponds to the "freeze-out" RL case. Since the dominant lepton flavour related asymmetries are η Bµ and η Bτ and η B ∼ = η Bµ +η Bτ , the flavour effects are not significant in this scenario. However, it is quite remarkable that in the case of the same initial condition -zero initial abundance of N i , the "freeze-in" mechanism of baryon asymmetry generation can transform into "freeze-out" mechanism for different choices of the parameters.
In Fig. 7, we show in the ξ − M plane the region of successful RL in the VIA case. With respect to the TIA case, the region of successful RL extends to lower masses. In particular, we find that the minimal lower bound on the mass is reached for for x (0) ≈ 10 6 and reads M ≈ 0.3 GeV.
As in the TIA case, we do not find any significant dependence of the results shown in Fig.  7 on α 23 when α 23 is varied in the interval (0, 4π]. More specifically, i) the lower bound on M is always at M 0.3 GeV and is reached for x (0) = 10 6 ; ii) max(ξ) 3.3 takes place at M = 100 GeV and x (0) = 1 − 10 for any choice of α 23 ; iii) at M = 1 GeV, max(ξ) 2.5 practically for any α 23 ; iv) at M = 20 GeV, max(ξ) 2.55 (2.63) for α 23 = 0 (2π). We notice that the curves corresponding to fixed values of x (0) 10 4 show some dependence on α 23 . However, this has a minor effect on the full white region in Fig. 7, corresponding to values of the parameters for which we have successful leptogenesis.
The values of the observables max( l,i |(RV ) li | 2 ) and min( l,i |(RV ) li | 2 ) for M = 10, 20, The allowed region (white) in the ξ-M plane for which we can always have successful resonant leptogenesis by varying x (0) , ω and/or the CPV phases in the PMNS matrix, in the case of VIA. For values in the red region the baryon asymmetry is always too small compared to that observed today. The solid contours corresponds to the maximal value of ξ and x (0) = 10 7 (blue), 10 6 (red), 10 5 (green) 10 4 (magenta), 10 3 (orange), 10 2 (cyan), 10 (brown), and 1 (black), for which the predicted asymmetry coincides with the observed one. The figure is obtained for δ = 3π/2. See text for further details. 30,40,50,60,70,80, 85 GeV we find in the VIA case are the same as those reported for the TIA case at the end of Section 3.3. This is in agreement with the fact that in the freeze-out leptogenesis we are analysing, in which the observed BAU is generated in the strong wash-out regime, there is no dependence on the initial conditions. We emphasise that our analysis does not include the N 1 − N 2 oscillation mechanism proposed in [35,36], and thus our results rely purely on 1 ↔ 2 decay and 2 ↔ 2 scattering processes in the RL case.

Conclusions
In the present article we have considered sub-TeV scale flavoured resonant leptogenesis (RL) within the minimal type-I seesaw scenario with two (RH) singlet neutrinos N 1,2 forming a pseudo-Dirac pair. We concentrated on the case when the masses of the pseudo-Dirac pair have values M 1,2 ∼ < 100 GeV, (M 2 − M 1 ) ≡ ∆M M 1,2 , and have considered temperatures in the interval T ∼ (100 − 1000) GeV. The change of the baryon asymmetry η B during the generation process "freezes" at the sphaleron decoupling temperature T sph = 131.7 GeV and the value of η B at this temperature should be compared with the observed one. We have not analysed in this study the scenario [35,36] in which the BAU in leptogenesis is generated via N 1 ↔ N 2 oscillations, which has been extensively studied by many authors.
We have investigated and presented results for two possible scenarios by which the BAU can be produced. They correspond to two different "initial conditions", i.e., N 1,2 initial abundances at T 0 >> T sph : i) N 1,2 thermal initial abundance (TIA), and ii) N 1,2 vanishing (zero) initial abundance (VIA). In our analyses we took into account both the relevant 1 ↔ 2 decays and inverse decays and 2 ↔ 2 scattering processes including the thermal effects.
In the case of TIA the baryon asymmetry is produced via the so-called "freeze-out" mechanism, i.e., by out-of-equilibrium processes when N 1,2 essentially coincides with their respective thermal abundances (Figs. (2) and (3)). We find that in this case successful RL is possible for M 1.2 as low as 5 GeV, M ∼ > 5 GeV. The flavour effects are not so relevant in the generation of the baryon asymmetry leading only to a moderate enhancement approximately by a factor of 2-3 of the asymmetry.
Our results show that in the VIA case one can have a successful leptogenesis for M 1,2 lying in the interval M 1,2 ∼ = (0.3 − 100) GeV (Figs. 4, 5, 6 and 7). For the values of M 1,2 in this interval and depending on the values of the other leptogenesis parameters, the baryon asymmetry can be generated i) either during the production of the heavy Majorana neutrinos N 1,2 corresponding to the "freeze-in" leptogenesis scenario (Figs. 4 and 5), or ii) after the produced N 1,2 abundance reached the thermal equilibrium value, corresponding to the "freeze-out" leptogenesis scenario (Fig. 6). It is quite remarkable that in the case of the same initial condition -vanishing (zero) initial abundance of N i (VIA), one can have a successful leptogenesis by both the "freeze-in" and the "freeze-out" mechanisms, each of the two scenarios being operative in different regions of the relevant parameter space.
We have shown also that in the VIA case the flavour effects, in general, play a particularly important role for having successful leptogenesis. In certain cases they can lead to an enhancement of the baryon asymmetry by a factor of a few hundred (Fig. 4) with respect to asymmetry produced by neglecting the flavour effects, i.e., within the unflavoured leptogenesis scenario. The magnitude of the flavour enhancement exhibits strong dependence on the value of the Dirac CPV phase δ.
Our results show further that, depending on the values of the leptogenesis parameters, the dominant among the three flavour components of the baryon asymmetry, η Be , η Bµ and η Bτ , generated by the corresponding CP violating (CPV) e-, µ-and τ -flavour (lepton charge) asymmetries, could be the η Be , or the sum η Bµ + η Bτ or else the contribution from all three components can be significant. Thus, the total baryon asymmetry, η B can originate either from the CPV e-flavour asymmetry, η B ∼ = η Be , or from the CPV µ-and τ -flavour asymmetries, η B ∼ = η Bµ + η Bτ , or else from all three the CPV flavour asymmetries, η B = η Be + η Bµ + η Bτ .
To summarise, we have shown that RL can be successful across the whole of the experimentally accessible region of M 1,2 ∼ = (0.3 − 100) GeV. Furthermore, we have found that leptogenesis at the considered sub 100 GeV scales is compatible with values of the charged and neutral current couplings of N 1,2 in the weak interaction Lagrangian, whose squares for, e.g., M 12 = (10 − 85) GeV are in the range of 3 × 10 −10 − 10 −12 . This may pose challenges for testing the considered leptogenesis scenario in low-energy experiments. A large part, if not all, of the indicated range, however, can be probed in the experiments at the discussed future FCC-ee facility [73,74].