Leptogenesis from heavy right-handed neutrinos in CPT violating backgrounds

We discuss leptogenesis in a model with heavy right-handed Majorana neutrinos propagating in a constant but otherwise generic CPT-violating axial time-like background (motivated by string theory). At temperatures much higher than the temperature of the electroweak phase transition, we solve approximately, but analytically (using Padé approximants), the corresponding Boltzmann equations, which describe the generation of lepton asymmetry from the tree-level decays of heavy neutrinos into Standard Model leptons. At such temperatures these leptons are effectively massless. The current work completes in a rigorous way a preliminary treatment of the same system, by some of the present authors. In this earlier work, lepton asymmetry was crudely estimated considering the decay of a right-handed neutrino at rest. Our present analysis includes thermal momentum modes for the heavy neutrino and this leads to a total lepton asymmetry which is bigger by a factor of two as compared to the previous estimate. Nevertheless, our current and preliminary results for the freezeout are found to be in agreement (within a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 12.5\%$$\end{document}∼12.5% uncertainty). Our analysis depends on a novel use of Padé approximants to solve the Boltzmann equations and may be more widely useful in cosmology.


Introduction and motivation
A plethora of cosmological measurements, especially those associated with observations of the cosmic microwave background radiation (CMB) in the universe [1,2], lead to an estimate of the observed asymmetry between matter (mostly baryons) and antimatter of order: in the early stages of cosmic expansion, i.e. for times t < 10 −6 s and temperatures T > 1 GeV. In the above formula n B (n B ) denotes the (anti-) baryon density in the universe, and s is the entropy density of the universe. Moreover, the observations of the CMB background, indicate that at present the temperature of the universe is T 0 = 2.727 K = 0.235 meV and the ratio of baryons over photons is where n γ is the density of photons in the universe. At first sight, the asymmetry (1) (and the result (2)) appears to be in conflict with fundamental properties of relativistic quantum field theories, which form the basis of our phenomenology of elementary particles. Specifically, in flat space-time, any unitary and local Lorentz invariant quantum field theory, which respects unitarity and locality, should be described by a Lagrangian that is invariant under CPT transformations where C denotes charge conjugation, T denotes reversal in time and P denotes parity (spatial reflection) transformations. This is the celebrated CPT theorem [3]. For the physics of the early universe based on any Lorentz invariant quantum field theory, such a theorem implies that matter and antimatter should be created in equal amounts after the Big Bang. If such is the case, the universe today would be filled with radiation, as a result of matter-antimatter annihilation processes, in conflict with (2).
Within the context of our current understanding of fundamental physics, A. Sakharov [4][5][6][7], postulated the following three necessary conditions for the dominance of matter over antimatter (baryon asymmetry in the universe (BAU) (1)), and hence for our very existence: • Baryon (B) number violation.
• Charge (C) and charge-parity (CP) symmetries need to be broken.
• Chemical equilibrium does not hold during an epoch in the early universe, since chemical equilibrium washes out asymmetries.
In fact there are two types of non-equilibrium processes in the early universe that can produce asymmetries between particles and antiparticles: the first type concerns processes generating asymmetries between leptons and antileptons (leptogenesis) [8][9][10][11][12], while the second produces asymmetries between baryons and antibaryons directly (baryogenesis) [13][14][15][16]. Unfortunately, within the framework of the Standard Model (SM), although Sakharov's axioms can be qualitatively reproduced, especially because one has both B and CP violation in the quark sector, the resulting baryon asymmetry is several orders of magnitude smaller than the observed one (1) [17][18][19]. There are several ideas that go beyond the SM (e.g. grand unified theories, supersymmetry, extra dimensional models etc.) and provide extra sources of CP violation, necessary for yielding the observed magnitude for the asymmetry. Some of these attempts, involve the elegant mechanism of baryogenesis via leptogenesis, in which a lepton asymmetry is generated first, by means of decays of right handed sterile neutrinos to SM particles; the lepton asymmetry is subsequently communicated to the baryon sector by means of sphaleron processes which violate both Baryon (B) and Lepton (L) numbers, but preserve the difference B-L [20][21][22][23][24][25][26]. Heavy sterile neutrinos, through the seesaw mechanism [27][28][29][30][31], play another essential rôle in particle physics, since they provide a natural explanation for the existence of three light neutrinos with masses small compared to other mass scales in the SM), as suggested by observed neutrino oscillations [32,33]. Fine tuning and some ad hoc assumptions are involved though in such scenarios, especially in connection with the magnitude of the CP violating phases and the associated decay widths. Consequently the quest for a proper understanding of the observed BAU still requires further investigation.
In the scenario of Sakharov it is assumed that CPT symmetry holds in the very early universe and this leads to the equal production of matter and antimatter. CPT invariance is regarded as fundamental since it is a direct consequence of the celebrated CPT theorem [3]. However, it is possible that some of the assumptions in the proof of the CPT theorem do not hold in the early universe, leading to violations of CPT symmetry. Sakharov has stated that non-equilibrium processes are necessary for BAU in CPT invariant theories. If the requirement of CPT is relaxed, the necessity of nonequilibrium processes can be dropped . In a low-energy version of quantum gravity Lorentz invariance and unitarity are likely to emerge since not all degrees of freedom are accessible to a low-energy observer. Lorentz invariance violation has been singled out in Ref. [34] as a fundamental reason for inducing CPT violation (CPTV) and vice versa. (However, such claims have been disputed in [35,36], through counterexamples of Lorentz invariant systems, which violate CPT through relaxation, for example, of locality.) In our work we will consider Lorentz invariance violating (LV) backgrounds in the early universe as a form of spontaneous violation of Lorentz and CPT symmetry.
If LV is the primary source of CPTV, then the latter can be studied within a local effective field theory framework, which is known as the Standard Model Extension (SME) [37]. The latter provides the most general parametrization for studying the phenomenology of Lorentz violation in a plethora of physical systems, ranging from cosmological probes, to particle and precision atomic physics systems. For the current era of the universe [38][39][40][41] very stringent upper bounds on the potential amount of Lorentz and CPT violation have been placed by such systems. However, under the extreme conditions present in the very early universe, such violations could be significantly stronger than in the present era (where they could be extremely suppressed (or absent), in agreement with current stringent constraints). 1 In a previous work [43] we presented a phenomenological model for generating a lepton asymmetry via CPTV in the early universe. The model was based on a specific extension of the SM, involving massive Majorana right-handed neutrinos (RHN), propagating on a Lorentz and CPTV, constant in time, axial vector background coupling to fermions. The latter could be traced back to a specific configuration of a cosmological Kalb-Ramond (KR) antisymmetric tensor field [44] that appears in the gravitational multiplet of string theory [45][46][47][48][49], and plays the rôle of torsion in a generalised connection, although we do not restrict ourselves to such an identification. 2 The involvement of sterile RHN in the model is physically motivated primarily by the need to provide a natural explanation for the light neutrino masses of the SM sector. The lightest RHN may also have a potential role as (warm) dark matter candidates [25,26,[51][52][53]. However, in our CPTV models sterile neutrinos responsible for leptogenesis have masses in the 10 5 GeV range or higher [43] and so cannot be considered as dark matter.
In [43] we only gave a qualitative and rather crude estimate of the induced CPTV lepton asymmetry, based on the decaying right handed Majorana neutrino being at rest. In this way it was possible to estimate the lepton asymmetry, without following the standard procedure of solving the appropriate Boltzmann equation that determines correctly the asymmetry value at decoupling of RHN. In the early universe the heavy right-handed neutrinos are not at rest but have a Maxwell-Boltzmann momentum distribution. The purpose of this article is to properly take into account this momentum distribution in the calculation of the lepton asymmetry.
The structure of the article is as follows: in the next Sect. 2 we review the model of [43] and an earlier estimate of the CPTV-background induced lepton asymmetry, which shall be compared with the much more accurate result of the present article, obtained by solving the appropriate Boltzmann equations analytically. We commence our analysis by considering the lepton asymmetry associated with the decays of the RHN into charged leptons. In Sect. 3, we construct the appropriate system of Boltzmann equations in the presence of a weak CPTV axial background involved in the problem, and compare it with the standard CP violating case [20][21][22][23][24][25][26]. In Sect. 4, we solve the Boltzmann equations using Padé approximants [54], which is an approximation popular in several fields of physics, ranging from statistical mechanics to particle physics and quantum field theory [55][56][57][58][59][60]. In this way, we manage to compute the induced lepton asymmetry at RHN decoupling analytically, avoiding numerical treatment. It should be remarked, that setting up and solving such a system of differential equations is a highly non-trivial and algebraically complicated task. Our analytical results agree (within ∼ 12.5% accuracy) with our earlier preliminary estimates of the freezeout point, as outlined, in [43]. In view of 2 In four space-time dimensions, the field strength of the KR field is dual to a massless pseudoscalar (axion-like) field. In the recent literature [50] axion-based approaches to leptogenesis, involving an effective CPT-violating coupling between the (temporal component of the) lepton number current and the time derivative of a (time-dependent) axion field (which is quite different from our KR axion which couples to axial fermion currents), have been proposed. This interaction breaks time translation invariance and, thus, generates an effective chemical potential which differentiates between leptons and anti-leptons. The presence of this effective chemical potential allows for the generation of a lepton asymmetry by means of RHN-mediated L = 2 scattering processes in that model. this, we consider our system of Boltzmann equations as providing another efficient use of Padè approximants, this time with relevance to cosmology. The lepton asymmetry that we find in our analytic treatment is slightly larger (by a factor of about 2) than the estimate of [43]; this is to be expected, since non-zero momentum modes of the RHN have been included. In Sect. 5 we complete our analysis by including the contributions to the Boltzmann equations and the lepton asymmetry coming from the decays of the RHN into the neutral Higgs and active neutrinos. Our calculations show that the resultant lepton asymmetry increases by a factor ∼ 2 compared to the one based on the RHN decays to charged leptons only. Conclusions and outlook are given in Sect. 6. A review of the formalism and derivations of the corresponding decay amplitudes and thermally averaged rates used in the Boltzmann equations, are presented in several Appendices.

Review of the CPT violating model for leptogenesis
It will suffice for our purposes to consider a single species of RHN as in [43]. If the phenomenology is required to include the seesaw mechanism it is necessary (and possible) to add more species of RHN. The option of using a single species of RHN is not available within the standard CPT conserving but CP violating scenario, where to obtain a lepton asymmetry one needs more than one species of RHN [20,21,24]. Our Lagrangian is given by [43]: where N is the Majorana field,φ is the adjoint (φ i = ε i j φ j ) of the Higgs field φ, and L k is a lepton (doublet) field of the SM sector, with k a generation index. y k is a Yukawa coupling, which is non-zero and provides a non-trivial interaction between the RHN and the SM sectors via the Yukawa type interaction ("Higgs portal"): In our case of a single Majorana neutrino species we take k = 1 to label the first generation, and from now on we set Since in SM the leptons have definite chirality, the Yukawa interactions L Y U K can be rewritten as where in the last equality we used the properties of the charge conjugation matrix and the Majorana condition N c = N . The two hermitian conjugate terms in the Yukawa Lagrangian are also CPT conjugate. This is to be expected on the basis of the CPT theorem. In fact CPT violation is introduced only by interactions with the background field. The background field / B ≡ γ μ B μ is assumed at most a function of the cosmic time, so as to respect the isotropy and homogeneity of the early universe, where such backgrounds are non-trivial. We note at this point that, if the axial background field B μ is to be identified [43] with the totally antisymmetric field strength (H μνρ = ∂ μ B ρσ + cyclic permutation of indices) of the Kalb-Ramond [44] spin-one field B μν , that appears in the massless gravitational multiplet of string theory [45][46][47], then the latter is viewed as part of a torsion background [48]: In such a case one should also consider the coupling of the axial field B μ to all other fermions of the SM sector, ψ j ( j = leptons, quarks) via a universal minimal prescription, with the coupling with all fermionic species ψ being the same : In four space-time dimensions the H νρλ field is dual to a pseudoscalar field b(x) [48,49]: There is an exact cosmological solution in the bosonic string theory [49], in which the H -torsion background is identified with a homogeneous and isotropic cosmological Kalb-Ramond axion, linearly dependent on the cosmic time [49]. The solution satisfies the corresponding conformal invariance conditions of the associated σ -model, thus constituting a consistent background of strings. The resultant axial backgrounds are constant in time and have non-trivial temporal components only In [43] we have generalised the above solution (7) in theories with fermions, in which the latter condensed in the early universe. Such backgrounds can then be viewed as spontaneously breaking Lorentz and CPT symmetry in the system and are consistent with isotropy and homogeneity of the early universe. In what follows we shall consider the Lagrangian (4) in the generic background (7), without specifying further its microscopic origin. The form of the Lagrangian coincides with one of the simplest forms of the so-called Standard Model Extension (SME) [37], namely that in which the temporal component of the so-called b μ coefficient assumes a constant value.
There are stringent constraints [38][39][40][41] (coming from a plethora of measurements ranging from astrophysical to laboratory precision tests of Lorentz and CPT symmetries) for today's value of b 0 ≤ 0.02 eV (and much suppressed spatial components b i < 10 −32 GeV). Although in our model in the frame of Robertson-Walker (cosmic microwave background) the axial background is assumed to have only the temporal component (7), nevertheless the slightest motion of the observer with respect to that frame will generate a spatial component by means of a Lorentz transformation. It is therefore essential that any current value of B 0 is severely suppressed today, and also during the nucleosynthesis era. In [43] we have provided arguments in favour of scenarios in which the universe undergoes a phase transition soon after the decoupling of heavy neutrinos, so that the background B 0 ceases to be a constant, and decreases with the temperature according to the scaling law T 3 . The qualitative estimates of [43], have indicated that for Yukawa couplings y k of order 10 −5 (assumed in [43]), the decoupling temperature of the heavy neutrino T D of order T D m N ∼ 100 TeV, implies a phenomenologically consistent leptogenesis for B 0 ∼ 1 MeV at T T D . The cooling law implies for the present era a negligible B 0 = O(10 −44 ) meV today, and also a very small value during the nucleosynthesis era.
As we shall be interested in high temperatures T T D ∼ 100 TeV, which are much higher than the electroweak phase transition, the SM fields are treated as massless, while the heavy RHN can still be assumed to be massive. 3 In such a case, the Higgs field does not develop a vacuum expectation value; consequently the charged Higgs (denoted by h ± ) and neutral Higgs (h 0 ) play a rôle in the physical spectrum. From the form of the interaction Lagrangian in Eqs. (4), and (6), it is straightforward to obtain the Feynman rules for the diagrams giving the decay of the Majorana particle in the two distinct channels: The neutral channel decay N → ν h 0 , where ν are the SM neutrinos, would not lead to any lepton asymmetry, if the active neutrinos ν were purely Majorana; this follows directly from the Yukawa term (6), when expressed in terms of Majorana fields for the neutrinos. However, in standard see-saw scenarios [27][28][29][30][31], the Lagrangian contains Dirac mass terms for the active neutrinos and schematically ν = ν; so there would be additional contributions to the lepton asymmetry from the tree level decays in the channels In the absence of the background, the squared matrix elements obtained from tree level diagrams for the two decays in Eq. (9) (cf. Fig. 1), and also in Eq. (10), would be the same [20,21,24,62]. In such a case, a lepton asymmetry is generated due to the CP violation present in the neutrino mixing matrix in the pertinent one loop diagrams, and hence require more than one species of right-handed neutrinos. In the presence of the background B 0 = 0, however, there is a difference in the decay rates of the tree level processes (9), and this leads to CPTV-induced lepton asymmetry. 4 In what follows, we shall first calculate the lepton asymmetry based only on the decay channels (9), involving charged leptons in the final stage. In Sect. 5 we shall include the neutral decay channels (10), into active neutrino and neutral Higgs. As we will demonstrate, the complete lepton asymmetry is increased by a factor of 1.98 as compared to the contribution from the charged channels (9) alone (the case considered in the estimate of [43]). It will turn out that the estimate of [43]) for the lepton asymmetry is of the same order of magnitude as the one derived in our current accurate treatment, thus providing an a posteriori justification of the simplified analysis of [43].
In [43], by assuming the heavy Majorana neutrino at rest, we estimated the lepton asymmetry induced by the (Lorentzand-CPT-violating) background B 0 . We assumed one single 4 Scattering processes l l →hh or l h →lh, are of higher order in the Yukawa coupling y and hence are suppressed in our case, although such processes are equally important in standard CPT invariant, CP violating leptogenesis, with more than one species of right-handed neutrinos, as they are of the same order as the CP violating one-loop graphs [24].
Majorana neutrino N with the corresponding Yukawa coupling for the Higgs portal y. For N , the tree-level decays (cf. Fig. 1) for the two channels (9), in the presence of the background B 0 , yields in that case: The decay process goes out of equilibrium when the total decay rate drops below the expansion rate of the universe. Assuming standard cosmology [43] during the decoupling period, 5 which is also hypothesised to coincide with the radiation-dominated era of the Universe, this expansion rate is given by the Hubble constant [63] where N is the effective number of degrees of freedom of all elementary particles and M pl is the Planck mass. For a minimal extension of the SM, with only right-handed neutrinos and the background B 0 , we may estimate N = O(100) at temperatures higher then the electroweak transition [64]. From the last equation one can estimate the right-handedneutrino decoupling temperature T D , in terms of the phe- 5 Such an assumption is non trivial and depends on the microscopic model considered. For instance, in terms of brane-world scenarios for the background B 0 [43], where the latter is derived from a cosmological Kalb-Ramond axion field b(t), such an assumption is justified by requiring a cancellation of the constant in time kinetic energy density of the field b by the (negative) dark energy of the higher-dimensional bulk. After the decoupling, where the string/brane Universe undergoes a phase transition, the dark energy falls off with the temperature sufficiently rapidly, so as today it reaches the value measured by cosmological observations. We shall not discuss such details in the current article.
nomenological parameters , |y| and B 0 [43] T D 6.2 · 10 −2 |y| Imposing a delayed decay mechanism, as for the standard leptogenesis [20,63,65], leads to the further requirement that T D ≤ leading to: In [43] we demanded that saturation of this inequality be satisfied for all values of the background field B 0 , which implies On assuming for the (phenomenological) coupling y the value |y| ≈ 10 −5 , we then obtain an order of magnitude estimate m N for the heavy neutrino mass In [43] we estimated the lepton number density by assuming that all the right-handed neutrinos were at rest before the decay; hence with branching ratios of the decays given by r = 1 and 1 − r , the decay of a single neutrino produces the lepton number Multiplying this quantity by the initial abundance of righthanded Majorana neutrinos N D at the temperature T D (averaged over the respective helicities), one gets a cr ude estimate of the lepton number density. Also, in [43] we assumed that the right-handed neutrino density distribution follows closely the equilibrium distribution for T ≥ T D and drops rapidly to zero at lower temperatures T ≤ T D ; furthermore the density of the sterile neutrino (normalised to the entropy density) is well approximated by a step-function. This implies that the total lepton asymmetry (normalised over the entropy density) produced in the full decay of the right-handed neutrino is given by [43] where is the total entropy density (assuming, for temperatures higher than the electroweak phase transition, SM-like values for the effective degrees of freedom N ∼ 100). For the non-relativistic right-handed neutrino, the Fermi-Dirac equilibrium densityn eq N is well approximated by the Maxwell distribution, yielding in the presence of the background B 0 : where g N = 2 is the effective number of degrees of freedom of the right-handed neutrino, and we assume that B 0 /m N 1, an assumption that proves to be self consistent. The lepton asymmetry L T OT s has not been measured directly, hence it can -depending on the theory -be different from the baryon asymmetry. However in theories with sphaleron transitions that preserve Baryon-minus-Lepton (B−L) number, such as minimal extensions of the SM with right-handed neutrinos, as the ones we are interested in [43] and here, L T OT /s is expected to be of the same order of magnitude as the baryon asymmetry (20), where n B (nB) is the number density of baryons (antibaryons) in the universe, provided it is communicated to the baryon sector by Baryon and Lepton number violating but Baryonminus-Lepton (B−L) conserving sphaleron processes in the SM sector. An order of magnitude estimate of the ratio B 0 m can be found making use of the approximation T D m N and retaining only first order terms in B 0 m 1. Equating the expression for the lepton asymmetry with the phenomenological value (20), and expanding (17) to first order in B 0 /m N , we obtain (for g N = 2) which implies [43] The small value of this ratio also allows us to justify a posteriori the neglect of higher powers of B 0 in the formulae above.
For the case where y = O(10 −5 ) and from the lower bound for m N of 100 TeV found in (14), we get an approximation for the smallest possible magnitude of the background field required in order for this mechanism to be effective: B 0 1 MeV. If other mechanisms contributed to the lepton asymmetry in the universe, or the Yukawa couplings assume smaller values, the minimum value of B 0 would be smaller than the one given here. Baryogenesis is then assumed to proceed via B-L conserving processes in the SM sector of the model. In order to get a physically correct and more accurate estimate of the induced lepton asymmetry, the relevant Boltzmann equation needs to be studied in detail, since the heavy right-handed neutrinos are not at rest, but characterised by the Maxwell-Boltzmann momentum distribution in the early universe. This requires a good approximation for the thermally averaged decay rates (9) of all the relevant processes and will be the subject of the current article. As the Boltzmann equations associated with the leptogenesis scenario advocated here and in [43] involve appropriately averaged thermal rates of the decays (9), we develop in Appendix 7.D the relevant formalism (for B 0 /m N 1); the formalism will be used in the next Sect. 3 to set up the pertinent system of Boltzmann equations. We shall often borrow methods and techniques from the standard case of CPT conserving RHN-induced leptogenesis, where the CPTV background B 0 is absent, but there is CP violation in the lepton sector [20,21,24]. In the current article we shall closely follow the formalism outlined in [24].

