Dynamical Evolution of Gravitational Leptogenesis

Radiatively-induced gravitational leptogenesis is a potential mechanism to explain the observed matter-antimatter asymmetry of the universe. Gravitational tidal effects at the quantum loop level modify the dynamics of the leptons in curved spacetime and may be encoded in a low-energy effective action Seff. It has been shown in previous work how in a high-scale BSM theory the CP odd curvature-induced interactions in Seff modify the dispersion relations of leptons and antileptons differently in an expanding universe, giving rise to an effective chemical potential and a non-vanishing equilibrium lepton-antilepton asymmetry. In this paper, the CP even curvature interactions are shown to break lepton number current conservation and modify the evolution of the lepton number density as the universe expands. These effects are implemented in a generalised Boltzmann equation and used to trace the dynamical evolution of the lepton number density in different cosmological scenarios. The theory predicts a potentially significant gravitationally-induced lepton-antilepton asymmetry at very early times in the evolution of the universe.


Introduction
The origin of matter-antimatter asymmetry in the universe is one of the most important outstanding issues in cosmology. Radiatively-induced gravitational leptogenesis (RIGL) [1][2][3] is a particularly elegant and economical mechanism to generate a lepton number asymmetry in the early universe with the potential to explain the presently observed baryon-to-photon ratio η = (n B − nB)/n γ 6 × 10 −10 .
The fundamental idea is that tidal gravitational effects at the quantum loop level can induce an asymmetry in the propagation of leptons and antileptons. This may be interpreted as the generation of an effective chemical potential for lepton number, which in the high temperature environment of the early universe produces a non-vanishing lepton number asymmetry in (quasi-)equilibrium.
The three basic conditions for the generation of a baryon or lepton number asymmetry have been known since the early work of Sakharov [4], viz. (i) a source of B or L violation; (ii) C and CP non-conservation; (iii) out-of-equilibrium B or L violating reactions. Gravitational leptogenesis [5] circumvents the third condition by exploiting the time dependence of the gravitational background in the expanding early universe, allowing a lepton asymmetry to be induced essentially in equilibrium. This may subsequently be converted to a baryon asymmetry at lower temperatures through electroweak sphaleron processes [6]. Beyond the Standard Model (BSM) physics remains necessary in this scenario to satisfy the first two criteria.
The RIGL mechanism is most easily understood in terms of an effective action for the light neutrinos ν L in curved spacetime [7]. This extends the original effective action for tidal gravitational effects in QED [8] to the neutrino sector [9], To leading order in the gravitational field strength, this may be written as suppressing flavour labels. The direct curvature couplings in this Lagrangian effectively violate the strong equivalence principle, allowing gravity to distinguish between matter and antimatter. They are generated perturbatively at one and two-loop level in the standard model or, as required for realistic leptogenesis, a BSM theory characterised by a high mass scale M and exhibiting CP non-conservation, and can be understood qualitatively as the result of tidal gravitational forces acting on the extended cloud of virtual particles in the quantum loops. The effective action (3.1) is a faithful description of the dynamics for weak gravitational fields R/M 2 < 1 (where R denotes a typical curvature) and for sufficiently low energies, as we specify carefully later.
The physics implications of these new gravitational interactions is quite different depending on whether the corresponding operators are CP even or CP odd. The only CP odd operator is the one with coefficient b and this is responsible for generating a non-vanishing equilibrium lepton number density n eq L . Noting that J µ = ν L γ µ ν L is the lepton (neutrino) number current, its µ = 0 component is the corresponding charge and the operator in (1.1) may be understood as the introduction of an associated chemical potential µ = bṘ, which is non-vanishing in a time-varying gravitational field. At finite temperature, this gives rise to an asymmetry n eq L = 1 3 bṘ T 2 . The two-loop BSM calculation of this operator, its interpretation in a modified Boltzmann equation, and its implications for cosmology have been described in detail in [1][2][3].
The new feature in this paper is an analysis of the role of the CP even operators with coefficientsã, c andd in driving the evolution of the lepton asymmetry at very early times. One way to understand these effects is to note 1 that the effective action (1.1) implies that the lepton number current is not conserved. In fact, where a =ã − 1 2d andb = 1 2ã + c − 1 2d . As we show in section 5, this implies the following equation for the dynamical evolution of the lepton number density n L in a FRW universe, where H is the Hubble parameter. These new curvature terms crucially modify the evolution of n L at very early times. Whether this acts to amplify or suppress the magnitude of n L depends on the signs and relative magnitudes of the coefficients, especially a, which are not arbitrary in RIGL but are determined by the fundamental BSM theory.
The combined effect of these two distinct mechanisms -the generation of an equilibrium asymmetry n eq L by the CP odd b operator and the modified evolution of n L by 1 Non-conservation of the lepton number current in a Lagrangian with an arbitrarily included operator of type R ν L γ.Dν L was observed in an interesting recent paper [10], which was an important motivation for the present work. The application to leptogenesis proposed in [10] is however very different from the mechanism developed here. the CP evenã, c,d operators -in a cosmological setting is clearly expressed in terms of a generalised Boltzmann equation. This is traditionally written for the lepton-to photon ratio N L = n L /n γ expressed as a function of inverse temperature z = M/T . In section 5, we derive this new Boltzmann equation for gravitational leptogenesis as, Here, W(z) reflects the curvature-induced evolution terms in (1.3), while W (z) is determined by the L-violating interactions in the BSM theory. Normally interpreted (see for example [11]) as a 'washout' term in leptogenesis mechanisms for which N eq L is zero, in our theory it plays the quite different role of driving the lepton asymmetry to its equilibrium value. Putting everything together, the entire evolution predicted by the RIGL Boltzmann equation is shown in Fig. 11. Essentially we find three stages: a very early high temperature phase in which the new evolution term W(z) keeps N L below equilibrium; followed by a phase where W (z) > W(z) and N L is driven to N eq L ; and finally decoupling when W (z) becomes too weak to hold N L to N eq L (i.e. the L-violating reactions are too slow compared to the Hubble expansion to maintain equilibrium) and it decouples leaving the final constant asymmetry predicted by this theory.
Clearly the transitions between these stages depend on the dynamical balance between the curvature and temperature dependent rates W (z) and W(z), and N eq L (z), as the universe expands and cools. In particular, the sign of W(z), controlled by the coefficient a in the effective action, is key to whether the new evolution term amplifies or suppresses the lepton asymmetry at early times.
We emphasise that these mechanisms for the generation and evolution of leptogenesis are very general and depend only on the SEP violating effective action (1.1). The philosophy of RIGL is that this effective action is generated automatically, and necessarily, in any given BSM theory incorporating CP violation. Nothing is added by hand, and all the coefficients are calculable in the fundamental theory. The overall picture of gravitational leptogenesis presented here is, however, not dependent on any particular choice of BSM theory.
The remainder of the paper is organised as follows. For clarity, we present our general theory using a particularly well-motivated BSM theory [12] in which the standard model is augmented with heavy right-handed sterile neutrinos ν R . This model provides an explanation, via the see-saw mechanism, for the light neutrino masses. As we see in section 6, where we discuss the implications for cosmology, the parameter bounds set by the experimental values for the neutrino masses severely constrain the resulting gravitational leptogenesis predictions in this model. The BSM Lagrangian is introduced in section 2, where explicit calculations of the one-loop diagrams determining the coefficientsã, c,d are presented. The two-loop coefficient b was previously calculated in [2]. These Feynman diagrams for the neutrino self-energies are matched to the effective action in section 3.
The new dynamical evolution mechanism is presented in sections 4 and 5. First, in section 4, the wave solutions of the equation of motion derived from the effective action are found using the eikonal approximation [7]. At leading order, the phase is determined by the CP odd b operator -this modifies the dispersion relation in a different way for the light neutrinos and antineutrinos, giving the essential asymmetry which ultimately leads to a non-vanishing n eq L in thermal equilibrium. The sub-leading eikonal order describes the time dependence of the wave amplitude, which is interpreted as lepton number density. This provides an independent derivation of the evolution equation (1.3). The derivation from current non-conservation is given in section 5. Section 5 also contains the derivation and initial interpretation of the gravitationally modified Boltzmann equation.
Finally, in section 6, we return to the full BSM theory and study gravitational leptogenesis quantitatively in different cosmological scenarios [3,5] -first in the standard radiation-dominated era where FRW cosmology is well-established, then during conventional reheating where the effective equation of state is characterised by 0 < w < 1/3, and finally in a more speculative post-inflationary era in which the usual reheating phase is replaced by a period in which the expansion is dominated by inflaton dynamics with w > 1/3. We track the gravitationally-induced evolution of the lepton number density from very early times through to decoupling and discuss the conditions and parameter choices for which the observed value of the baryon-to-photon ratio η may be realised.
Analytic results complementing the numerical plots of the evolution of N L (z) in section 6 are given in Appendix A.

Fundamental BSM Theory and Leptogenesis
The requirements for a fundamental theory which could realise our mechanism of radiatively-induced gravitational leptogenesis include CP violation active at a high energy scale where the spacetime curvature is sufficiently large, and a mechanism for lepton number violation. These are naturally incorporated in a BSM theory with three sterile neutrinos (i.e. having no interactions with the SM gauge fields) with a hierarchy of large Majorana masses, coupling to the usual left-handed leptons via the SM Higgs field. Independently of the question of leptogenesis, there are compelling reasons to augment the standard model in this way, with the see-saw mechanism then explaining the existence of the non-zero masses of the light neutrinos.

