On different approaches to freeze-in and freeze-out leptogenesis with quasi-degenerate neutrinos

We compare two approaches for determining the generation of lepton asymmetry during production and decay of quasi-degenerate neutrinos, namely the density matrix formalism and a recent proposal which does not involve any counting of neutrino number densities and is based on plugging the resummed propagator in a quantum field theory model for neutrino oscillations. We show numerically and analytically that they are almost equivalent for small mass splittings and also discuss the hierarchical limit. The comparison, performed in a simple scalar toy model, helps to understand several issues that have been discussed in the literature.


Introduction
Quasi-degenerate exotic neutrinos allow for low-scale leptogenesis mechanisms that may be tested in current and near future experiments.More specifically, in the type I seesaw model the baryon asymmetry of the universe can be generated in the freeze-out of Majorana neutrinos with typical masses around O 10 2 − 10 3 GeV , and also in the freezein of much lighter neutrinos, with masses even below the GeV scale.The former mechanism is commonly known as resonant leptogenesis [1] and the latter as baryogenesis via neutrino oscillations (or ARS leptogenesis) [2,3].Both mechanisms have been analyzed with different formalisms, see e.g.[4], [5], and section II of [6] for comprehensive reviews and references.
Indeed the treatment of CP violation in baryogenesis models with quasi-degenerate neutrinos and its implementation in transport equations is not trivial and has actually been discussed over several decades.Some points of concern include the proper calculation of the CP asymmetry in neutrino decays to avoid a divergence in the limit of equal masses (and particularly in the double degenerate limit of equal masses and couplings), the distinction and interplay of two sources of CP violation (mixing and oscillations), and the joint analysis of freeze-in and freeze-out leptogenesis with a single set of kinetic equations valid for any spectrum of the heavy neutrino masses (cf.[6][7][8][9][10][11][12][13][14][15][16][17][18][19]).As a matter of fact, there are differences in the literature regarding these issues, like in the regulators for the decay asymmetry (reviewed e.g. in [4]), in the relative sign of the contributions from mixing and oscillations (compare [11][12][13][14] with [16][17][18]), in the existence of an interference term between mixing and oscillations (compare [11][12][13] with [14,16,18]), and in the set of kinetic equations (cf.[11,15]).Note in particular that for a hierarchical spectrum of neutrino masses, the standard classical approach requires the inclusion of real intermediate state (RIS) subtracted rates to be consistent with unitarity, which are present in [11] but not in [15] (although the correct hierarchical limit is claimed to be obtained in both works, see also [16][17][18]).
Motivated by these considerations we compare in this work two different approaches to freeze-in and freeze-out leptogenesis, namely the one based on quantum kinetic equations for matrices of densities and the one proposed in [16].The density matrix approach (DMA), like the one used in [6,15] to study jointly freeze-in and freeze-out leptogenesis, is based on a generalization of the Sigl and Raffelt formalism [20] to include heavy neutrinos, and can also be derived under some approximations from non-equilibrium quantum field theory (see e.g.[8,10,17,21]).In the approach proposed in [16], an expansion around the complex poles of the resummed one-loop propagator is plugged into a quantum field theory model of neutrino oscillations, in order to obtain time dependent probabilities for lepton number violating processes.The CP asymmetry obtained from these probabilities is suitably integrated over time in order to get a source term for the evolution of the lepton asymmetry.Given that in [16,18] we have employed an external wave packet model for neutrino oscillations [22,23] (following [24]), we will refer to it as the external wave packet approach (EWPA).The EWPA is particularly transparent regarding unitarity (which is a key principle to guide understanding), because no counting of the heavy-unstable-neutrino states is performed.
The issues we are interested in discussing can be captured in a simple scalar toy model, consequently we organize this work as follows: In section 2 we review the main results from [16,18] and work out the kinetic equations in the DMA for the scalar toy model, then in section 3 we compare both approaches and discuss the results, finally we conclude in section 4.