Setting up the Boltzmann equations for leptogenesis in the presence of CPTV backgrounds
In the presence of the weak background B 0 the following Boltzmann equation for the number density n r of a fermion species χ of mass m χ and helicity λ r , has been derived in the Appendix of [43]: where g denotes the number of degrees of freedom, and f is the Fermi Dirac distribution of a relativistic fermion assuming zero chemical potential: The B 0 dependent energy-momentum dispersion relation (cf. Appendix 7.B) should be used and an expansion up to and including first order terms in the background B 0 /(m N , T ) is performed for our weakly CPTV background. The term C[ f ] denotes the appropriate thermally averaged decay or interaction rates involving the species χ [64]. In practice, it is convenient when calculating the lepton asymmetry, to consider the number densities normalised over the entropy density of the universe (18) [64]: In the problem at hand, we consider a system of Boltzmann equations, associated with the heavy neutrino N , as well as the lepton l ± abundances. The Boltzmann equation (23) applies to both a relativistic (massless) neutrino as well as a heavy right-handed neutrino, upon using the appropriate dispersion relation (25). We shall follow the standard analysis in constructing the relevant equations [24], with the important difference being that the energy momentum dispersion relation and the interaction rates C[ f ] term involve now the LV and CPTV background B 0 .
In terms of the abundances (26), the Boltzmann equations associated with the interactions (9) of a RHN with a given helicity λ take the form: where Y N is the heavy neutrino abundance, and the superscript eq denotes thermal equilibrium quantities. The equilibrium abundances Y eq N are discussed in detail in Appendix 7.C; the γ eq,(λ) (N ← → ± h ∓ ) denote the appropriate thermally averaged decay rates, discussed in Appendices VIIB and VII D. We shall use their explicit expressions later on, in order to construct the final form of the Boltzmann equations. The term λ I in (27) is a generic notation for an appropriate integral stemming from the terms proportional to the CPTV background B 0 and the helicity λ on the left-hand-side of (23). Such terms vanish when we average over helicities, since r λ r = 0. The reader should notice that apart from the λ I term, the rest of the structures in (27) are the same as in conventional CPT invariant but CP violating cases for leptogenesis [24]; but, as already mentioned, the relevant dispersion relations (25) are modified by the CPTV background From the expressions for the relevant amplitudes in Appendix 7.B, we know that, on account of helicity conservation, for the processes N ← → l − h + we only have one helicity λ = −1 and for the processes N ← → l + h − we only have λ = +1. Following standard treatments [24], we also take the charged Higgs boson as well as the charged leptons to be roughly in equilibrium; hence we set Y l,h Y eq i,h for the corresponding abundances in (27), and find: z Hs dY Next we will generate the lepton and anti-lepton Boltzmann equations, which are needed in the calculation of the lepton asymmetry. As there is only one forward and reverse process for a lepton l − with a definite helicity λ = −1, the corresponding Boltzmann equation obtained from (23), reads z Hs dY Again we take the Higgs particle to be in equilibrium Y h + Y eq h + [24]. Moreover, from the relevant discussion in Appendix 7.B, we know that we only have one helicity (λ = −1) for the processes concerning the leptons l − , which implies that the Boltzmann equation for the lepton becomes z Hs dY Applying a similar analysis, but now concentrating on the opposite helicity λ = + 1, we arrive at the Boltzmann equation for the anti-lepton l + : z Hs dY In the specific leptogenesis scenario of [43], the leading contributions to the lepton asymmetry (as far as the small Yukawa coupling (5), y ∼ 10 −5 1, is concerned) come from the tree level decays (9) and their reverse processes. As already mentioned in the previous section, the additional interactions lh →lh and ll → hh, involving a tree-level heavy neutrino exchange, are both of higher order in y and suppressed by the heavy mass m N , hence they will be ignored in our case. (It should be remarked that these latter interactions yield contributions comparable to the one loop order graph of Fig. 1 and hence play an important rôle in CPT invariant, conventional leptogenesis scenarios [24]).
From now on, we shall concentrate on constructing the system of Boltzmann equations associated with: (i) the heavy neutrino abundance in units of entropy density (cf. (26)), and averaged over helicities λ = ± 1: and (ii) the lepton-asymmetry for the processes (9), defined in terms of the lepton abundances: where we took into account that the asymmetry is generated between the leptons of helicity λ = −1 and the anti-leptons of helicity λ = +1, since these are the only decays for the heavy neutrino (9), for each of which helicity is conserved. There will be no asymmetry between leptons of helicity λ = +1 and anti-leptons of helicity λ = −1 and so Y Moreover, all of the negative helicity lepton abundance Y l − comes from the decay of the negative helicity heavy neutrino. The same argument for the anti-lepton positive helicity abundance generated by the positive helicity heavy neutrinos. These imply the second of the relations (33).
The total observable lepton asymmetry, which we want to compute, and compare the result with the estimate (21), is defined with respect to the corresponding abundances (averaged over helicities) in units of the entropy s, as follows: on account of (33). In what follows we proceed with the construction and solution of the Boltzmann equations that correspond to the quantitiesȲ N and L.
To obtain a Boltzmann equation, summed up over helicities, for the averaged RHN abundanceȲ N (33) from the system (28), we sum up these equations, to obtain: 2z Hs dȲ N dz The asymmetry (34) will be evaluated at decoupling temperatures by solving explicitly the appropriate system of Boltzmann equations for L andȲ N and the result will be compared with the estimate (21) of [43]. In solving the equations we shall approach decoupling by starting from high temperatures T and gradually approaching decoupling T → T D by making use of appropriate approximations (Padé approximants [54][55][56][57][58][59][60]), which will allow for analytic expressions for the lepton asymmetry.
In this high-temperature (relativistic) regime, the entropy density of the Universe scales with T as s ∼ 14T 3 , whilst the Hubble parameter behaves as [64], H ∼ 6T 2 /M pl , with M pl the Planck mass. Using these relations, we can write The terms λI that appear on the left hand side of the Boltzmann Eqs. (28), (30), (31), in the high-temperature regime T m χ for a generic fermion of mass m χ , and degrees of freedom g χ , can be written as: We only have to consider the (massless) lepton case and expand the series upto second order, The integral I l can therefore be expressed as, where the integration variable was changed to |p l |/T = x.
The expression for I l up to second order is given by, is the relativistic energy of the lepton and is taken to be independent of B 0 , since in our analysis we are only considering terms of linear order in B 0 T, m N [43]. All series expansions are taken to second order in the appropriate small parameters, for reasons that will become clear below, when we consider the Padé approximated analytic solution for the Boltzmann equations extrapolated to the RHN decoupling temperature T D m N (13), (15).
The integral I χ , in the lepton case, can be approximated by Hence, from (35), (30), (31), (36) and (41), we observe that the Boltzmann equations for the heavy neutrino abundance and lepton/anti-lepton asymmetry L, averaged over helicities, in the high temperature regime, acquire the form (we reminder the reader that the leptons l ± are strictly massless, m l ± = 0, in the high temperature regime, above the electroweak phase transition): and 84 with the definitions We next proceed to solve the above equations which, since they are linear and first-order, can be in principle exactly solved. However, for the exact solutions to be amenable to analysis, approximations will need to be made; the goal is to find an analytic expression for the lepton asymmetry.