BSM Lagrangian and lepton number violation
The fundamental action, including the right-handed neutrinos, is then: Here, ν α R (α = 1, . . . 3) are the right-handed neutrinos, with Majorana mass matrix M αβ , which we take to be diagonal. i L (i = 1 . . . 3) are the SM lepton doublets and φ is the Higgs field. 2 The complex Yukawa couplings λ iα introduce CP violation into the theory, as required by the second Sakharov condition for leptogenesis.
With a Higgs VEV v, the usual see-saw mechanism gives rise to three light Majorana neutrinos with mass matrix along with the heavy sterile Majorana neutrinos with masses M α . This is diagonalised by the PMNS matrix U such that U T m ν U = m diag , with the corresponding relation between the mass and flavour eigenstates. Except where explicitly mentioned, we neglect these light neutrino masses in the rest of this paper.
The coupling to gravity is through the connection alone, as required for a Lagrangian satisfying the strong equivalence principle, with the covariant derivative acting on spinors being D µ = ∂ µ − i 4 ω µab σ ab , where ω µab is the spin connection and To see how the gravitational interactions induce an asymmetry in the propagation of leptons and antileptons, it is convenient to write S explicitly in terms of the fields 2 Our notation here is φ r = rsφ * s whereφ is the usual Higgs doublet giving mass to the lower fields in the SU (2) lepton doublets. 3 We use Greek indices µ, ν, . . . to refer to coordinates in the curved spacetime and Latin indices a, b, . . . to refer to the local Lorentz frame at each point. Our metric and curvature sign conventions are [S1][S2][S3] = − + + in the terminology of [13] and we use the Dirac matrix definitions of [14]. L , ν R and their charge conjugates c L ≡ ( L ) c and ν c R . Note that c L is a right-handed field, since in general (ψ L ) c = (ψ c ) R . We then have: The corresponding propagators are denoted by while for the right-handed fields we have both 'charge conserving' and 'charge violating' propagators, Note that in curved spacetime, translation invariance is lost and the propagators are not simply functions of the coordinate difference |x − y| as in flat spacetime. This becomes crucial below. In flat spacetime, the momentum space propagators are (neglecting the light neutrino masses, and writing G(p) for the physical Higgs component only) and (2.9) With a non-vanishing Majorana mass, the propagator S × α allows reactions which violate lepton number by 2 units, viz. ν L H ↔ ν c L H and ν L ν L ↔ H H, illustrated in Fig. 1. Both diagrams depend on the Yukawa coupling factor through λ S × λ T = α λ iα S × α λ T αj . This is the source of the lepton number violation which is required by the first Sakharov condition.

Lepton-antilepton asymmetry in curved spacetime
To implement the mechanism of radiatively-induced gravitational leptogenesis, we need to show that the propagation of leptons and antileptons is different in a gravitational field. Specifically, we find that at loop level, the self-energies Σ and Σ c for the leptons and antileptons differ when translation invariance no longer holds, leading to distinct dispersion relations. Together with the lepton number violating reactions in Fig. 1, this enables a lepton-antilepton asymmetry to be generated in thermal equilibrium in an expanding universe.
One-loop self-energy diagram for the light ν L neutrinos with an intermediate charge-conserving ν R propagator S. The fields at the vertices are shown explicitly for comparison with the Lagrangian (2.3).
At one loop, there is a single self-energy diagram involving the right-handed neutrinos, shown in Fig. 2. Evidently, at this order there is no corresponding diagram involving the charge-violating propagator. For clarity, we explicitly show the fields at the vertices in the diagram to allow the Yukawa couplings to be easily read off from the Lagrangian in the form (2.3). The self-energy is therefore, Note that depending on the Yukawa couplings, this can induce lepton flavour-changing processes dependent on the gravitational field.
The corresponding self-energy diagram for the antileptons is evidently given by Since we are interested in the violation of total lepton number, we can trace over the light lepton flavours, leaving As we calculate below, these one-loop self-energies contribute to the CP even terms in the effective Lagrangian (3.1), though not to the CP odd term generating the matterantimatter asymmetry. Figure 3. Two-loop self-energy diagrams for the light ν L neutrinos giving rise to a leptonantilepton asymmetry in curved spacetime.
The situation is different at two loops [1,2]. Here, we have the three self-energy diagrams illustrated in Fig. 3. The corresponding self-energies are: from the 'nested' diagram with two S propagators, and and for the 'nested' and 'overlapping' diagrams with two S × propagators. Here, G(x, y) generically denotes the appropriate component of the Higgs field propagator. Note that there is no overlapping-type diagram with two S propagators -this follows from the complex SU (2) L doublet nature of the Higgs propagator.
In curved spacetime, all these diagrams produce an asymmetry in Σ ij − Σ c ij but, as we now show, only the diagrams with two S × propagators produce a total leptonantilepton asymmetry, unsurprisingly since these are the diagrams involving the chargeviolating propagators S × . Tracing over the light flavours as before, and focusing on the Yukawa couplings, we have, with the integral factor read off from (2.13), and tr Σ Noting further that Im (λ † λ) βα (λ † λ) βα is antisymmetric in α, β, we can replace the integral factors above by I [αβ] and I [αβ] .
This brings us to a key point [1,2]. To obtain a lepton-antilepton asymmetry, the integral factors in (2.17), (2.18) must have a contribution which is antisymmetric in α, β. Now, in flat spacetime, translation invariance implies that the propagators are all functions of the difference of coordinates, i.e. ∆(x, y) → ∆(x − y), etc. But then, making a suitable change of variables on z, z , we can readily show 4 that the factors I (2) αβ (x, y) and I (3) αβ (x, y) are symmetric in α, β. This is, however, special to flat spacetime. It accords with the general theorems presented in [2] that CPT invariance together with translation invariance ensures that the propagation of particles and antiparticles is identical. It is no longer true in curved spacetime. In this case, we have found that at two loops there is an asymmetry in the self-energies of the light leptons and antileptons given by [αβ] + I . (2.19) Provided that the Majorana masses of the sterile neutrinos are non-degenerate, we therefore have a gravitational mechanism for establishing a lepton-antilepton asymmetry through two-loop radiative corrections to the light lepton propagators. In subsequent sections, we show how this leads to a non-vanishing net lepton number density in thermal equilibrium in an expanding universe.

Weak gravitational field expansion
The gravitational dynamics of the light neutrinos in this model is captured by the effective action S ef f discussed in section 3. The nature of the interactions in S ef f is constrained by general principles including local Lorentz invariance and hermiticity, with the leading terms at low energy and first order in the curvature shown in (3.1). This includes both CP conserving and violating operators, reflecting the potential breaking of CP invariance by the Yukawa couplings in the fundamental theory.
It remains only to determine the coefficients of these operators in terms of the parameters and couplings of the fundamental theory (2.1). The simplest method is to 4 For example, for diagram (2), translation invariance implies under the change of dummy variables u = x + y − z and u = x + y − z.
compare the predictions of the two theories for the light neutrino-graviton vertex in the weak gravitational field limit. In fact, as explained in [7], it suffices to consider only conformally flat metrics for this purpose, but note that this is not a restriction on the validity of S ef f for general curved spacetimes.
So, writing the metric in the weak-field approximation as for which the expansion of the curvatures at O(h) is together with conformally rescaled fields (in n dimensions) where we have taken the Higgs field to be a conformal scalar, i.e. ξ = 1/6.  From this, we can read off the elementary vertices for the h couplings to the left and right-handed neutrinos and the Higgs field. A representative set are shown in Fig. 4. Note particularly the factor of O(n − 4) entering all the Yukawa vertices -these will give a finite contribution when inserted into UV divergent diagrams. The next step is to incorporate these into the one and two-loop self-energy diagrams for the light neutrinos, effectively integrating out the heavy fields, and subsequently comparing to the equivalent vertices derived directly from the effective Lagrangian. This enables us to match coefficients and determine the values of the effective couplings in S ef f in this model. Consider first the one-loop self-energy diagram in Fig. 5, without an h insertion. Evaluating in flat spacetime, the self-energy is 5 Figure 6. One-loop self-energy diagrams with h insertions on the Yukawa vertices. 5 Evaluating the finite part of the self-energy diagram gives