Implementation of two approaches in a scalar toy model
We will consider a simple scalar toy model which has often been used to study basic aspects of leptogenesis with quasi-degenerate neutrinos.The particle content is given by one complex and two real scalar fields, denoted by b and ψ i (i = 1, 2), respectively.The Lagrangian, in the basis where the mass matrix of the real scalars is diagonal, can be written as The b-particles play in this toy model the analogous role that leptons play in standard leptogenesis, hence they will be referred to as "leptons", whose mass m will be neglected for simplicity.Lepton number is violated by the cubic Yukawa interaction terms involving the ψ i , to be called "neutrinos" henceforth.The last term is a quartic interaction which does not violate lepton number, therefore it will not appear explicitly in our analysis, although we note that it could be used as a way to localize the leptons and thus satisfy the conditions to have oscillations [24].
Next we describe two approaches that can be implemented to study the production of lepton symmetry in this model.

External wave packet approach
The EWPA, proposed in [16] and extended to highly degenerate states in [18], proceeds in three steps: (1) expansion around the complex poles of the resummed propagator, (2) calculation of a time dependent CP asymmetry between lepton number violating processes using a quantum field theory model for neutrino oscillations, and (3) integration over time of this CP asymmetry to obtain a source term for the evolution of the lepton asymmetry.Below we extract the main formulae from [18] that will be needed in this work.
The resummed one-loop propagator matrix, G, is given by iG −1 (p 2 ) = p 2 1 − M 2 (p 2 ), with and M 2 a,b the poles of the propagator given by the roots of the determinant of G −1 .At O h 2 (with h representing any of the Yukawa couplings), the elements of the Z matrix are equal to , and the O h 4 terms are finite in the limit → 0. Taking The last equalities in each of the above equations define the real quantities M1 , M2 , Γ1 and Γ2 .
Using the approximate propagator in an external wave packet model of neutrino oscillations and neglecting the factors related to coherence and localization (which could destroy oscillations), one arrives at these time-dependent probabilities for the lepton number violating processes (see [18]): where , N is a normalization constant, t is the time between production and decay of the neutrinos that mediate these processes, γ ≡ E/M is the Lorentz factor, E is the average energy of the neutrinos, and Finally, a proper integration over time of the CP asymmetry |A|2 − Ā 2 yields a source term, S EWPA , for the evolution of the lepton density asymmetry n L ≡ n b −nb, where n b (nb) is the number density of leptons (antileptons).For our purposes it is sufficient to consider a static universe, a single average momentum p for the neutrinos (to avoid unessential integrals over momentum), and neglect finite density effects.Hence the normalization constant N in eq.(2.6) is the same as in [18] The source, S EWPA (t), contains the terms in the kinetic equation which may be non-null in the absence of a lepton density asymmetry and will be the focus of this work, while the washout part, W (t), will not be considered.Moreover, n eq (t) is the equilibrium number density of a scalar particle with mass M .Note that in an expanding universe, n eq (t) varies according to the time-dependent temperature, while for our simplified study in a static universe n eq will be varied artificially, with equilibrium corresponding to constancy over time.

Density matrix approach
Following [20] we can obtain a kinetic equation for the density matrix ρ of the neutrinos and write it as where each term depends on the momentum k of the neutrinos, Γ d is a decay rate matrix, γ p is the production term and we will work in the basis where the Hamiltonian is diagonal, . Note that for simplicity we have used Maxwell-Boltzmann statistics and we are not including terms proportional to the lepton asymmetry.In the scalar model the decay rate matrix is given by with h = (h 1 h 2 ) the matrix of Yukawa couplings.Note that decay into both, leptons and antileptons (with total energy E), have been included in this expression.As explained e.g. in [3], the production and decay terms should be related in order for the right-hand side of eq.(2.8) to vanish in equilibrium, namely with ρ eq = n eq (t) I and I the identity matrix, so that the kinetic equation for ρ becomes Then the equation for the lepton density asymmetry, neglecting washouts, can be obtained taking into account that for each production (destruction) of a neutrino state two leptons or antileptons are destroyed (produced): where Tr (.) denotes the trace of the corresponding matrix.Given that (h † h) T = h T h * , the first term, which comes from the destruction of leptons and antileptons, is actually zero, therefore only the production of lepton asymmetry from the decay of the neutrinos contributes to this equation: (2.12)