Heavy-right-handed-neutrino abundance Boltzmann equation
We commence our analysis with the heavy-RHN-Boltzmann equation (42). The equilibrium populations are calculated in Appendix 7.C. The corresponding thermally averaged decay rates read (see Appendices VIIB and VII D, and in particular Eq. (188)): The reader should notice the "reciprocity" equalities even in the presence of the CPTV background B 0 = 0. These are consequences of the equality of the corresponding amplitudes (140) and energy conservation, as explained in Appendix 7.D. Also, it is immediately seen from (45) that it is only in the presence of the CPTV background B 0 = 0 that a lepton asymmetry is generated at tree level between the decay channels (9) (see Fig. 1), as a consequence of the pertinent differences in (45) and (140). In this respect, the similarity of the rôle of the CPTV ε 1 parameter with the corresponding one, ε, of conventional leptogenesis [24] should be noticed.The important difference is that, in contrast to our CPTV case, conventional lepton asymmetry occurs at one loop level for the decays of Fig. 1 and requires more than one flavour of the RHN. After substitution of the relevant expression for the thermally-averaged quantities γ eq , we have the following intermediate results (for details see Appendix 7.C), from which it follows where to obtain the last expression of (50) we have expanded the function in the round brackets in the definition of Y (λ),eq N up to second order in z < 1, neglecting terms of order The remaining terms in the Boltzmann equation (42) become: We now evaluate the sum and difference of the abundances normalised to their respective equilibrium values, Substituting these expressions in (51), we obtain where again the term involving the differences of the abundances will be of order B 2 0 since ε 1 (z) is already linear in B 0 and so is neglected. We may write the right-hand-side of the heavy neutrino Boltzmann equation (42) as: Upon substitution of the relevant expressions, the heavy neutrino Boltzmann equation at high temperatures becomes: which can be finally written as: We stress once more that this equation is derived in the high temperature regime in which m N < T .