Feynman diagrams
where we have expressed the result in terms of the dimensionally renormalised couplings and µ is the corresponding mass scale.
We can then immediately read off the contribution of the diagrams with h insertions on the Yukawa vertices, shown in Fig. 6. These give 6 Now consider the h insertion on the right-handed neutrino propagator, Fig. 7. From the sum of the two diagrams, we find (2.26) and evaluating in the standard way introducing Feynman parameters, we find At the mass scales of interest for leptogenesis in this model, m 2 H M 2 α , so provided the diagram is IR safe we can immediately set m H → 0 at this point. This is indeed the case. Next, in order to match coefficients with those in the effective Lagrangian, it is sufficient [7] to consider the leading terms in an expansion in terms of p 2 , p.q and q 2 , neglecting terms of O(p 4 /M 2 α ) etc. After some calculation, we find Note that the term linear in momentum in (2.29) precisely cancels the contribution (2.25) from the Yukawa insertion diagrams. For the Higgs insertion diagram shown in Fig. 8, we have Here, we do have to be careful since the integral is divergent as m H → 0. Evaluating the divergent terms, and cancelling against the overall m 2 H vertex factor, we find that this diagram gives a non-vanishing finite contribution in the m H → 0 limit, viz.
where again we are keeping only terms of O(p 2 /M 2 α ) etc.
Finally, collecting results, we find the total contribution of the h insertions into the one-loop light neutrino self-energy, at leading order in the momentum expansion, is 33) This allows the coefficients of the three CP conserving operators in the effective Lagrangian to be determined by matching the h ν L ν L vertex. Figure 9. Two-loop self-energy diagrams for the light leptons with an h insertion on the sterile neutrino α propagator, corresponding to the two lower diagrams in Fig. 3. Similar diagrams with the h insertion on the sterile neutrino β propagator are also to be included, along with insertions on the Higgs propagator and Yukawa vertices [2]. The diagrams with two S α propagators have an h ν R ν c R vertex while those with two S × α propagators have an h ν c R ν R vertex, as seen from (2.23) and Fig. 4.
In order to generate the vital CP violating operator in (3.1), we need to go to two loops, where we have already demonstrated a difference in the propagation of leptons and antileptons. The key diagrams here are shown in Fig. 9, corresponding to the two lower diagrams in Fig. 3. Note that the h insertion on a sterile neutrino propagator S × α gives rise to two contributions, one from the vertex −iM α h ν R ν c R associated with two S α propagators and another from −iM α h ν c R ν R with two S × α propagators. These insertions can be made in both nested and overlapping diagrams and on the α and β propagators. In addition, we need to consider insertions on the Higgs propagators and on the Yukawa vertices, which, as we have seen at one loop, can in principle contribute despite the h vertex being of O(n − 4) .
The resulting calculations are extensive and perhaps surprisingly stretch the limits of known analytic methods. The hardest diagram is essentially a two-loop 3-point triangle diagram with arbitrary external momenta. These calculations were carried out in detail in [2], although in that reference we were unable to complete the evaluation of this final diagram. This leaves some uncertainty in the heavy mass hierarchy dependence of the final result.
To extract the relevant coefficient ('b') in the effective Lagrangian, we need only isolate the term involving q 2 / q in the momentum expansion of Σ 2 loop ij . In [2] we found, where the integral factor I αβ was indeed shown to have an antisymmetric part. The result in the large hierarchy limit M β M α (see [2] for the full result) is up to an O(1) numerical factor. This shows that the gravitationally induced leptonantilepton asymmetry discussed above is realised explicitly at two-loop order, and consequently generates a non-vanishing coefficient for the CP violating operator in the effective Lagrangian describing the low-energy dynamics in this model.
The question of whether the hierarchy parameter p is 0 or 1 was left unresolved in [2]. The most natural result, which holds in all the diagrams we calculated to a conclusion, is p = 0. This would accord with the expected decoupling of heavy mass intermediate states in the Feynman diagram. 7 Nevertheless, in calculating the arbitarymomentum triangle diagram, we found contributions with p = 1 which would require a remarkable cancellation if they were to be absent in I [αβ] . For this reason, we retain the possibility that p could be 1 in what follows since, as we see in the cosmological scenarios, the presence of a large sterile mass hierachy dependence would significantly affect the ultimate prediction for the baryon-to-photon ratio in this particular BSM theory. 7 However, note that because of the external factor of M β from the h coupling, in order to find a p = 1 hierarchy dependence we only require the relevant diagram to tend to a constant or grow no faster than log M β in the large M β limit. As an illustration that such a logarithmic dependence on the mass of a heavy neutrino propagator may arise in an individual Feynman diagram, note the asymptotic behaviour of the UV divergent one-loop self-energy diagram in Fig. 5, where we find a log M α dependence in the M 2 α p 2 limit (see footnote 5).

Effective Lagrangian for Gravitational Leptogenesis
In this section, we construct an effective action which encodes the tidal curvature effects on the propagation of the light neutrinos and fix the coefficients of the induced operators by comparing with the explicit calculations of the self-energy loop diagrams in section 2.

Effective Lagrangian
The gravitational dynamics of the light neutrinos, to leading order in the curvature, is encoded in the general effective action, With real coefficientsã, b, c,d, each term in (3.1) is individually hermitian. We have also allowed for gravitationally-induced flavour-changing effective interactions, with the coefficients carrying flavour labels.
This effective Lagrangian arises quite generally as the low-energy limit of a fundamental UV-complete theory, characterised by some high mass scale M . The coefficients a, b, c,d are then of O(1/M 2 ), so the Lagrangian is a weak gravitational field expansion in R/M 2 , where R denotes a typical curvature. Taking the BSM theory of section 2 as the fundamental theory, the mass M is identified with the Majorana mass M α of the heavy sterile neutrinos.
It is also a low-energy expansion, since we have retained only leading-order terms in derivatives. In the papers [15,16], which analysed the full energy dependence of UVcomplete quantum field theories coupled to gravity, the relevant expansion parameter was shown to be E √ R/M 2 , for a typical particle energy E. This will be an important factor in the application of our results for leptogenesis in realistic cosmological settings, and we will return to a careful discussion of this point later.
The presence of direct couplings to the curvature means the Lagrangian (3.1) violates the strong equivalence principle and cannot be regarded as a fundamental theory in its own right. In particular, in the absence of an embedding into a UV-complete theory, the theory described by this effective Lagrangian is not causal. This issue has been explored extensively in the series of papers [16][17][18][19][20] on the realisation of causality and unitarity in QED in curved spacetime. 8 The properties of the operators in (3.1) under the discrete symmetries C, P, T is important, and has been described carefully in [7]. Note that in curved spacetime, the symmetries P, T are defined only in the local Lorentz frame, or tangent space, at each point in spacetime. Crucially, the operatorsã, c,d are all CP even, while only the b operator is CP odd. All are CPT even. 9

Weak gravitational field expansion
The necessary formalism to match the coefficients in the effective Lagrangian to the couplings and masses of a fundamental theory was developed extensively in [7] and we need only summarise the essential results here. In [7], we considered the low-energy limit of the standard model coupled to gravity; here, we are interested in the effective Lagrangian for the BSM model of section 2.
To facilitate comparison with [7], and to simplify the application to the eikonal expansion in section 5, it is convenient to rewrite the effective action (3.1) in the form, To relate S d and Sd we use the identities,

4)
8 See also [21] and the discussion of causality in explicitly CPT violating theories in [22]. 9 This is in contrast to a model where instead of the curvatures, the bilinear neutrino operators are multiplied by constant parameters, as in the Lorentz and CPT violating standard model extension introduced in [23]. Importantly, the Lorentz vector and tensor operators i ν L γ a ← → D b ν L and ν L γ a ν L are CPT even and odd respectively.
As discussed in section 2, the simplest approach to determine the couplings in the effective Lagrangian is to compare with the fundamental theory in the weak gravitational field approximation. Using the weak field expansion of the curvatures in (2.21), we can show that the operators in (3.3) give rise to an effective h ν L i ν j L vertex V ij as shown in Fig. 10 which, after some calculation [7], can be expressed in momentum space as 10 (3.6) 10 The vertex is subject to an important constraint from unitarity. Note first, suppressing flavour labels, that V(p, q) = iM(p, q), where This amplitude is required to satisfy the unitarity relation, Writing the general expression, this implies the following relations (assuming all but φ are purely real), which are satisfied by the expansion (3.6).

Matching conditions
We can now determine the effective Lagrangian coefficients for the fundamental BSM theory in section 2 by comparing the effective vertex V ij (p, q) in (3.6) with the explicit Feynman diagram results in (2.33) and (2.34). With the identification, we can read off the following values of the coefficients: It is easily checked that with these identifications, (2.33) satisfies the unitarity conditions in footnote 10. Two features of these coefficients deserve comment. First, note that d = −4c. As we see below, this will imply the cancellation of these coefficients from the light neutrino equation of motion in the effective theory. Remarkably, precisely the same relation is found from the light neutrino self-energies in the standard model, as shown in [7]. Second, again as in the standard model, we find the sign of the coefficient a is negative. This translates directly into the evolution of lepton number in the cosmological models.
As already noted, there is no contribution to the coefficient b ij of the CP violating operator at one loop. From (2.34) we find its leading-order contribution, with This completes the identification of the coefficients in (3.1) for the specific BSM theory described in section 2. In the following sections, we develop the theory of gravitational leptogenesis based on the general effective Lagrangian.

Lepton Number Evolution and the Eikonal Expansion
We now come to the key theoretical development in this paper, the investigation of gravitational leptogenesis in the radiatively-induced effective Lagrangian (3.1). We do this in two complementary ways -first, from a detailed analysis of the eikonal expansion of the loop-corrected light neutrino equation of motion, then, in the following section, using operator methods and non-conservation of the lepton number current.
The four operators in S ef f play different roles in leptogenesis. The b operator ∂ µ R ν L γ µ ν L acts as a chemical potential for lepton number and allows a leptonantilepton asymmetry, i.e. a net lepton number, to arise in thermal (quasi-)equilibrium. The remaining operators modify the evolution of the lepton number in time, giving a dynamical amplification or damping depending on the signs of the coefficients a, c, d.