Comparison and discussion
To perform an analytical comparison between the EWPA and the DMA, we start by noticing that from Then, expanding M 2 (p 2 ) around M 2 1,2 we get This is the square of the relation we are interested in, namely (using that Z Z T = I +O h 4 ) where we have introduced the effective Hamiltonian, This relation can be used to obtain an approximate expression for the probabilities of lepton number violating process given in eq.(2.6) that are used in the EWPA (again taking into account that Z Z T = I + O h 4 ): (3.3) In the same way, for the CP-conjugate process, On the other hand, the formal solution to the eq.(2.11) of the DMA can be written as with ρ p ≡ h † h + h T h * .Replacing this expression in the source term of the DMA (i.e. in eq.(2.12)), some terms cancel and the remaining ones read where we have used the cyclic property of the trace.Therefore, comparing this last equation with eqs.(2.7), (3.3) and (3.4), we conclude that The solid red lines correspond to the EWPA approach, the dashed black lines to the DMA, the blue line in the top right plot to the mixing contribution of the EWPA, and the dotted green lines in the bottom plots to the standard classical approach valid in the hierarchical limit.In the top plots we have chosen Γ 1 /M 1 = 1/100, Γ 2 /M 1 = 1/120 and ∆M = Γ 1 , while for the bottom plots, Γ 1 /M 1 = 1/100, Γ 2 /M 1 = 1/120 and ∆M/M 1 = 4.In all cases we have taken h 1 = |h 1 | and h 2 = |h 2 | e iφ , with the phase φ = π/4, and n eq (t) = e −M 1 t/(10000 γ) .The lepton asymmetry has been obtained integrating only the source term (i.e.ignoring washouts) and the scale on the vertical axis of the plots is not relevant (note that after a change of variables in the integration over time, the factor M 1 /γ becomes part of the normalization chosen for the lepton asymmetry).
i.e. S DMA and S EWPA are equal up to higher order terms in h 2 and/or .
To corroborate and exemplify this result, we plot in figure 1 the evolution of S DMA and S EWPA (left plots), as well as the corresponding lepton asymmetries (right plots), for two different values of the mass splitting ∆M ≡ M 2 − M 1 .The lepton asymmetries have been obtained by integrating the corresponding source term over time (ignoring washouts) and in both cases we have taken n eq (t) = e −M 1 t/(10000 γ) , in order that the time scale of the evolution of n eq be much larger than the lifetimes of the neutrinos.
The top plots illustrate the highly degenerate case by taking ∆M = Γ 1 .There are two main things to notice.First, the left plot shows that S DMA and S EWPA are-almostequal at all times.Second, from the right plot it can be seen that there actually small differences between both sources, coming from the O h 8 , h 6 , h 42 terms in eq.(3.7), which accumulate over time and result in noticeable differences between the lepton asymmetries at late times.In particular, while the final lepton asymmetry in the DMA is null, it is not so in the EWPA.This is interesting because unitarity and CPT invariance imply that when washouts are negligible and the initial abundances of neutrinos and lepton asymmetry are null, the final asymmetry must also be null.Therefore the discrepancy between the final lepton asymmetries reflect the non-trivial differences in the-crucial-implementation of unitarity constraints, as explained next.
On one hand, in the kinetic equations of the DMA, as derived above or e.g. in [3,20], the right-hand side of eq.(2.8) is made zero in equilibrium by relating the production and decay rates as in eq.(2.10).Moreover, the equation for the lepton asymmetry is inferred from some particle number conservation condition.It is then easy to show analytically that the final lepton asymmetry in the DMA is exactly null.Namely, the final lepton asymmetry in the DMA, n f L , can be obtained integrating the source S DMA given in eq. ( 2.12), where we have taken n L (t = 0) = 0 and for convenience the null contribution from the terms proportional to ρ eq has been included.Now, from eq. ( 2.11) and using that [H, ρ] = [H, ρ − ρ eq ], it is clear that ∞ 0 (ρ − ρ eq ) dt = 0 when the initial and final densities of neutrinos are null, and therefore that n f L is zero in the DMA.This procedure to obtain the kinetic equations in the DMA has sometimes been justified by detailed balance conditions.However, in a time non-invariant theory, as is the case when CPT is conserved but CP is violated, detailed balance may not hold and the canonical equilibrium distributions result from a cyclic balance involving several processes, which is precisely ensured by unitarity (see e.g.[26,27] and note that this is at the root of the RIS subtraction procedure to be discussed below).Moreover, this procedure does not clarify whether the rate in eq.(2.9) should be calculated at tree level (which leads to an agreement with the EWPA) or involve effective couplings to account for CP violation in production and decay (due to the presence of oscillatory and secular terms, as well as enhancements inversely proportional to the mass splitting, higher order corrections to the rates could be non-negligible) 2 .
On the other hand, in the EWPA there is no counting of neutrino abundances, avoiding unitarity problems that may arise when using perturbation theory with unstable particles in the initial or final states.Instead, the source term is built from the probabilities in the eqs.(2.6), which only involve the stable leptons and antileptons as asymptotic states, and have been obtained after a perturbative expansion of the propagator.Then unitarity constraints are satisfied as a matter of course up to higher order terms in the expansion.Notably, as explained in [16,18], there is a non-trivial cancellation between three different terms contributing to the final lepton asymmetry.These terms (dubbed mixing, oscillation, and interference) arise when calculating the CP asymmetry |A| 2 − Ā 2 from the eqs.(2.6) and will be discussed below in greater detail.To illustrate the cancellation we have plotted one of these contributions in figure 1 (the mixing term represented by the blue curve in the top right panel), that should be compared with the net lepton asymmetry (red curve), whose final value is almost null-but not exactly null due to higher order corrections.
In the bottom plots we have taken a very large mass splittings, ∆M/M 1 = 4, to illustrate the hierarchical limit.As can be seen, there are differences of around 50% in the lepton asymmetries obtained from the DMA and EWPA.We have verified numerically that the main difference between the left-hand side of eq.(3.1) and H ef lies in the imaginary part of the diagonal elements.Overall, however, the conclusion is that differences are not very significant and only appear at very large mass splittings for which, actually, decoherence effects not included in either approach may become dominant and destroy oscillations (see e.g.[24]).That being said, it is interesting to note that the lepton asymmetry from the EWPA exactly matches the one obtained with the standard classical Boltzmann equations (green curve in figure 1 ), except at early times (of the order of the oscillation period), and late times for the same reason as explained above, namely that the classical Boltzmann equations exactly satisfy the unitarity condition of null final asymmetry.As is well known, in order for the latter to happen it is key to include the rates of some processes subtracting RIS contributions.Therefore it is reasonable to ask whether this type of RIS-subtraction procedure should also be performed somehow in the DMA in order to be valid in the hierarchical limit (cf.[11]).Indeed, as we have shown, such procedure is not required since both, the DMA and the EWPA, which are equivalent up to higher order terms (eq.(3.7)), basically yield the same hierarchical limit as the classical Boltzmann equations (disregarding the minor differences mentioned above).
In [16,18] we identified two different types of CP even phases in eqs.(2.6), one independent of time in θ , and an oscillating one in the exponentials e −i Mj t/γ .Consequently the source term could be written as a sum of three contributions, one from mixing (involving only θ ), one from oscillations (involving only e −i Mj t/γ ), and interference terms (involving both θ and e −i Mj t/γ ).This result went along the lines of previous findings [11][12][13][14], although differences remained regarding the relative sign between the mixing and oscillation contributions and the existence or not of an interference term.However, in the DMA no such separation into CP violation from mixing and oscillations appears (see e.g.[6,15]).Now this can be understood from eqs. (3.3) and (3.6): if the effective Hamiltonian H ef present in the source term of the DMA (eq.(3.6)) is diagonalized by an orthogonal matrix, the time exponential involving H ef can be diagonalized by the same matrix, leading to an expression like the first line of eq.(3.3),where it is possible to identify a time independent CP even phase in the orthogonal matrix and oscillating CP even phases associated to the eigenvalues of H ef .Indeed, in the hierarchical limit this decomposition into time independent and oscillating CP even phases becomes quite meaningful, since the oscillation contribution can be identified with the contribution from RIS subtracted rates to the source of the classical Boltzmann equations, as shown in [16], while the mixing contribution matches the contribution from decays to the classical source (see also, e.g., [6,29] for related limits).
Moreover, the EWPA is valid for relativistic as well as non-relativistic neutrinos and does not involve any counting of neutrino states, therefore it can clearly be applied whether the lepton asymmetry is mostly generated during production or decay.Therefore we think the equivalence we have demonstrated between this approach and the DMA, provides further support to the validity of the DMA for describing jointly freeze-in an freeze-out leptogenesis, although issues related to the helicity degree of freedom cannot be captured in the scalar toy model (see [6,15,17,[30][31][32][33][34][35][36][37]).Finally, note that in the EWPA, as in other approaches, it is quite involved to obtain consistent results in the degenerate limit, and even more in the double degenerate limit of equal masses and couplings.Therefore the agreement between the EWPA and the DMA, which does not suffer from this problem, supports the treatment of the degenerate limit in [18].

Conclusions
We have shown analytically that the EWPA and the DMA are equivalent up to higher order corrections in the couplings and mass splitting (see eq. (3.7)).The fact that there are differences coming from higher order terms, actually makes the comparison between both approaches more meaningful, and in particular they reflect the different implementation of unitarity requirements, as discussed in section 3.For large mass splittings the lepton asymmetry obtained with the EWPA matches the one from the standard classical Boltzmann equations (except at very early and late times, as explained in section 3), while differences of some tens of percent arise with the DMA.Overall, however, we noted that these differences are not very significant and only appear at very large mass splittings (∆M/M 1 = 4 in figure 1), when anyway decoherence effects not included in any of the approaches could become dominant.
The connection between both approaches established in this work, helps to understand several issues that have been discussed quite extensively in the literature.In particular, we have shown that the production and decay rates in the DMA should be calculated at tree level (i.e.without involving effective couplings to account for CP violation in production and decay), in order to get the agreement given by eq.(3.7) with the EWPA, which does not bring this type of doubts.Regarding the meaningfulness or existence of different types of contributions to the generation of lepton asymmetry associated to CP violation in mixing, oscillations and interference, we have shown how to get this decomposition from a solution to the kinetic equations of the DMA.In turn, given that in the hierarchical limit the oscillation contribution can be identified with the contribution from RIS subtracted rates to the source of the classical Boltzmann equations (see [16]), this provides further insight into how the hierarchical limit can be obtained from the DMA without including RIS subtracted rates.Moreover, the validity of the EWPA in the degenerate limit and the double degenerate limit of equal masses and couplings is confirmed by its concordance with the DMA.Finally, the agreement between both approaches also supports their use for a joint analysis of freeze-in and freeze-out leptogenesis, although issues related to the helicity degree of freedom cannot be captured in the scalar toy model we have employed and might be interesting to address in future works.Another possibility worth considering for future research along the lines of this work is the tri-resonant heavy neutrino system [38], consisting of three nearly degenerate neutrinos.

Figure 1 :
Figure1: Absolute value of the source term (left plots) and lepton asymmetry (right plots), as a function of time normalized to γ/M 1 .The solid red lines correspond to the EWPA approach, the dashed black lines to the DMA, the blue line in the top right plot to the mixing contribution of the EWPA, and the dotted green lines in the bottom plots to the standard classical approach valid in the hierarchical limit.In the top plots we have chosen Γ 1 /M 1 = 1/100, Γ 2 /M 1 = 1/120 and ∆M = Γ 1 , while for the bottom plots, Γ 1 /M 1 = 1/100, Γ 2 /M 1 = 1/120 and ∆M/M 1 = 4.In all cases we have taken h 1 = |h 1 | and h 2 = |h 2 | e iφ , with the phase φ = π/4, and n eq (t) = e −M 1 t/(10000 γ) .The lepton asymmetry has been obtained integrating only the source term (i.e.ignoring washouts) and the scale on the vertical axis of the plots is not relevant (note that after a change of variables in the integration over time, the factor M 1 /γ becomes part of the normalization chosen for the lepton asymmetry).