Lepton asymmetry Boltzmann equation
We proceed now to study the equation for the lepton asymmetry (43) at high temperatures. Concentrating on the first two terms on the right hand side, which involve the heavy neutrino abundances, and substituting in the expressions for the thermally-averaged γ eq integrals (cf. Appendix 7.D), we obtain after some straightforward manipulations: where we have substituted in the expressions for the sum and difference of the heavy neutrino abundances from the previous section. The final two terms on the right hand side of the lepton asymmetry Boltzmann equation (43) can be expressed as: We next evaluate the sum and difference of the lepton and anti-lepton abundances normalised to their respective equilibrium values, that is, the quantities Using the explicit expressions for the equilibrium abundances for leptons and anti-leptons (cf. Appendix 7.C), we obtain Then (58) yields where the reader should recall that ε 1 (z) is already linear in The final form for the lepton-asymmetry Boltzmann equation at high temperatures, then follows: As with the equation for the RHN abundance, the reader should bear in mind that the lepton asymmetry equation above is derived in the high temperature regime m N < T .

Solutions to the system of Boltzmann equations
In this section we derive approximate analytic solutions of the system of Boltzmann equations (56), (62), which will allow us to compute the lepton asymmetry induced by the CPTV background in our model. So far we have derived equations for the RHN and lepton asymmetry (cf. (56) and (62) respectively) for high temperatures, z < 1. However, we are eventually interested in solutions of the corresponding Boltzmann equations at the RHN decoupling temperatures (13), (15), where z ∼ 1 [43]. We shall attempt to extrapolate our results above to this case, by performing a Taylor expansion of the series solutions to these differential equations. The expansion takes place around an arbitrarily chosen point in the interval 0 < z < 1, where the solution is valid, taking proper account of the (thermodynamic equilibrium) boundary conditions for the abundances as z → 0 (see Appendix 7.C), which fixes the integration constants characterising the solutions. In our analysis below, we take, as a Taylor expansion point, the mid-point of the interval (0, 1), z = 0.5 .
To extrapolate the solutions to the regime z 1, we shall use a Padé approximation [54]. As well known, a Padé expansion can accelerate the convergence of an asymptotic expan-sion or, for a series, turn a divergence into a convergence. It is widely used for producing in solving approximately complicated problems in several fields of physics, ranging from statistical mechanics to particle physics and quantum field theory [55][56][57][58][59][60]. Here we present another useful application of the method in cosmology. We outline the general concepts of the Padé approximants method and the specific algorithm used in our computation in this work in Appendix 8.