Eikonal expansion for the light neutrinos
The propagation of the light neutrinos in a background gravitational field is governed by the equation of motion derived from the effective Lagrangian. 11 This is written most easily with the parametrisation in (3.3) and we find We are viewing S ef f as being generated by radiative corrections, so we work to consistent perturbative order. Then, since γ.Dν L = O(λ 2 ), to this order we must omit the pre-factor (2cR + dD 2 ) which gives terms of O(λ 4 ) which we have not computed in the other coefficients. The perturbatively consistent equation of motion therefore reduces to The eikonal expansion consists of writing the field as the product of a slowly-varying amplitude A and a rapidly varying phase Θ. In order to keep track of the relative orders in the eikonal expansion, we temporarily introduce a formal counting parameter and write the ansatz where u L is the appropriate spinor wave function. The wave vector is defined as p µ = ∂ µ Θ. 11 For clarity, we present the formal developments in sections 4 and 5 for a single light neutrino flavour. Alternatively, the formulae given here may be viewed as matrix expressions with the flavour indices suppressed in the coefficients a ij etc. and currents J µ ij = ν L i γ µ ν j L . We have also taken the light neutrino mass to zero in this subsection.
The inclusion of an O( ) correction in the phase Θ is novel and is required to accomodate the CP violating b operator. As explained in reference [7], this operator induces a phase modulation of the wave solution, which in the particle interpretation produces a linear shift in the energy in the dispersion relation, with opposite sign for the neutrinos and antineutrinos. This difference in the propagation of particles and antiparticles in a gravitational field is the origin of gravitational leptogenesis.
The simplest route to a complete solution of the eikonal equation is to first act on (4.2) with γ.D to remove the explicit dependence on the gamma matrices and produce a second-order wave equation. After some manipulation of covariant derivatives acting on spinors, using the identities (3.4) and (3.5), and again keeping terms to consistent perturbative order only, we eventually find Inserting the ansatz (4.3) and collecting terms of the same order in then gives, at and at O(1/ ), where we have defined k µ = ∂ µ θ. Higher orders in give sub-leading corrections to the amplitude which do not concern us here.
The first equation determines the one-loop curvature induced modification to the dispersion relation when re-expressed in terms of the true wave vector p µ = k µ + ∂ µ α. For the second, note that the imaginary terms cancel at O(λ 2 ) with the choice α = −bR. In fact, the entire dependence of the wave operator on the coefficient b is accounted for by this addition to the phase factor. Now, in the absence of the loop effects, the dispersion relation from (4.5) is simply k 2 = 0, from which we deduce k.D k ν = 0. This is just the geodesic equation, showing that the tangent vector k µ is parallel transported along its integral curve. The O(1/ ) terms may be separated into equations for the amplitude and spinor wave function, viz.
This states that the wave function u L is parallel transported along the geodesic, whereas the amplitude varies as the 'expansion'θ = − 1 2 D.k of the geodesic congruence, one of the optical scalars in the Raychoudhuri equations. This has the clear interpretation of the amplitude increasing as the congruence focuses. In the particle interpretation of the eikonal formalism, the number density of particles comprising the wave is proportional to the square of the amplitude, so the light neutrino number density varies as n ν ∼ A 2 .
Putting all this together, we can finally write the eikonal solution to the loopcorrected equation of motion in the form, where the dispersion relation is given in (4.5), while the amplitude and spinor wave function satisfy Neither k µ nor p µ satisfies the geodesic equation with the loop effects included. The evolution equations for both the amplitude and the wave function are also modified as shown, the latter depending on the spin generator σ ρσ .
In terms of the true wave vector p µ = ∂ µ (θ − bR) = k µ − b ∂ µ R, and including a light neutrino mass, the dispersion relation is This is the form appropriate for an interpretation in terms of energy and momentum, with the ± sign for the neutrino ν L and antineutrino ν c L respectively.

Interpretation in FRW spacetime
We now show how this eikonal solution may be interpreted in terms of the evolution of lepton number in a FRW spacetime.
The spatially flat Robertson-Walker metric, has non-vanishing Christoffel symbols, where H is the Hubble parameter. The FRW Ricci curvatures are where ρ is the energy density and w is the equation of state parameter, p = wρ. The energy conservation equation,ρ Also note the Friedmann equation gives the Hubble parameter from The geodesic equations for a massive free particle in this metric are simply solved by x i = constant, i.e. the particle is co-moving with the cosmological expansion. Letting k µ be the corresponding momentum, and with dispersion relation k 2 = m 2 , this implies k 0 = m, k i = 0.
The corresponding solution of the loop-corrected dispersion relation (4.5) Along this trajectory 12 the amplitude equation (4.9) evaluates as 12 This is not in general the tangent vector to a geodesic since, writing K µ = k µ + 2aR µν k ν so the dispersion relation is simply K 2 = m 2 , we can show The time evolution of the amplitude is therefore Recalling that in the eikonal formalism, the particle number density n ∼ A 2 , we therefore find the time evolution of the light neutrino number density n ν in a FRW spacetime is given by

Gravitational Leptogenesis and the Extended Boltzmann Equation
In this section, we provide an independent derivation of the evolution equation (4.21) for the lepton number density by showing that the lepton number current is not conserved when the CP even curvature interactions are included in the effective action. We then formulate a new generalised Boltzmann equation incorporating these gravitational effects.
However, in the FRW metric and taking K 0 = m, K i = 0 as in (4.18), this does simplify to Also note that in terms of K µ , equation (4.9) for the amplitude becomes where we have made use of the Bianchi identity D µ R µν = 1 2 ∂ ν R. This form is especially close to the current non-conservation equation derived in section 5.

Current non-conservation and lepton number evolution
For the free Dirac Lagrangian, the lepton number current J µ = ν L γ µ ν L is conserved, i.e. D µ J µ ∼ 0 (where ∼ 0 indicates zero up to terms vanishing by the equation of motion). In the presence of the curvature terms in the effective Lagrangian, however, the current is no longer conserved, implying a non-trivial time dependence of the lepton number.
The current non-conservation identity is derived using standard methods. Under a variation ν L → e iθ ν L , ν L → e −iθ ν L , the effective action transforms as where the total derivative term defines the current J µ , while for a non-conserved current the remainder ∆ is non-zero. Taking the terms in the effective action (3.3) in turn, we find Collecting terms, the current non-conservation equation is therefore , the pre-factor 2cR of the D µ J µ term must be omitted to consistent perturbative order. So, up to terms vanishing by the equation of motion, we find simply, Now in the FRW spacetime, using the formulae (4.13), (4.14), this current nonconservation equation becomes The charge density corresponding to the current is J 0 , which we identify as the lepton number density n L . In an isotropic universe, the spatial gradients in (5.5) vanish. We therefore find the time evolution equation for the lepton number density, This reproduces (4.21), derived using the eikonal formalism, with n L = n ν − n ν c .
To develop this, noting that dn L /dt = −3Hn L + O(λ 2 ), we can rewrite the prefactor term in (5.6) so that to O(λ 2 ), Substituting for the curvatures finally gives, Whether the radiative curvature corrections amplify or reduce the lepton number density as the universe evolves therefore depends on the sign of the combination of coefficients 2a − 2b(1 − 3w) .
For a radiation dominated FRW universe, w 1/3, with a deviation arising purely from the beta functions characterising the trace anomaly in the energy-momentum tensor, T µ µ = 0. With standard model fields, this gives (1 − 3w) 0.1. Since this is small, the dominant factor in (5.8) is the coefficient 2a. A negative a corresponds to damping of the lepton number with time, while a positive a implies an amplification. Recall that in the particular BSM model described in sections 2 and 3, we found in (3.8) that a is negative.
Other cosmological scenarios with (1 − 3w) = 0 are considered in section 6. At this point, however, note that in the BSM model we have d = −4c, so the coefficient b = a/2. So in this particular model, the combination 2a − 2b(1 − 3w) = a(1 + 3w). Curiously, this is negative unless (1 + 3w) < 0, which is precisely the condition for an accelerating, or inflationary, universe.
We return to (5.8) in section 5.3 below, where we incorporate it into the new Boltzmann equation for lepton number. First, we consider the CP violating b operator and the mechanism for leptogenesis.