Solution to the heavy-neutrino Boltzmann equation
The heavy neutrino Boltzmann equation (56) decouplesȲ N from L so the former can be obtained by solving this equation with an appropriate integrating factor [66,67]. We therefore commence our discussion with a sketch of the solution of Eq. (56). Calling the equation reads The integrating factor for this differential equation is given by, Multiplying through the differential equation by the integrating factor gives where c 1 is the constant of integration and will be determined using the boundary condition (cf. (159) in Appendix 7.C), where for heavy right-handed neutrinos g N = 2, and we used (63). In our case 0 < z < 1, as we are interested in non-trivial populations in the phase where T > T D (for T < T D the populations drop sharply, this is our basic assumption [43]). From the qualitative analysis of [43], reviewed in Sect. 2, the freezeout temperature T D is expected to be of order (cf. (13), (15)): T D m N so z D 1. This is why it is important to give formal solutions first, before any expansion. Notice that in arriving at the system of Boltzmann equations forȲ N and L, we did not make more assumptions on the magnitude of z other than it belongs to the interval 0 < z < 1.
We now make some approximations in order to obtain a solution for the heavy neutrino abundance. We can write the integrating factor I N (x) as in order to simplify this expression we only take the first two terms in the series S n S 0 + S 1 .
where we have expanded again to first order the (upper) incomplete Gamma functions [68] that arise in this integration, The boundary condition (67) determines the value of the constant of integration: c 1 = 399.1256b 2 = 2.2351. After taking the inverse of the integrating factor (keeping first order terms), we obtain the expression for the abundance of the heavy neutrino in the interval 0 < z < 1, where any exponential factors that remain after multiplying by the inverse of the integrating factor have been expanded to first order. Also any terms of higher order than z 32/3 have been neglected from the expression due to the restriction 0 < z < 1 and any terms with factors of order 10 −5 or smaller have also been neglected.

Solution to the lepton asymmetry Boltzmann equation
In this subsection, we proceed with the substitution of the previous result onto the Boltzman equation (62) and proceed with its solution, which will allow for a determination of the lepton asymmetry.
Similarly to the previous case, the integrating factor I L for the lepton asymmetry Boltzmann equation is given by with the lepton asymmetry itself, being expressed as where c 2 is the constant of integration, determined by using the thermal equilibrium boundary condition (c.f. Appendix 7.C, Eq. (159)), After substituting in the solution for theȲ N (z) in the interval 0 < z < 1 the formal lepton asymmetry solution is given by, As in the previous case we make some simplifying approximations to obtain a solution for the lepton asymmetry. The integrating factor is approximated by the expansion of the series up to first order, Now that an approximate solution forȲ N (z) is known we may express the coefficient H (x) as, where we have neglected terms of higher powers then x 32/3 . We then have to solve the integral below, To determine the constant of integration c 2 we use the boundary condition L(z → 0) → L eq (z → 0) = 0 which yields c 2 = −0.8177 B 0 m N . Now multiplying by the inverse of the integrating factor (to first order) we obtain an expression for the lepton asymmetry in the interval 0 < z < 1.
similarly we have neglected terms of higher order powers than z 32/3 and any terms with factors of order 10 −5 or smaller. Now we want to estimate the lepton asymmetry at freeze out where T D ≤ m N corresponding to z ≥ 1.
To this end we Pade expand [54] (cf. Appendix 8) the expressions for L(z < 1) andȲ N (z < 1) around the point z = 0.5 in order to make the expressions for the abundances valid beyond the interval 0 < z < 1. We require a positive asymmetry L, as this is the only physically relevant solution for dominance of matter over antimatter, for our fixed sign of the background B 0 > 0. From (81) we observe that L(z < 1.44) < 0, hence we must have z = z = 1.44 as a critical value in our approximate treatment below which the lepton asymmetry switches sign. We interpret this as determining the freezeout point, after which (T < T D ) the asymmetry freezes out to a positive value. For this value we have and thus the observable lepton asymmetry (34) is given by, The reader should compare this result with that obtained in [43], see Eq. (21) above. Our result (84) yields a lepton asymmetry proportional to B 0 /m N as in (21), but with a proportionality coefficient which is 1.94 times larger. The fact that it is larger may be attributed physically to the fact that here we considered the non zero momentum modes of the heavy neutrino in estimating the asymmetry, which were neglected in [43]. Nevertheless, we consider this a good agreement between the two results. We have shown above that this lepton asymmetry can be generated at the freeze out point z = 1.44 (in order for a positive asymmetry) using first order approximations to the formal solutions of the abundances, this still satisfies the condition that freeze out should occur at T D ≤ m N . It is important to notice that the order of magnitude estimate for the Yukawa coupling |y| ∼ 10 −5 in earlier work [43], which was used throughout our previous calculations, providing numerical input (eg. (63)) into the approximate solutions, remains unchanged, and this provides a posteriori a self-consistency check of our approximation. The decoupling (82) now occurs at 1.44 T D = m N instead of the assumed one in [43] at T D m N , but this does not alter the order of magnitude of the Yukawa coupling. However, we believe that the fact that the asymmetry turns negative for z < 1.44 is an artefact of the approximations used. Full numerical analysis may lead to a freezeout point z D 1 as in [43]. To check on the stability of the freezeout value, we present next an alternative approximate derivation.

Series solutions of the Boltzmann equations
Here we present another method of obtaining the (approximate) solutions to the differential equations, in an attempt to get an idea on the stability of the freezeout point. Starting with the heavy neutrino Boltzmann equation we can Taylor expand the variable coefficients P(z), Q(z) around the point z = 0.5 and the solutionȲ N (z), On substituting these series into the differential equation we obtain We can then see a recurrence relation for the coefficients of the solution forȲ N (z) in terms of the coefficients of the P(z) and Q(z) series, Using this recurrence relation, the first few coefficients are: The Taylor expansion around the point z = 0.5 of the heavy neutrino abundance is then, We now take the limit z → 0 in such a way that the boundary condition (67) is satisfied, that is,Ȳ N (z → 0) →Ȳ eq N (z → 0) = 0.0335. This places the final constraint in order to obtain the value for the last remaining coefficient c 0 = 0.0339. The final expression for the heavy neutrino abundance around z = 0.5 is given by, We proceed with the analogous calculation for the lepton asymmetry Boltzmann equation, The recurrence relation is similar to (87) under the change p n → j n , q n → h n , c n → l n where l n are the coefficients in the lepton asymmetry Taylor expansion, The coefficients for J (z) and H (z) are: The coefficients l n are given below using the recurrence relation, We use the boundary condition (cf. (159) in Appendix 7.C) L(z → 0) → L eq (z → 0) = 0 to find the last coefficient l 0 = −0.0189 B 0 m N and the final expression for the lepton asymmetry is given by, We now perform a Padé expansion [54] (cf. Appendix 8) around the point z = 0.5 to be able to use the solutions outside the interval 0 < z < 1. In order to obtain a positive asymmetry, we observe from (94) that we must have z ≥ 1.62, thus in this approximation the critical point appears to be at z * = 1.62. This is identified with the freezeout, which, upon substitution into the Padé approximant for the lepton asymmetry, yields L(z = 1.62) = 0.0005 B 0 m N , with the corresponding heavy neutrino abundance at this point is Y N (z = 1.62) = 0.0307. The observable lepton asymmetry (34) in that case is found to be We see that the series solutions yield a similar answer to the method using an integrating factor. The point of decoupling z D = 1.62 still satisfies T D ≤ m N ⇒ z ≥ 1 and the order of magnitude estimate for the Yukawa coupling |y| ∼ 10 −5 is unchanged. Comparing with (21), we see that the result (96) is in excellent agreement with the lepton asymmetry estimated in [43].
From either (82) or (96), we obtain that phenomenologically relevant leptogenesis in our system, in the sense of (21), is achieved for B 0 /m N = O(10 −9 − 10 −8 ), which is in the same approximate range as the estimate of [43], but here the result includes all the non-zero momentum modes of the heavy neutrino. This implies that for m N = O(100) TeV, we must have a B 0 in the range B 0 ∼ 0.1 − 1 MeV for leptogenesis to lead to the observed baryogenesis via the B-L conserving sphaleron processes.
Comparing the freezeout points between the two approximate methods (82) and (95), we observe agreement with only 12.5 % uncertainty, indicating stability of the freezeout point in the region around one. This completes our analysis. Perhaps as we mentioned earlier, a full numerical solution will yield a freezeout point closer to the qualitative value of [43], although we should emphasize that the above approximate analyses have yielded results in this respect that are of the same order of magnitude. This adds confidence to the efficient application of Padé approximant method to our cosmological problem.

Inclusion of the neutral Higgs portal
In the system of Boltzmann equations we will now include the contributions from the decays of the RHN into a neutral Higgs field h 0 and an (active) neutrino ν of the SM sector (see Eq. (10)). As we shall demonstrate below, the computed (complete) lepton asymmetry has an approximate factor of two compared with the one based only on the charged lepton decay channels (9) .

The complete heavy neutrino Boltzmann equation
Upon considering the additional contributions of the decay channels N → νh 0 and N →νh 0 to the Boltzmann equation for the abundance of a RHN with helicity λ, Y We start by considering only the negative helicity λ = −1.
The reader should recall, from our analysis in the previous sections, that only the decay channel N → l − h + and its inverse process yield non trivial contributions to the Boltzmann equation. Adding the decay channel N → νh 0 (10) and its inverse, for helicity λ = −1, to the Boltzmann equation, yields, z Hs dY We now make the approximation that the Higgs and lepton fields are in equilibrium Y Y eq . We also have the reciprocity of the thermal decay rates into charged leptons γ eq,(−) (N → l − h + ) = γ eq,(−) (l − h + → N ) and the same will be true for the light neutrino decay channels γ eq,(−) (N → νh 0 ) = γ eq,(−) (νh 0 → N ). Thus, the Boltzmann equation for the abundance of the negative helicity RHN becomes, z Hs dY Similarly for the decay into anti-particles we only have the positive helicity λ = +1 and the resulting Boltzmann equation reads, z Hs dY The thermal equilibrium decay rates of RHN into charged leptons and neutral leptons will be identical, since the magnitude of the electric charge was not specified, when considering the charged leptons. Only the relativistic nature of the corresponding dispersion relations was important. Therefore we have γ eq,(−) (N ↔ l − h + ) = γ eq,(−) (N ↔ νh 0 ) and γ eq,(+) (N ↔ l + h − ) = γ eq,(+) (N ↔νh 0 ). Taking this into account, and combining the two Boltzmann equations, expressed in terms of the averaged-over-helicities abundancesȲ N , we obtain, The heavy neutrino Boltzmann equation can then be expressed as, This is the complete Boltzmann equation for the abundance of the heavy right-handed Majorana neutrino, averaged over helicities, obtained by considering the complete set of decay channels of the heavy neutrino into charged and neutral leptons and anti-leptons (9), (10). The reader should notice that the constant factors a 2 and b 2 appearing in (101) are twice as large as compared with those in the case where only the RHN decays into charged leptons and anti-leptons were considered, cf. Eqs. (63) and (64).