Gravitational leptogenesis
The CP violating operator in the effective Lagrangian (3.1) may be written in FRW spacetime as since R = R(t). In this form, given that J 0 is the lepton number density n ν , we see that this may be interpreted as introducing a chemical potential µ = bṘ for lepton number.
Since the neutrinos interact via the Higgs field with the finite temperature medium in the early universe, this chemical potential biases the net lepton number n L = n ν −n ν c to produce a non-vanishing value n eq L in thermal equilibrium, with n eq L ∼ µT 2 .
The equilibrium is maintained by the lepton number violating reactions described in section 2.1. As well as the ∆L = 2 reactions ν L H ↔ ν c L H and ν L ν L ↔ HH considered there, there are further ∆L = 1 reactions involving other standard model fields which are described in detail elsewhere (see for example [11]. Of course, since we need the Ricci scalar to be varying in time to produce the chemical potential, we are not in true thermal equilibrium. However, provided the ∆L = 0 reactions are faster than the rate of changeṘ, the neutrinos and antineutrinos are in quasi-equilibrium and the formalism here is an accurate description. As discussed in section 2, the time-dependence inherent in R(t) is necessary to satisfy the third Sakharov condition and allow leptogenesis to occur.
To see this in detail, recall the dispersion relation (4.11). Rewriting in components, this reads and identifying energy as E = p 0 , we find From standard statistical mechanics, the equilibrium lepton number density is then given, for small µ/T , by Since this is already O(λ 4 ), we may neglect the O(a) corrections to E(|p|) in the integral in (5.13), and find to a good approximation, This may be compared to the photon number density n γ = 2ζ(3) π 2 T 3 or the entropy density s = 2π 2 45 g * s T 3 , where g * s is the effective number of degrees of freedom at the energy scale T , both of which are commonly used to normalise the lepton asymmetry.
To summarise, we have shown that radiative corrections in the BSM model of section 2 induce a non-vanishing lepton number density asymmetry in thermal equilibrium, given by with the coefficient b = tr b ij = O(λ 4 ) given by (3.9). This is the mechanism of radiatively-induced gravitational leptogenesis, first proposed in [1].

Boltzmann equation
The final step is to incorporate these two new effects into the Boltzmann equation describing the time, or temperature, dependence of the lepton number density.
We begin with the modification to the evolution term (5.8). It is traditional in discussions of leptogenesis to give results in terms of the ratio of lepton to photon number, N L = n L /n γ . Since n γ ∼ T 3 and T ∼ 1/a, where a(t) is the FRW scale parameter, we have We also usually express the evolution in terms of temperature, or more specifically z = M 1 /T , where in our BSM model we can take the mass scale M 1 to be the mass of the lightest sterile neutrino. Then, d/dz = (1/Hz)d/dt. Recalling that 8πG = 1/M 2 p defines the reduced Planck mass, we can finally re-express (5.8) in the form, Next, consider the modification due to the non-vanishing equilibrium number density n eq L . For this, we need some kinetic theory. We only sketch the key results we need, referring to standard texts [24] for the general theory. We also need to step back from only using the effective Lagrangian here, since we need to consider the specific lepton number violating reactions described in section 2 in the BSM model itself.
For a ∆L = 2 reaction such as ν L H ↔ ν c L H, the number density of neutrinos is given by the Boltzmann equation [25], with the corresponding result for dn ν c /dt.
The cross-section term above is re-expressed in terms of the thermally-averaged reaction rate for particle ν, viz. Γ ν ∆L=2 , given by Γ ν ∆L=2 = 1 n eq ν γ ∆L=2 = 1 n eq ν n eq ν n eq H σ|v| in terms of the u-average of the amplitude M(u, s) for the ∆L = 2 reaction. (See for example [3,11] and references therein for details and notation.) We may take n H n eq H and, using Maxwell-Boltzmann statistics for the neutrinos (i.e. neglecting the 'Fermi blocking' effect), approximate n eq ν n eq ν c e 2µ/T , (5.20) from (5.13) due to the non-vanishing chemical potential. It then follows that, with Next, since in the free theory, n ν and n ν c take their common equilibrium values, and referring to (5.13) for n eq L , we see that to O(λ 4 ), µ T n ν + n ν c n eq L .

(5.22)
We therefore find, Despite the technicalities of its derivation, (5.23) has a very straightforward physical interpretation, viz. that provided the ∆L = 2 reactions are active, the lepton number density n L is driven towards its equilibrium value n eq L .
Finally, re-expressing in terms of dN L /dz as before, we find where W = Γ/zH is essentially the ratio of the reaction rate to the expansion rate of the universe. In the full theory, W ∆L=2 is replaced by the complete rate term W incorporating the inverse sterile neutrino decay and ∆L = 1 reaction rates as well as the ∆L = 2 reactions considered in detail here.
The full Boltzmann equation is therefore given by combining (5.17) and (5.24) and we find with W given by and, using the expression (5.15) for n eq L together with (4.16) for the curvatureṘ, This is the key equation of the paper. The new features compared to conventional leptogenesis models are the radiatively-induced curvature dependent terms -the nonvanishing equilibrium asymmetry N eq L already presented in [1][2][3], and the new evolution term WN L derived here.
In the next section, we study this equation in detail in different cosmological settings using parametrisations of the reaction rate W and other couplings expressed in terms of neutrino mass parameters in the specific BSM model in section 2. Equation (5.25) is, however, far more general than this particular model and it is interesting at this point to give a general estimate of the temperature dependence of the various terms in a radiation-dominated cosmology, with w 1/3. The energy density ρ in this case is given by ρ = σT 4 , where σ = π 2 g * /30 with g * the effective number of degrees of freedom. Then, recalling that the loop coefficients in (5.26) are Note that in this case, both the curvature-dependent terms N eq L and W fall off rapidly as T 5 as the temperature falls and the universe expands. Conversely, they become increasingly important in the high-temperature, strong-curvature regime of the early universe.
The lepton number violating rate term W is conventionally known as the 'washout' factor, since in the absence of a non-vanishing N eq L this would dampen out any preexisting lepton asymmetry. In the gravitational leptogenesis theory described here, however, it plays a quite different role, driving the asymmetry towards its equilibrium value N eq L . The new evolution term W, on the other hand, can act at early times either to slow the drive to N eq L or as a source of amplification of any existing lepton asymmetry, depending on the coefficients in the effective Lagrangian. When these are generated by loop corrections, the sign of these contributions is not arbitrary but is determined by the dynamics of the fundamental BSM model. In our particular model, we found W > 0, but it is not clear whether there is some general principle enforcing this sign.
The full evolution of the lepton asymmetry N L derived from the Boltzmann equation (5.25) is illustrated in Fig. 11, which extends the picture already found in [3] to higher temperatures. In the temperature region (b), with z < 1, the scattering factor W (z) dominates and drives N L to its gravitationally-induced equilibrium value N eq L . At still higher temperatures, the 1/z 5 dependence of W(z) means we reach a regime with |W(z)| > W (z), labelled as (a) in the Figure, where N L is forced away from the equilibrium value, either reducing from N eq L if W(z) > 0 as shown, or increasing very sharply if W(z) < 0. For lower temperatures (region (c)), as the universe expands and the temperature falls, the lepton number violating rate W (z) eventually becomes too slow to maintain equilibrium and N L decouples, becoming essentially constant. Finally, around z 1, that is T M 1 , the rate W (z) has a resonance peak due to the sterile neutrino propagator in the diagrams of Fig. 1 and is sufficiently strong to force N L back towards N eq L (region (d)), before once more becoming negligibly small leaving N L to converge to its final, low-temperature, constant value as shown.
These features can be quantified using an exact analytic solution for the Boltzmann equation in the regions (a), (b), (c) for temperatures above the resonance regime T M 1 . This is given in Appendix A, where we show that in the ultra-high temperature region (a), the new evolution term induces a dependence N L ∼ 1/z 2 in the lepton asymmetry, moderating the sharp rise N eq L ∼ 1/z 5 in the equilibrium value.
The balance of these competing terms in the extended Boltzmann equation (5.26) depends on the detailed choice of masses and couplings in the original BSM model. In our case, these are constrained by the light lepton mass spectrum, since the theory also plays the role of generating neutrino masses via the see-saw mechanism with the heavy sterile neutrinos. In the following section, we study in detail how this works out in particular cosmological scenarios.

Evolution of Gravitational Leptogenesis in Cosmology
In this final section, we explore the consequences of the generalised Boltzmann equation (5.25) for the creation and evolution of the lepton number asymmetry in the early universe. As in [3], we consider two scenarios in detail -first, when leptogenesis occurs in a conventional radiation-dominated FRW spacetime and second, in a post-inflationary era preceding radiation dominance when the expansion is driven by a source with effective equation of state w > 1/3, including the extreme 'kination' scenario with w = 1.
As usual in leptogenesis models, we assume that the lepton asymmetry generated in the early, post-inflationary universe is subsequently converted to a baryon asymmetry by sphaleron processes at the electroweak scale [6,26]. The required relation for the baryon-to-photon ratio η is (see for example [27]) where C sph is the fraction of the lepton asymmetry converted into a baryon asymmetry by sphaleron processes, and f is a dilution factor accounting for photon production between leptogenesis and recombination. In this model, C sph = (8n + 4)/(22n + 13) with n = 3 fermion generations and f = 2387/86, leaving η 0.02|N L |. To achieve the observed value of η = 6 × 10 −10 , we therefore require a leptogenesis mechanism to yield a lepton asymmetry |N L | 10 −8 .

Radiation-dominated FRW cosmology
Our focus in this paper has been to derive and study the two radiatively-induced gravitational terms in the generalised Boltzmann equation (5.25). We have developed this in terms of the effective Lagrangian (3.1) which describes the dynamics of the light neutrinos in curved spacetime, including the contributions from loop diagrams with virtual sterile neutrinos ν α R . To describe the complete dynamics of leptogenesis in the BSM model (2.1), however, we also need to take into account the lepton number violating decays of the sterile neutrinos themselves.
The Boltzmann equations describing leptogenesis through out-of-equilibrium decays of the sterile neutrinos are well-known (see [11] and references therein for reviews), this being the original Fukugita-Yanagida model [12]. We can therefore easily include these along with our gravitational terms, giving the complete Boltzmann equations, Here, N ν R is the ratio of the number density of the sterile neutrino ν 1 R to photons (since it is the lightest of the ν α R that gives the biggest contribution to the Boltzmann equation for N L ) and N eq ν R (z) is its equilibrium value at temperature T = M 1 /z. The decay parameter is given as and ε 1 is a CP-violating parameter of O(λ 4 ) characterising the asymmetric contributions of the ∆L = 1 decays of the sterile neutrino ν 1 R to ν L H and ν c L H. Together with W , and unlike N eq L and W, these terms are non-zero in flat spacetime and we can neglect any gravitational contributions to them in (6.2).
A natural question at this point is why it is consistent to use expressions derived from the 'low-energy' effective Lagrangian, incorporating the effects of integrating out the virtual heavy sterile neutrinos in loop diagrams, at an energy scale where the sterile neutrinos are themselves dynamical. The resolution depends on the extra (curvature) scale in the gravitational theory.
Conventionally, effective Lagrangians in flat spacetime are valid for momenta such that p 2 < M 2 , where M is the characteristic mass scale of the fundamental theory, and the low-energy expansion is in the parameter p 2 /M 2 . As noted in section 3, however, it has been shown [15,16] that the relevant expansion parameter in the curved spacetime theory is instead E √ R/M 2 . Equivalently, using the typical energy E ∼ T , the lowenergy expansion here is valid for z > √ R/M 1 . This means that for R/M 2 1 1, as required for the validity of the weak gravitational field approximation, we may legitimately use the effective Lagrangian at scales T significantly above the sterile neutrino mass M 1 . Recalling that for a radiation-dominated FRW universe, R ∼ T 4 /M 2 p , a reasonable estimate for the validity of the effective Lagrangian in this case is therefore 13 z 3 10 −2 M 1 /M p . The corresponding weak field condition is a lesser constraint, z 2 M 1 /M p .

Neutrino parameters and the Boltzmann equations
Since the BSM model incorporates the see-saw mechanism for generating the light neutrino masses, the parameters are constrained and can be re-expressed in terms of the experimentally known neutrino masses. Here, we follow [11] and references therein, summarised for our purposes in [3], so we only quote a few essential relations without further motivation.
Next, we introduce the key parameter K characterising the strength of the Yukawa interactions 14 , To accommodate the conventional physics of see-saw neutrino masses, we choose M 1 ∼ 10 10 GeV here, with low values of K between around 1 and 5.
The sterile neutrino decay parameter is readily expressed in terms of K and Bessel functions as 14 K is equivalently defined as the ratio of the zero-temperature decay rate to the Hubble parameter at T = M 1 , which controls whether the ν 1 R decays are in equilibrium. Figure 12. Diagrams for the decay ν α R → ν i L H which contribute to the decay rate asymmetry factor ε α . At O(λ 4 ) the relevant contribution to the Γ(ν R → ν L H) decay rate arises from the interference of the tree and one-loop diagrams shown. (Two further loop diagrams, with a self-energy insertion on the incoming ν R line or on the outgoing ν L line, do not contribute to the asymmetry.) Similar diagrams, where the S α and S × α type ν R propagators are switched compared to the figure, give the decay rate for ν α R → ν i c L H. Reading off the vertices from the action (2.3) shows that the asymmetry depends on the combination is a kinematical factor from evaluating the diagrams.
The CP violating decay parameter controlling the contribution to the lepton asymmetry from the out-of-equilibrium ν α R decays is given by [28] ε α 3 16π This arises from the interference of the tree and one-loop diagrams for the decays ν α R → ν i L H and ν α R → ν i c L H shown in Fig. 12. Here, we only require the parameter ε 1 for the dominant decays of ν 1 R . Note that (6.8) then involves a linear combination of both Im(K 2 12 ) and Im(K 2 13 ). A reasonable range of values for ε 1 , consistent with the constraints from neutrino masses, is ε 1 ∼ 10 −5 − 10 −10 .
The final element of the conventional Boltzmann equation is the factor W (z) including the ∆L = 2 scatterings, ∆L = 1 reactions and inverse decays ν L H → ν R . In the usual theory this is known as the 'washout' term, since it removes any pre-existing lepton asymmetry. However, as we have seen, once we have a non-vanishing N eq L (z) it plays a quite different role, driving N L (z) towards this equilibrium value.
The full derivation of W (z) is quite lengthy and involves a number of subtleties which need not concern us here. A summary of the necessary results is given in [3]. The form of W (z) is especially simple for large and small z. In terms of the parameters defined above, we have with W (z 1) satisfying the same formula with the K 2 term omitted. Sincem 2 /m 2 * 2150, the effect of the K 2 term is negligible for the values considered here, and the coefficient of the asymptotic 1/z 2 dependence of W (z) is essentially the same for small and large z. In the region z 1, W (z) displays a resonance arising from the ν 1 R intermediate state in the scattering diagrams (see Fig. 1).
The complete temperature dependence of W (z) for low values of K is shown in Fig. 13, illustrating the characteristic 1/z 2 behaviour anticipated in section 5 together with the resonance enhancement at T M 1 . This brings us to the new gravitationally-induced terms in the Boltzmann equation (6.2). The equilibrium lepton number asymmetry N eq L has been derived and discussed in detail in our earlier work [1][2][3]. Here, from (5.15) and (3.9) we have With a sterile neutrino mass hierarchy M 1 < M 2 < M 3 , the dominant term is α = 1. Then, substituting forṘ from (4.16) and with the equation of state w 0.3 for a radiation-dominated background (including the trace anomaly), we find where we have used (2.35) for I [αβ] and collected the numerical terms into a single pre-factor.
In writing (6.11), we have retained the possibility of a hierarchy enhancement with p = 1 or the, perhaps more likely, form of I [αβ] with p = 0. With p = 1, the dominant contribution comes from the heaviest sterile neutrino ν 3 R , whereas for p = 0 it is ν 2 Notice that in this parametrisation, the small value of the Yukawa coupling constants are exchanged for an extra power of M 1 /M p compared to the estimates in (5.28), assuming K is chosen of O(1) as we do here.
As already observed, W(z) has a 1/z 5 dependence, characteristic of the gravitational terms. It therefore rises sharply for temperatures T > M 1 . Comparing W(z) with the coefficient W (z), which has a milder 1/z 2 temperature dependence, we can determine the value of z at which W(z) begins to dominate and force the lepton asymmetry away from its equilibrium value. With parameters chosen conventionally to reproduce light neutrino phenomenology, and with M 1 10 10 GeV and K 1 − 10, this crossover value is around z 10 −6 , which unfortunately lies outside the range of z where we can unambiguously rely on the low-energy expansion of the effective Lagrangian. It is, however, very model dependent and for that reason, and with this caveat, we include its effects in the Figures below where we illustrate the evolution of N L (z) with temperature.
Note especially that this crossover value of z is particularly sensitive to the observed light neutrino masses, which control the scale of W (z) in (6.9). This emphasises again the strong constraints on the gravitational leptogenesis mechanism when implemented in this particular BSM model, which is playing the dual role of generating the light neutrino masses through the see-saw mechanism. In section 6.2, we show how this value also depends on the equation of state parameter in generalised cosmological scenarios, and occurs for significantly lower temperatures for w > 1/3.

Evolution of the lepton number asymmetry
We now present a variety of plots of the evolution of the lepton number asymmetry N L (z) with temperature, found by numerically solving the coupled Boltzmann equations (6.2) with different assumptions and parameters in the BSM model.
To begin, in Fig. 14, we show the sterile neutrino density N ν R (z) and the lepton asymmetry N L (z) in the conventional model neglecting the gravitational contributions. The two options of starting at high temperature with N ν R (z) taking its equilibrium value, or with N ν R (z 1) 0, are shown separately. The effect of the decay term proportional to N ν R − N eq ν R in the Boltzmann equation for dN L (z)/dz is readily understood, becoming appreciable only when the decay rate D(z), which grows with z, becomes sufficiently large. The above equilibrium abundance of ν 1 R in the upper righthand plot in Fig. 14   The left-hand diagram corresponds to an initial condition with N ν R (z) = 0 while the righthand plot starts with N ν R (z) already at its equilibrium value. The lower figures show the corresponding absolute value of the lepton asymmetry induced by the out-of-equilibrium decays of ν 1 R . Notice the cusp in the left-hand plot, indicating that the asymmetry N L (z) changes sign as the sterile neutrino density N ν R (z) switches from under to just over its equilibrium value. The parameters here are M 1 = 5 × 10 10 GeV and K = 5, with ε 1 = 10 −6 . drives a positive value for N L (z) before switching, at the cusp in |N L |, to a negative value as N ν R (z) becomes temporarily over abundant. This describes the well-known mechanism of generating the lepton asymmetry from the out-of-equilibrium decays of the sterile neutrinos.
Our main interest here is of course the gravitational contributions, N eq L and W. Consider first the scenario in which N eq L has a hierarchy enhancement, i.e. p = 1 in (6.11). In this case, the dominant contribution to N eq L comes from the heaviest sterile neutrino, so β = 3 in (6.11). A reasonable choice of BSM parameters then gives the evolution shown in Fig. 15, resulting in a final value of |N L | 10 −8 as required to give the observed η. Figure 15. Dynamical evolution of the lepton number asymmetry N L (z) in the case of a hierarchy enhancement, p = 1, with parameters chosen so that the gravitational leptogenesis mechanism dominates over the out-of-equilibrium ν 1 R decays. Here, M 1 = 5 × 10 10 GeV, K = 1, with M 3 = 10 16 GeV, ImK 13 /(4π) 2 = 5 × 10 −4 and ε 1 = 10 −10 .
Here, we have chosen the CP violating decay parameter ε 10 −10 , sufficiently small that the effect of the out-of-equilibrium ν R decays makes a negligible effect on the final lepton asymmetry. N L (z) therefore evolves just as described at the end of section 5. At high temperatures z 1, N L (z) is driven to its equilibrium value N eq L ∼ 1/z 5 , moderated to ∼ 1/z 2 at ultra-high temperatures (with these parameters only for z 10 −6 , which lies outside the range of validity of the effective Lagrangian) by the damping effect of the new evolution term W(z). After a period following its gravitationally-induced equilibrium value, N L (z) decouples as the scattering rate maintaining equilibrium falls below the Hubble expansion rate, i.e. W (z) < 1, and becomes essentially constant. A final dip occurs for z 1 when the resonant increase in W (z) gives a further temporary push towards N eq L , which is falling away very rapidly as 1/z 5 . These features are quantified in the analytic solution in Appendix A.
We therefore see that given the hierarchy enhancement p = 1, the radiativelyinduced gravitational leptogenesis mechanism can reproduce the observed lepton asymmetry for conventional choices of the BSM parameters.
Without a hierarchy enhancement, p = 0, the largest contribution to N eq L arises from the lighter sterile neutrino ν 2 R (that is β = 2 in 6.12), but is still significantly suppressed relative to p = 1. In this case, however, the gravitational terms alone do not give a big enough final value for N L to reproduce the observed η, given the parameter constraints in this particular BSM model from the light neutrino masses in setting the magnitude of W (z). So here we rely on the conventional decay mechanism to produce the final value of N L at low temperatures. Nevertheless, the gravitational effects still produce a radical change in the high temperature evolution of N L (z), illustrated in Fig. 16. Figure 16. Dynamical evolution of the lepton number asymmetry N L (z) in the case of no hierarchy enhancement, p = 0. Here, the final lepton asymmetry is determined by the out-of-equilibrium ν 1 R decays, with the gravitational mechanism producing a steep rise in the asymmetry at earlier times for temperatures z 0.01. Here, M 1 = 10 10 GeV, K = 5, with M 2 = 10 12 GeV, ImK 12 /(4π) 2 = 5 × 10 −4 and ε 1 = 10 −7 .
Following this evolution back in time towards the early universe, the post-leptogenesis value N L 10 −8 arises around z 1 due to out-of-equilibrium decays of sterile neutrinos. As the temperature rises, N L (z) initially falls as in the conventional model. However, this drop is then overtaken by the gravitational effects (from around z ∼ 10 −3 in Fig. 16) and N L (z) rises sharply, driven to its equilibrium value N eq L . Eventually, though outside the range of validity of L ef f in this model, N L (z) evolves away from N eq L due to the new gravitational term W(z) in the Boltzmann equation. Far from reducing to zero for high temperatures above the sterile neutrino mass scale M 1 , the lepton asymmetry drops then rises steeply, its dynamical evolution being driven by the radiatively-induced gravitational tidal effects.
The physics applications of a substantial lepton number asymmetry in the early universe above the conventional leptogenesis scale remain to be explored. In the following section, we generalise the discussion to allow for a post-inflationary era preceding radiation dominance, where the universe expansion is driven by a source with an effective equation of state w = 1/3. As we shall see, this has a very significant effect on the dynamical evolution described here.

Gravitational leptogenesis in the pre-radiation era
This study of the effects of gravity on the dynamical evolution of lepton asymmetry in the radiation-dominated era motivates an exploration of how our gravitational leptogenesis mechanism would be modified in an earlier, post-inflationary era where the universe expansion is controlled by a source with w = 1/3.
An economical way to illustrate how such an equation of state can arise is to consider a single (inflaton) scalar field with different potentials. In an isotropic spacetime, the energy density and pressure are given by For a potential V ∼ φ 2n , it is easy to show [29] that in a (reheating) phase where the scalar field is oscillating around its minimum, the average value of ρ and p over a cycle are related by p = w ρ with w = (n − 1)/(n + 1). For example, a non-interacting scalar field, with a φ 2 potential, has w = 0. The conformal case V (φ) ∼ φ 4 gives w = 1/3, just as for radiation. Steeper potentials give higher values of w, while in the limiting case of 'kination' scenarios, where the kinetic term dominates completely over the potential, the equation of state is w = 1.
In the conventional picture of the early universe, radiation dominance is preceded by a post-inflationary reheating era during which the inflaton field oscillates around the minimum of a shallow potential (such that V ∼ 1 2 m 2 φ 2 ) while decaying into relativistic standard model particles at temperature T . This requires a non-vanishing but very weak coupling of the inflaton to the standard model fields. The process ends when the inflaton decay rate falls below the Hubble parameter and the universe enters the radiation-dominated era at the reheating temperature T rh . During reheating, the effective equation of state parameter therefore evolves in a model-dependent way from an initial value w 0 towards 1/3 as a successively greater fraction of the initial energy density is in the form of radiation. If T rh M 1 , then the gravitational effects described here would be active during reheating.
Gravitational leptogenesis in this era would therefore be characterised by 0 < w < 1/3, a lower value than for the radiation dominance scenario just described. Now, as can be seen from the analysis below, for an equation of state w < 1/3, the Hubble parameter is decreased relative to w = 1/3 resulting in an increase in the factor W (z) = Γ/zH in the Boltzmann equation. Compared to the evolution shown in Fig. 16, this means the lepton asymmetry N L (z) is driven more strongly to the equilibrium value N eq L , reaching it at an earlier time and decoupling later at a higher value of z. This results in a lower final asymmetry for N L due to the gravitational effects prior to the ν R decays. Moreover, this value of N L is further diluted by entropy production during the reheating phase. Overall, therefore, the gravitational leptogenesis mechanism is less effective at generating the required final lepton asymmetry during the reheating phase than it is during radiation dominance. The sharp rise in the asymmetry to N eq L for temperatures above M 1 shown Fig. 16 would however still occur.
Instead, we consider here a scenario in which at early times the universe comprises a non-thermalised 'exotic' component with w > 1/3 which drives the expansion, together with an initially sub-dominant radiation component at temperature T [3,5]. As the universe expands, its energy density naturally decreases faster than the energy density of the radiation, becoming equal at some critical temperature T * = M 1 /z * . Below this temperature, the universe becomes radiation dominated as before.
There are several ways in which such a scenario could arise [5]. An attractive picture [30,31] is that instead of conventional reheating, the relativistic particles giving rise to the entropy of the universe arise through gravitational particle creation at the transition from the de Sitter inflationary vacuum to the FRW vacuum at the end of inflation. This obviates the need for a direct coupling of the inflaton to the standard model fields. Immediately after inflation, the energy density is dominated by the inflaton component with an equation of state w > 1/3, which gradually dilutes relative to the thermalised relativistic particles as the universe expands leaving a conventional radiation-dominated FRW era. If this phase is characterised by a temperature T M 1 then our gravitational leptogenesis mechanism could take place at this time.
The energy density of the radiation component is ρ γ = σT 4 , where T ∼ 1/a(t). From the conservation equation (4.15), the energy density of the exotic component varies with temperature as ρ w = σ w T 3(1+w) , where the constant σ w sets its overall magnitude. We can trade this for the critical temperature T * by setting σ w T 3(1+w) * = σT 4 * . It follows immediately that the ratio of the energy densities of the exotic and radiation components is The total energy density ρ = ρ w + ρ γ .
To see how the various terms in the Boltzmann equation are modified in this new background, we track back through their derivation writing the curvatures and Hubble constant in terms of the new energy density ρ. For the decay rate D(z) and W (z), the only change is due to the factor 1/H in their derivations, so we find where here we define D γ , W γ etc. as the formulae given above for a pure radiationdominated background. For the gravitational terms, there is also a direct dependence on the curvature, and here we find 17) and N eq L (z) = N eq L γ (z) 1 + . (6.18) where w γ = 0.3.
The dependence of W (z) on the parameter w is plotted in Fig. 17 for the two cases where T * is greater or less than M 1 . The key feature is that for z < z * , W (z) becomes smaller as w is increased from its radiation-dominance value 0.3 towards the kination limit w = 1.
We can also show how the relative strength of the evolution factor W(z) and W (z) depends on w. From Fig. 18, we see how the crossover point where W(z) W (z) occurs at successively lower temperatures as w is increased.
These features explain the evolution of the lepton asymmetry N L (z) shown in Fig. 19. The qualitative features are the same as described at the end of section 5, and exact analytic expressions for N L (z) in the different regions are derived for arbitrary w in Appendix A. At ultra-high temperatures the evolution term W(z) dominates  in the Boltzmann equation and gives the universal behaviour N L (z) ∼ 1/z 2 for all w. Beyond the crossover point in Fig, 18, N L (z) is driven to the equilibrium value N eq L , which in itself increases with w, by the lepton number violating rate term W (z). However, since W (z) is weaker for successively bigger values of w, the asymmetry N L (z) decouples sooner from N eq L and consequently at a higher value. Indeed, as we approach the kination limit, N L (z) is so weakly driven to N eq L that after the initial ultra-high temperature phase when W(z) > W (z), it remains essentially constant and the whole Boltzmann equation is dominated by the new evolution term W (see Appendix A). Overall, the point of decoupling and the final asymmetry are clearly very sensitive to the value of w.
We see, therefore, that for the same BSM parameters, the final asymmetry produced by the gravitational leptogenesis mechanism increases with the equation of state parameter w. An optimal case is N L (z) with w = 0.5 in Fig. 19, yielding a final asymmetry N L 10 −8 . This shows how with physically reasonable BSM parameters, the observed baryon asymmetry η can be entirely generated by radiatively-induced gravitational leptogenesis, even without a hierarchy enhancement, in a FRW spacetime characterised by an effective equation of state w 0.5. This could readily be realised in the post-inflationary scenario of [30,31].

Summary and Outlook
In this final section, we highlight some of the key features of radiatively induced gravitational leptogenesis and discuss potential future developments.
First, we again emphasise the generality of the leptogenesis mechanism, independent of its realisation in the particular BSM theory analysed here. Whatever its origin, the effective action (3.1) summarises the gravitational interactions of the light neutrinos and implies the picture of the origin and evolution of the lepton number density described above. We saw there was a clear distinction between the roles of the CP odd and CP even operators in (3.1). The coupling ∂ µ R ν L γ µ ν L to the CP odd neutrino operator modifies the light neutrino dispersion relation and gives rise to an effective chemical potential for lepton number. This is the origin of the non-vanishing equilibrium lepton number density n eq L . On the other hand, the gravitational couplings to the CP even operators, especially R µν ν L γ µ ← → D ν ν L , modify the dynamical evolution of the lepton number density at early times.
To complete the gravitational leptogenesis mechanism, these two tidal curvature effects encoded in the effective action need to be augmented by a lepton number violating reaction to drive the lepton number density n L towards its equilibrium value n eq L and maintain it at this value until decoupling. These non-gravitational ∆L = 0 reactions are required to satisfy the first Sakharov condition and are provided naturally here by the fundamental BSM theory.
Together, these effects give rise to the generalised Boltzmann equation (5.25) which implies the picture of the dynamical evolution of the lepton asymmetry summarised in Fig. 11. An initially zero lepton-to-photon ratio N L is driven rapidly towards N eq L , its approach moderated by the new factor W(z) in the Boltzmann equation. N L then tracks N eq L until the ∆L = 0 reaction rate drops towards the Hubble expansion rate and can no longer maintain N L in quasi-equilibrium. At this point, N L decouples giving (up to the final dip in region (d) of Fig. 11) the final gravitationally-induced lepton asymmetry.
As we have seen, the final value for the lepton-to-photon ratio N L depends very sensitively on the temperature at which decoupling from equilibrium occurs. This is because N eq L (z) is falling very sharply with temperature, as 1/z 5 for the radiationdominated background. This behaviour is inherited from the original curvature dependence of n eq L in (5.15), from which n eq L ∼ ρ H T 2 . In turn, the point of decoupling is determined by W (z), so depends on the Hubble parameter (as exploited in the w > 1/3 backgrounds described in section 6.2) but also on the strength of the ∆L = 0 reactions, which is dependent on the parameters of the BSM theory.
Turning now to the features specific to the see-saw BSM theory we have used here to illustrate RIGL, note that it has three closely-linked properties relevant to the analysis of leptogenesis. First, it exhibits the lepton number violating interactions in Fig. 1, which depend on the exchange of virtual ν R neutrino propagators with S × propagators. Reorienting these diagrams and replacing the Higgs fields by their VEV gives the seesaw mass generation mechanism for the light neutrinos. This explains the parameter constraint linking the magnitude of W (z) to the observed neutrino masses. Third, the same interaction vertex necessarily allows the ν R → ν L H and ν R → ν c L H decays of real sterile neutrinos. At one-loop level, these decays exhibit CP violation as described in Fig. 12. As they can occur when the ν R are out of equilibrium around T ∼ M 1 (see Fig. 14), this realises the Fukugita-Yanagida leptogenesis mechanism [12]. We noted that the dynamical weighting of the CP-violating combination of Yukawa couplings at O(λ 4 ) governing these decays is similar, but not identical, to the combination arising in the gravitationally-induced lepton number density n eq L . This difference allows us to choose parameters such that one or other mechanism for leptogenesis dominates the final asymmetry. Both mechanisms are however necessarily present in the seesaw model, or any related BSM theory with L-violation occurring through Feynman diagrams of the type in Fig. 1.
It would be interesting to find alternative BSM theories in which this tight link between ∆L = 0 reactions, out-of-equilibrium ν R decays, and neutrino masses could be relaxed, giving a greater flexibility in the parameter choices controlling the magnitude and evolution of the lepton asymmetry in RIGL.
Further exploration of the evolution of the lepton number density n L at very early times with T M 1 requires going beyond the 'low-energy' approximation used in the effective action. Recall that this is constructed only to leading order in derivatives, and implies the restriction T √ R/M 2 1 < 1 on its range of validity. To complete the description of the evolution of n L beyond this region, we need more powerful techniques to evaluate the one and two-loop self-energy diagrams in section 2 beyond the smallmomentum expansion.
At one loop, the required methods have been developed in the series of papers [16][17][18] where the high-energy behaviour of photon vacuum polarisation diagrams in QED has been analysed. The motivation in these papers was to examine fundamental issues involving causality, analyticity and unitarity in quantum field theories in curved space-time. The same methods, exploiting the Penrose limit of the background spacetime, may readily be adapted to the one-loop self-energy diagrams here, using the generalisation to fermions described in [18], so in principle the evolution term W(z) can be found for arbitrarily high temperatures, provided only that R/M 2 1 < 1. The remaining technical difficulty is in extending these methods to the two-loop diagrams necessary to find n eq L . Two loops seems to be necessary since, at least in the BSM theory considered here, CP violation requires going to fourth order in the Yukawa couplings. Clearly, the ∼ 1/z 5 rise in N eq L will not continue indefinitely for very small z, and the evidence from QED vacuum polarisation strongly suggests that all these gravitational loop effects revert to zero in the ultra-high energy limit. This was essential in reconciling causality with the apparent superluminal propagation in the low-energy effective action for QED in curved spacetime. It would be very interesting either to overcome this difficulty (which seems unlikely given the complexity in evaluating all the required two-loop diagrams even in the low-energy approximation) or to find an alternative BSM theory in which n eq L is generated at one loop so that the entire dynamical evolution of n L (z) from the earliest times could be quantitatively traced.
In conclusion, in this paper we have extended the theory of radiatively-induced gravitational leptogenesis developed in [1][2][3]7] to give a fuller picture of the dynamical evolution of the lepton number asymmetry in the early universe. Whether the final value of the lepton asymmetry leading to the observed baryon-to-photon ratio η is determined by the conventional out-of-equilibrium sterile neutrino decay mechanism, or entirely through the generation of the equilibrium asymmetry n eq L by tidal curvature effects at the quantum loop level, we have seen how at earlier times the universe will experience a phase with a significantly higher lepton-antilepton asymmetry determined purely by gravitational effects, as illustrated in Fig. 11. It remains an open and interesting question what other physical consequnces may follow from the existence of such a non-vanishing matter-antimatter asymmetry in these very early times in the evolution of the universe.