The complete lepton asymmetry Boltzmann equation
The inclusion of the extra decay channels (10), leads to two additional Boltzmann equations when calculating the lepton asymmetry, Upon using the approximation that the Higgs is in equilibrium and the fact that for decays into particles we only have helicity λ = −1, whereas for decays into anti-particles we have λ = +1, as well as the reciprocity of the thermal decay rates γ eq,(λ) (N → A, B) = γ eq,(λ) (A, B → N ), the above equations become, The complete lepton asymmetry is given by the difference between the lepton abundances and the anti-lepton abundances, Y ν . By combining the last four equations we arrive at the following result, Using the same arguments as in the case of the heavy neutrino Boltzmann equation above, we realise that the abundances of the SM (active) neutrinos are identical to that of the charged leptons Y l − = Y ν , Y l + = Yν. The thermally averaged decay rates will also be equal γ eq,(−) (N ↔ νh 0 The Boltzmann equation for the complete lepton asymmetry (L = Y This is identical in form to the previous asymmetry Boltzmann equation (62) (for the case where only the decays of the heavy neutrino into charged leptons and anti-leptons were considered). However, in view of the difference of the corresponding Boltzmann equation forȲ N , we expect some differences in the solutions for L, which we now proceed to evaluate.

Integrating factor solutions of the complete Boltzmann equations
We consider here the solutions to the complete Boltzmann equations on including the additional decay channels into SM neutrinos and neutral Higgs. As seen in (101), the only difference is an increase of the constants a 2 and b 2 by a factor of two when considering the additional decay channels. The solution obtained using the integrating-factor method is completely analogous to that shown is Sect. 4.1; the complete RHN solution is found to bē Similar to the previous solution, terms of higher order than O(z 32/3 ) have been neglected as well as terms with coefficients O(10 −5 ). The complete RHN solution is then substituted into the right hand side of the complete lepton asymmetry equation (105) (which has the same form as the lepton asymmetry Eq. (62)). Following the calculation given in (4.2), the solution for the complete lepton asymmetry is given by, (z D ) = 0.0343 (3.3% increase from the previous solution). The complete observable lepton asymmetry is given below: On using the integrating factor method, we note that the observable lepton asymmetry has increased by a factor of 1.4 from the previous calculation (where RHN was considered to decay into charged leptons only); so the addition of the neutral decay channels increases the complete observable lepton asymmetry.

Series solutions of the complete Boltzmann equations
For concreteness we consider the series solution of the heavy neutrino abundance around the point z = 0.5. As we have seen above, (101), when the RHN decay mode into SM neutrinos is taken into account, the only difference from the previous case, where only RHN decays to charged leptons have been considered, is that the factors a 2 and b 2 in (101) are twice as large as compared with their counterparts in (64). This has the effect that the coefficients p ı and q j in the series solutions also increase by a factor of two, the rest of the calculation being identical. The Taylor expansion of the RHN abundance around z = 0.5 is, This solution must be inserted into the lepton asymmetry Boltzmann equation;Ȳ N only appears in the H (z) coefficient of (62) and so the coefficients h ı are the only quantities to change; these are given below: The coefficients l ı are calculated in the same way as before and the Taylor expansion for the complete lepton asymmetry around the point z = 0.5 is We now perform a Padé expansion of this expression around z = 0.5, we find that to generate a positive asymmetry we must have a later decoupling of z D = 1.77 where L pade (z D ) = 0.001B 0 /m N . (100% increase from the previous solution.) The RHN abundance at decoupling is Y pade N (z D ) = 0.0313 (less than a 2% increase from the previous solution). The complete observable lepton asymmetry is given by We thus note that the observable lepton asymmetry has increased by a factor of 1.98 from the previous calculation when only considering the RHN decay into charged leptons only, and so the addition of the neutral decay channels effectively doubles the lepton asymmetry. The fact that decoupling now occurs at z D = 1.77 does not affect the order of magnitude estimate for the Yukawa coupling y.

Conclusions and outlook
In this work we have completed the analysis presented in an earlier work [43] by computing the lepton asymmetry generated due to the decays of heavy right-handed neutrinos in the presence of a CPTV axial vector background with only temporal components B 0 = 0 in the early universe through an analytic (but approximate) solution of the corresponding algebraic system of Boltzmann equations. In [43] we only presented a heuristic estimate of the generated asymmetry. In order to facilitate the comparison of our detailed analysis with the heuristic order of magnitude estimates of [43] we have first assumed the active neutrino to be purely Majorana as in [43] and concentrated only on the decays of RHN into charged leptons. The inclusion of the decays to neutral Higgs and light neutrinos has also been done here, with the (expected) result that the total asymmetry is increased by a factor of about 2 as compared to the charged-lepton-case.
Let us now recapitulate the main result for the case where the neutral decay channel is ignored. The current solution of the Boltzmann equations that describe the leptogenesis in the model has been obtained through an appropriate Padé approximation around the point z = m N /T = 0.5. This procedure allowed the power series representation of the lepton asymmetry to be extrapolated outside the interval 0 < z < 1 and to be evaluated at the point z = 1.44 to show the positive lepton asymmetry. The obtained result for the asymmetry is in qualitative agreement with the estimate of [43], in that it is proportional to the small quantity B 0 /m N 1. However the proportionality coefficient in the case the solutions are evaluated using an integrating factor is found to be 1.94 times larger than in the case of [43]. On the other hand, in case one uses a series solution to the Boltzmann equations, the proportionality coefficient is in excellent agreement with the case of [43]. This implies that in our numerical treatment the lepton asymmetry can be estimated to be in agreement with the estimate (22) of [43]. In our analysis we assumed Yukawa couplings of order |y| ∼ 10 −5 in the Higgs portal term (6), that couple the right-handed neutrino to the SM sector of the model. This prompted us to ignore higher order terms of order |y| 4 ∼ 10 −20 B 0 /m N , which a posteriori was proved to be a self-consistent result, due to the smallness of the B 0 /m N (114), required for the observed baryon asymmetry today (through leptogenesis).
The inclusion of the neutral channel in the RHN decay yields an estimate of the lepton asymmetry which is (2.3 − 2.7) times as large as the one in [43] and between (1.4 − 2) times as large as the asymmetry calculated in (113) implying a ratio B 0 /m N which is (less than) an order of magnitude smaller than (22). The slight increase of the freezeout point does not affect the order of magnitude of the asymmetry nor that of the Yukawa coupling of the Higgs portal; hence in order of magnitude there is qualitative agreement with the estimates of [43]. Although our analysis has been generic in not specifying the microscopic origin of the CPTV background, nonetheless some microscopic scenarios originating from string theory have been presented in [43]. According to these scenarios the background is identified with the dual of the Kalb-Ramond antisymmetric tensor field strength, μνρσ H νρσ , which in a four-dimensional space-time is equivalent to the derivative of a pseudoscalar field b(x) (Kalb-Ramond axion), ∂ μ b. Nevertheless such an identification is not binding. However, if it is made, then within the context of realistic brane/string models the pressing question concerns the microscopic mechanism, for the transition from a relatively large value (in the Robertson-Walker frame) of the constant B 0 = 0 CPTV background in early eras of the (string) universe, required for leptogenesis, to a very weak background today, compatible with the very stringent limits of CPTV in the current era [38][39][40][41]. Some speculations have been presented in [43] but detailed microscopic mechanisms, compatible with the rest of astroparticle phenomenology, including the open issue of the smallness of the (observed) cosmological constant (or dark energy) today, are still lacking and will be the subject of future investigations.
Nevertheless, we believe that the scenario for baryogenesis through leptogenesis presented initially in [43] and completed here, is an attractive, relatively simple one, which deserves further investigations, within the context of appro-priate microscopic models (not necessarily within the framework of string/brane theory). We hope to come back to such studies in the near future. Another important aspect of our current work is the demonstration of the efficiency of the Padé approximant method [54] in solving Boltzmann equations, thus adding yet another successful example of this method, this time of relevance to cosmology.
Before closing, for the case that the background B 0 originates from microscopic string-inspired models [43,69], we will make some comments on the coupling of the torsion-like antisymmetric tensor Kalb-Ramond field strength to fermion species ψ i in our model. Due to its gravitational origin, this coupling concerns the axial fermion currents of all fermions, including quarks in the SM sector where b(x) is the Kalb-Ramond axion field, dual to H μνρ in four space-time dimensions, as discussed briefly in Sect.
2. In such a case, direct CPT induced baryogenesis occurs, given that the CPTV background B 0 =ḃ = const would imply effective chemical potentials that are different (by a sign) between particles and antiparticles (and left-and righthanded chiral fields). This, for the case of quarks of the SM sector, can lead to direct baryogenesis. However, for the scenario discussed in the present paper, the magnitude of the background B 0 =ḃ (where the overdot denotes time derivative) for temperatures at or below the electroweak phase transition T ≤ 10 2 GeV -pertinent for such a case, assuming more or less standard decoupling temperatures of quarkswould imply a contribution to the baryon asymmetry that would be much smaller compared to that induced by the lepton asymmetry for B 0 ∼ 1 MeV, at T D = 10 5 GeV. This is due to the scaling of B 0 with the cubic power of temperature T , as discussed in [43], and reviewed above (cf. (8)).
In general, though, several other contributions to the baryon asymmetry are expected, which depend on the microscopic model considered. Indeed, as discussed in [69], where various microscopic scenarios for CPTV induced matter-antimatter asymmetry in the universe have been studied, the induced asymmetry is of order B 0 (T D )/T D , where T D is a decoupling temperature depending on the model. For T D less than the electroweak phase transition, as mentioned above, such contributions are much smaller than the lepton asymmetry calculated here, due to the cooling law (8), and hence our mechanism of baryogenesis through leptogenesis advocated here would be the dominant one. However, this analysis holds only for the particular model of the CPTV background originating from Kalb-Ramond fields. In general for models with CPTV originating from gravitational space-time defects [69], the situation could be different and the study of such models constitute interesting avenues for future research.

Appendices
In the following Appendices we discuss in detail several technical aspects of our work, which have been used in various parts of the main text.

7.1: Notation and conventions
Throughout this work we use the following conventions. Our metric signature convention is: which implies The Dirac γ matrices have the properties (we use the symbol ı to denote the imaginary unit) The chiral representation for the Dirac matrices will be used throughout: with the 2 × 2 Pauli matrices In this Appendix we work out the amplitudes for the decay channels (9) in an arbitrary frame, where the decaying right handed neutrino N has a four-momentum p μ , μ = 0, . . . 3. This generalises the approximate treatment of [43], where the field N was assumed at rest. Our starting point is the Lagrangian for (Dirac) spinors in an axial Background B μ , which is taken to be purely along the temporal axis (B μ → B 0 ), with B 0 a small, positive (by convention), non zero constant, 0 < B 0 1: The corresponding (Dirac) equation of motion reads On assuming plane-wave solutions for the spinor ψ, corresponding to positive (ψ(x) = u( p)e −ı px ) or negative (ψ(x) = v( p)e +ı px ) frequencies, separately, and substituting in (120) we easily obtain [43] the pertinent polarization spinors u( p) (v( p)) for the positive-(negative) frequency solutions, of helicity λ r = ±, r = 1, 2, in the presence of the background B 0 are given by [43] with u (v) pertaining to the (anti) particle, respectively; ξ r are helicity eigenspinors, satisfying with the helicites λ 1 = −1, λ 2 = +1, and σ i , i = 1, 2, 3 the 2 × 2 Pauli matrices. In the expressions (121) we used the normalisation N ± = √ E r (|p|) ∓ (B 0 + λ r |p)|, with the (−) ((+)) sign referring to u (v) spinors, respectively. The eigenspinors ξ r satisfy the orthogonality condition The energy-momentum dispersion relation for a fermion of mass m in the presence of B 0 = 0 reads [43]: For the Majorana neutrino we have m = m N = 0; on the other hand, the leptons l ± in the early Universe, at temperatures much higher than the electroweak symmetry breaking, of interest here, are massless (m = m l = 0). Thus, the lepton and neutrino energies are explicitly written as: Working out the amplitude for the decay process N → l − h + we obtain where the outgoing lepton spinor isū l − ,s ( p l − ), and the incoming heavy neutrino spinor u N ,r ( p N ); the notation E χ,r (= E λ r χ ) indicates the energy of a spinor χ with helicity λ r . P R = 1 2 (1 + γ 5 ), y is the Yukawa coupling (5) and the orthogonality condition (123) forces the helicities of the incoming and outgoing particles to be the same (helicity conservation). After squaring the amplitude (126) and averaging over initial spins (S = 1/2) we obtain for a given helicity λ, We are now going to consider (127) for the two different helicities λ = ±1. The terms within the first bracket of the above expression take the form where we have substituted in the lepton energy E λ l − = |B 0 + λ|p l − ||. If we take λ = +1, then the above expression is zero (provided B 0 > 0 which is our initial assumption) and so the heavy neutrino of helicity λ = +1 can not decay into leptons.
We now consider the case of λ = −1 for the decay process N → l − h + . The terms in the first bracket of the right-handside of (127) become, We must examine separately the following two cases: (i) when |p l − | ≤ B 0 , the term (130) vanishes, whilst (ii) when |p l − | > B 0 , this term becomes 2(|p l − | − B 0 ). So for the decay process N → l − h + , the only way for the amplitude to be non-zero is when λ = −1 and |p l − | > B 0 . The expression for the amplitude squared for this process is then given by: where we have substituted in the expression for the relativistic heavy neutrino energy for λ = −1, expanded up to second order in small quantities, and neglecting terms of order O(B 2 0 ) (for our purposes, we assume relativistic regime of temperatures, such that 0 < B 0 m N p N ∼ T ): For the reverse process l − h + → N we have, where the outgoing heavy neutrino corresponds to P R u N ,r ( p N ), whilst the incoming lepton to u l − ,s ( p l − ). This process yields the same amplitude as for the decay process N → l − h + , along with the same constraints on the helicity and momentum, For our purposes in this work, we shall extend the range of the lepton momentum to cover all momenta |p l − | ∈ [0, ∞].
For the decay of the heavy neutrino into anti-leptons N → l + h − we have the outgoing anti-lepton spinor v l + ,s ( p l + ) and the incoming heavy neutrino spinorv N ,r ( p N ) with N being its own anti-particle. The amplitude for this decay is Again we square the amplitude and average over the initial spins of the heavy neutrino, to obtain Consider the energy of the anti-lepton for the possible helicities E λ l + = |B 0 + λ|p l + ||. We find that the only two non-zero amplitudes are We will neglect the contribution from the decay amplitude for negative helicity, as it requires |p l + | < B 0 . Then, for the decay process N → l + h − we have where we have substituted in the expansion of E λ=+1 N up to second order. We see that this decay amplitude differs from the previous process under a change of sign of B 0 . The amplitude for reverse process l + h − → N is and we can readily see that it is the same as that of the forward process. The squared amplitudes averaged over initial spins of all the processes are given below: where we see that the forward and reverse processes of each decay yield the same result and the difference between the two decay channels into leptons and anti-leptons is a difference in sign of B 0 .

7.3: Thermal equilibrium populations
The (thermal) equilibrium population of a particle species is given by [64] where f eq is the equilibrium distribution function given by Fermi-Dirac or Boson-Einstein statistics.
We proceed now to determine the equilibrium abundances of the heavy right-handed neutrino (RHN) and the leptons. In the high-temperature era of the universe that we are considering we have T > m N ∼ T D , |p N | > m N and so the particles behave relativistically. The dispersion relation for the heavy neutrino is given by (125) which has been expanded up to second order in small quantities, neglecting terms of O(B 2 0 ). From (141), then, the equilibrium population is given by Where we expand the series to second order, therefore the equilibrium distribution is approximated by.
the equilibrium population then becomes.
each of the above integrals is of the form, where n = 1, 2, 3 and we have expanded out the exponential to second order to record all necessary terms, x t s−1 e −t dt is the upper incomplete Gamma function [68], with the values For our analysis in this work we shall need the averaged over helicities heavy neutrino equilibrium abundanceȲ eq N , and the lepton asymmetry equilibrium abundance L eq , which are given by: For heavy right-handed neutrinos, we have g N = 2.

7.4: Thermally averaged interaction rates
To calculate the thermal equilibrium density integral for each decay process, which enters the pertinent Boltzmann equation, we must sum over the different helicities: where, as discussed previously in Appendix VII B, we will only have the λ = −1 case for the process N → l − h + and the λ = +1 case for the process N → l + h − . The interaction integral for the process N → l − h + is given by: with the equilibrium distribution f eq N = 1/(e E N /T + 1). Above, we have integrated over the momentum delta function, to perform explicitly the integration over d 3p h + , which enforces momentum conservationp h + =p N −p l − . The quantity f (|p l − |) is given below the root of f (|p l − |) = 0, |p l − | 0 , is: where we have only considered the leading term in the expansion of the denominator in the appropriate small quantities. We want to perform the d|p l − | integration in the integral above, which will force |p l − | → |p l − | 0 . The density integral (γ eq,(λ=−1) (N → l − h + )) then becomes We now wish to do the angular integration and change the variable sin(θ )dθ = −d cos(θ ).
where we have expanded the square root for |p N | > |p l − | 0 which is true for most angles and called cos(θ ) = u. Note that the denominator remains always positive. Relabelling v = 1 − u, the integral above becomes To simplify this integral we will split it up into two regimes where different terms in the denominator are dominant, where α 4/3 /2 denotes the point where v 3 starts to dominate over the other terms in the denominator. We then have Substituting this into the expression for the γ eq integral we obtain the expression for the inverse of the heavy neutrino energy is approximated below up to second order, keeping all necessary terms.
substituting this into the integral and multiplying out with the expression in the round brackets we obtain.
The equilibrium distribution can be expressed as we consider terms upto second order in the exponential series. The gamma integral becomes.
in order to simplify these integrals we notice that there is a general expression, γ eq,(λ=−1) (N → l − h + ) = β I 1 − I 2 + I 3 where n = 1, 2, 3 and again expanding the exponential to second order. We thus obtain an expression for the integral I n : where we have employed a change of variable |p N |/T = x. Substituting in the different values for n, we obtain the solutions for J n , after performing the appropriate integrations: with E n (x) = ∞ 1 dy y −n e −x y is the exponential integral function [68], and we have the values Therefore to obtain the expression for the γ eq integral we recall that, This implies that the γ eq integral (to linear order in B 0 ) becomes We next proceed to obtain the expression of the pertinent γ eqintegral for the reverse process l − h + → N . The steps will parallel those of the previous calculation, the only difference is that now one should make the substitution f eq N → f eq l − f eq h + in (164). We thus have On using energy and helicity (λ = −1) conservation in the reaction l − h + → N , we observe that the numerator of the fraction in the exponent in (185) can be replaced by the energy of the RHN E which, upon substitution in (184) and comparison with (164), implies the reciprocity (chemical equilibrium) relation for the thermally averaged decay rates, γ eq,(λ=−1) (l − h + → N ) = γ eq,(λ=−1) (N → l − h + ), (187) in the presence of CPTV background B 0 = 0. The results for the decay and reverse processes N ← → l + h − will be analogous to those of the previous calculations but with a change in the sign of B 0 , due to the opposite helicity λ = +1 involved in those processes. The results for all thermally averaged decay rates are summarized below: Equation (188) implies the generation of a lepton asymmetry between the decay channels (9) of Fig. 1 at tree level only when B 0 = 0, due to the difference in the respective decay rates.

Padé approximants method
It is often possible to increase our knowledge of a function f (z) beyond the region of convergence of its Taylor series using the method of Padé approximants. The Padé approximation [54] can be considered as follows: given a function f (z) (with a Taylor expansion around z = 0), and two nonnegative integers m, n ≥ 0, the Padé approximant [m/n] f (z) is provided by the function We can Taylor expand P n m (z) to order z n+m and equate the expression to T m+n (z). (As n, m → ∞ often P n m (z) → f (z) even when the Taylor series for f (z) is divergent.) Let us consider a m × m matrix A defined by A i j = c n+i− j (1 i, j m). From the matching of the series the b i satisfy the matrix equation The coefficients a j , b j in (189) are uniquely determined, provided we normalise the zeroth order term in the denominator to one. The coefficients a i are determined by the set of equations A common procedure is to examine the convergence of the sequence P J 0 , P 1+J 1 , P 2+J 2 , P 3+J 3 , . . . with n = m + J . We shall use the J = 0 sequence known as the diagonal sequence.
This method will be applied to our system of Boltzmann equations to extrapolate their solution from z 1, where the equations are derived analytically, to the z 1 case. It is understood that although above we considered a Taylor expansion about z = 0 (which was assumed to be in the region of analyticity of f (z)), the discussion can be straightforwardly extended for Taylor expansions about any other point inside the region of analyticity of f (z). The application of Padé approximants and justification of Padé approximants are well described in [70].