A Analytic results for the Boltzmann equation
The Boltzmann equation (5.25) for gravitational leptogenesis admits simple analytic solutions in the temperature regions (a), (b), (c) in Fig. 11. This provides additional insight into the numerical plots for the lepton number asymmetry N L (z) presented in the text.
Keeping only the terms relevant to the gravitational mechanism, i.e. neglecting the contribution from ν R decays, the Boltzmann equation is In the region z < 1 (and also z < z * in the cosmological scenario in section 6.2), this is well-approximated by simple power-law behaviours of the coefficients and we can write, Consider first the radiation-dominated era. Here, p = 2, q = 5, r = 5 and the coefficients α, β and γ can be read off from (6.9), (6.11) and (6.13). It is straightforward to show that the general solution of (A.2) is  16 It is clear that with this explicit solution, N L (z) rises almost instantly towards its value in region (a) of Fig. 11, driven by the large initial value of the W N eq L term in the Boltzmann equation for N L (z). Note, however, that this is an artifact of using the effective Lagrangian for the Boltzmann coefficients at extremely early times outside its range of validity. In reality, we expect from previous experience with QFT loop calculations in curved spacetime that for sufficiently high temperatures, where T √ R/M 2 1 1, the loop effects will revert smoothly to zero, so the transition of N L (z) from its initial zero value at z = z 0 towards its value in region (a) is expected to be gradual. For small z, i.e. such that γ/z 4 1, we can approximate this using the asymptotic expansion erfc(x) e (−x 2 ) (1/ √ πx) (1 + O(1/x 2 )) at large x, so we find This explains the slower rise of the lepton asymmetry in region (a) compared to the sharp 1/z 5 rise of the equilibrium value N eq L .
This matches smoothly to N L (z) ∼ β/z 5 in region (b) where N L (z) follows its equilibrium value. From (A.7), the crossover occurs at z 3 γ/α, which of course matches the point where W (z) W(z). This matches the smooth evolution of N L (z) from region (b) through decoupling and into region (c). The dominant term in region (c) is given by the n = 0 term in the sum, leaving [3] 17 N L (z) 120 β α 5 , (z > α) . (A.10) The same methods can be applied to the pre-radiation cosmological scenario in section 6.2 using the expressions (6.16), (6.18) and (6.17) for W (z), N eq L (z) and W(z) approximated for z z * . In this case, the Boltzmann equation (A.2) has the powerlaw dependence p = 1 2 (5 − 3w), q = 1 2 (7 + 9w) and r = 4 + 3w. The coefficients are w = 1, N L (z) never holds to its equilibrium value and its final constant value is given by the large z expansion of (A.12), viz.
For intermediate values of w, we need to keep the full Boltzmann equation (A.2) balancing both the W (z) and W(z) terms, with the numerical solutions for the lepton number asymmetry shown in Fig. 19.