Leptogenesis via Axion Oscillations after Inflation
IIPMU 14-0352
Leptogenesis via Axion Oscillations after Inflation
Alexander Kusenko,
1, 2
Kai Schmitz, ∗ and Tsutomu T. Yanagida Department of Physics and Astronomy, University of California, Los Angeles, California 90095-1547, USA Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan (Received 18 December 2014; revised manuscript received 18 February 2015; published 2 July 2015)Once a light axionlike scalar field couples to the electroweak gauge bosons, its classical motionduring reheating induces an effective chemical potential for the fermion number. In the presenceof rapid lepton number ( L )-violating processes in the plasma, such a chemical potential provides afavorable opportunity for baryogenesis via leptogenesis. We are able to demonstrate that L violationdue to the exchange of heavy Majorana neutrinos is sufficient for a successful realization of this idea.Our mechanism represents a novel and minimal alternative to thermal leptogenesis, which turns outto be insensitive to the masses and CP -violating phases in the heavy neutrino sector. It is consistentwith heavy neutrino masses close to the scale of grand unification and, quite complementary tothermal leptogenesis, requires the reheating temperature to be at least of order 10 GeV.
PACS numbers: 98.80.Cq, 11.30.Fs, 14.60.St, 14.80.Va
The Peccei-Quinn (PQ) solution to the strong CP problem [1] has led to the prediction of a light scalarfield called the axion [2], which appears generically in anumber of models [3, 4] and which has a characteristiccoupling to the gauge fields of the form f − a a F ˜ F , where f a is the scale of PQ symmetry breaking. However, themotivation for considering axionic fields extends well be-yond the context of the strong CP problem. Axions areubiquitous in string theory, where at least one such fieldis generically associated with the Green-Schwarz mecha-nism of anomaly cancellation [5] and the scale f a of whichlies a few orders of magnitude below the Planck scale [6];but multiple other axions can also appear. While themodel-independent axion is coupled to all gauge groupswith a universal coupling strength, the additional axionicfields can couple to different groups with couplings thatdepend on both the gauge group and the particle contentof the model [7]. The masses of these model-dependentaxions can be different as they arise from their couplingsto different anomalous groups. We will focus, in partic-ular, on the axion (or linear combination of axions) thathas a coupling to the electroweak SU (2) gauge fields.During inflation, light scalar fields develop large expec-tation values [8]. The relaxation of the axion field to theminimum of its effective potential begins once the Hub-ble rate becomes comparable to the axion mass. Whilethe field a ( t ) is slowly rolling towards its origin, its cou-pling to the SU (2) gauge fields and, via the anomaly,to the fermionic current j µ = ¯ ψγ µ ψ induces an effective CP T -violating term a ( t ) F ˜ F ∝ ( ∂ t a ( t )) j , which servesas a chemical potential for fields carrying nonzero baryonor lepton number. Then, in the presence of rapid lep-ton number-violating processes in the plasma—for exam- This illustrate that, while the axion could equally couple to thehyercharge gauge boson, a coupling to the standard model gluonswould, by contrast, not provide a sufficient basis for leptogenesis. ple, due to the exchange of virtual heavy right-handedneutrinos—the conditions for successful leptogenesis aresatisfied. Here, one important detail is that lepton num-ber violation needs to occur at a fast rate, Γ L (cid:29) ˙ a/a ,so that the axion field acts as an adiabatic backgroundduring leptogenesis. It is therefore essential that the gen-eration of the lepton asymmetry be driven by an externalsource, i.e., by scatterings with heavy neutrinos in the in-termediate state in our case. If the rate of lepton numberviolation was instead tied to the change in the axion fieldvalue, we would, by contrast, not be able to interpret theaxion velocity as an effective chemical potential [9]. Aswe will see, a consequence of the requirement of a fastrate Γ L is, in particular, that we will have to adopt largevalues for the reheating temperature, T rh (cid:38) GeV.A similar scenario was discussed in connection witha flat direction that carries no baryon or lepton num-ber [10]. Our treatment of the asymmetry is different,and we obtain very different results. Our scenario is alsosimilar to leptogenesis via Higgs relaxation [11], wherethe Higgs coupling to F ˜ F is assumed to arise from ahigher-dimensional operator, unlike in the present sce-nario, where the required coupling appears genericallyfor any axion coupled to the electroweak gauge fields.One can also draw an analogy to models of spontaneousbaryogenesis [12], in particular to realizations of sponta-neous baryogenesis during the electroweak phase transi-tion (EWPT) [13]. Here, the Higgs field in the expandingbubble wall generates an effective chemical potential forthe fermions. Our scenario is different in that the “wall”is represented by an axionic field moving in the timelikedirection uniformly in space, unlike the bubble wall mov-ing in a spacelike direction. Finally, we note that, in thecontext of electroweak baryogenesis, also the QCD axionmay be employed to generate an effective chemical po-tential [14]. Successful baryogenesis then requires thatthe EWPT be delayed below the GeV scale; for instance,due to some additional Higgs-dilaton coupling [15]. a r X i v : . [ h e p - ph ] J u l Provided that the axion a is to be identified with thepseudo-Nambu-Goldstone boson of a spontaneously bro-ken U (1) symmetry with a compact global topology, itsinitial value at the end of inflation is a = f a θ , where theangle θ takes a random value in the range θ ∈ [0 , π ).Assuming that the PQ symmetry is broken sufficientlyearly before the end of inflation (and is not restored dur-ing reheating), this initial value ends up being constanton superhorizon scales. For definiteness, we set a = f a and treat f a as a free parameter in the following. Antic-ipating that the final baryon asymmetry will depend on a , we require that the baryonic isocurvature perturba-tions induced by the quantum fluctuations of the axionfield during inflation be smaller than the observationalupper limit. This implies a constraint on the Hubblerate during inflation: H inf / (2 π ) /a (cid:46) − [16], or H inf (cid:46) × GeV (cid:18) f a GeV (cid:19) . (1)The evolution of the homogeneous axion field in itseffective potential V eff around the origin is described by¨ a + 3 H ˙ a = − ∂ a V eff , V eff ≈ m a a , (2)where we have neglected the backreaction of lepton num-ber generation on the evolution of the axion field. m a denotes the axion mass, which we assume to arise viadimensional transmutation, i.e., from an additional cou-pling of the axion to the gauge fields of some stronglycoupled hidden sector. Given a dynamical scale Λ H inthis hidden sector, the axion mass is then of O (cid:0) Λ H /f a (cid:1) .For consistency, we require m a to be smaller than H inf ,the Hubble rate at the end of inflation: m a (cid:46) H inf . (3)When inflation is over, the axion field remains practicallyat rest until the Hubble parameter drops to H osc = m a .Once the axion field is in motion, the effective Lagrangiancontains the term L eff ⊃ g π a ( t ) f a F ˜ F = − a ( t ) N f f a ∂ µ (cid:0) ¯ ψγ µ ψ (cid:1) (4)= ∂ t a ( t ) N f f a (cid:0) ¯ ψγ ψ (cid:1) + · · · = µ eff j + · · · , (5)with g being the SU (2) gauge coupling and N f = 3 thenumber of fermion generations in the standard model,where we have used the anomaly equation in Eq. (4), andintegration by parts in Eq. (5). In the following, we willabsorb N f in our definition of f a and simply determinethe effective chemical potential as µ eff = ˙ a/f a .Now the necessary conditions for generating a leptonasymmetry are satisfied. A nonzero effective chemicalpotential shifts the energy levels of particles as comparedto antiparticles. If lepton number is not conserved, the minimum of the free energy in the plasma is reached for adifferent number density of leptons than for antileptons,i.e., for n L ≡ n (cid:96) − n ¯ (cid:96) (cid:54) = 0. Instead, if the lepton numberviolation is very rapid, the minimum of the free energyis obtained for an equilibrium number density of n eq L = 4 π µ eff T . (6)Lepton number violation is mediated by the exchangeof heavy neutrinos. In contrast to thermal leptogen-esis [17], we will assume all heavy right-handed neu-trino masses to be close to the scale of grand unification(GUT), M i ∼ O (cid:0) − · · · (cid:1) Λ GUT ∼ · · · GeV,so that the heavy neutrinos are never produced ther-mally, i.e., T (cid:28) M i at all times. This assumption servesthe purpose to separate the mechanism under study fromthe contributions from ordinary thermal leptogenesis. Inthe expanding universe, the evolution of the L numberdensity n L is then described by the Boltzmann equation˙ n L + 3 Hn L (cid:39) − Γ L ( n L − n eq L ) , Γ L = 4 n eq (cid:96) σ eff , (7)where n eq (cid:96) = 2 /π T and with σ eff ≡ (cid:104) σ ∆ L =2 v (cid:105) denotingthe thermally averaged cross section of two-to-two scat-tering processes with heavy neutrinos in the intermediatestate that violate lepton number by two units,∆ L = 2 : (cid:96) i (cid:96) j ↔ HH , (cid:96) i H ↔ ¯ (cid:96) j ¯ H , (8) (cid:96) Ti = (cid:0) ν i e i (cid:1) , H T = (cid:0) h + h (cid:1) , i, j = 1 , , . We note that the term proportional to n eq L now acts as anovel production term for the lepton asymmetry, as longas the axion field is in motion. For center-of-mass ener-gies much smaller than the heavy neutrino mass scale, √ s (cid:28) M i , the effective cross section σ eff is practicallyfixed by the experimental data on the light neutrino sec-tor [18], assuming the seesaw mass matrix [19]: σ eff ≈ π ¯ m v (cid:39) × − GeV − , ¯ m = (cid:88) i =1 m i , (9)where v ew (cid:39)
174 GeV and where we have assumed thatthe sum of the light neutrino masses squared is of thesame order of magnitude as the atmospheric neutrinomass difference, ∆ m (cid:39) . × − eV [20].For a (cid:28) M Pl , and as long as H (cid:29) m a , i.e., priorto the onset of the axion oscillations, the axion energydensity ρ a is much smaller than the total energy den-sity ρ tot = ρ ϕ + ρ R + ρ a ≈ ρ ϕ + ρ R , where ρ ϕ and ρ R In the following, we shall approximate all number and en-ergy densities by their corresponding expressions in the classicalBoltzmann approximation. We will only take care of quantum-statistical effects when counting relativistic degrees of freedom. are the energy densities of the inflaton and of radiation.Reheating is described by a system of equations:˙ ρ ϕ + 3 Hρ ϕ = − Γ ϕ ρ ϕ , ˙ ρ R + 4 Hρ R = + Γ ϕ ρ ϕ , (10) H ≡ (cid:0) ˙ R/R (cid:1) = ρ tot M , ρ tot ≈ ( ρ ϕ + ρ R ) , (11)where Γ ϕ is the inflaton decay rate. The inflaton mustnot decay before the end of inflation, which impliesΓ ϕ (cid:46) H inf . (12)The rough temperature scale of leptogenesis as wellas the axion mass scale in our scenario are determinedby the requirement that the heavy neutrino-mediated∆ L = 2 interactions must be in thermal equilibriumbefore the onset of axion oscillations, Γ L (cid:29) H (cid:38) m a ,which yields T ∼ T L = g / ∗ / ( π σ eff M Pl ) ∼ GeVand m a ∼ σ eff T L ∼ GeV. Upon closer examination,the solution for the temperature, T ≡ π / /g ∗ ρ R , ac-cording to Eqs. (10) and (11) shows the following charac-teristic behavior: within roughly one Hubble time afterthe end of inflation, T quickly rises to its maximal value, T max (cid:39) × GeV (cid:18) Γ ϕ GeV (cid:19) / (cid:18) H inf GeV (cid:19) / , (13)after which the temperature decreases because the en-ergy density is dominated by the inflaton oscillations(which scale as matter). During reheating, the tempera-ture drops as T ∝ R − / until radiation comes to dom-inate at time t = t rh (cid:39) Γ − ϕ , when ρ R = ρ ϕ , and thereheating temperature is T rh (cid:39) × GeV (cid:18) Γ ϕ GeV (cid:19) / . (14)After the end of reheating, i.e., for t > t rh , the expansionis then driven by relativistic radiation and the tempera-ture simply decreases adiabatically, T ∝ R − . In the caseof a large axion decay constant, this phase of radiationdomination, however, does not last all the way to the timeof primordial nucleosynthesis. Instead, the axion comesto dominate the total energy density at some time priorto its decay, which marks the beginning of yet anotherstage of matter domination. The decay of the axion intorelativistic gauge bosons and the corresponding renewedtransition to radiation domination then represent a sec-ond installment of reheating, which can be described bythe same set of equations as the primary reheating pro-cess, cf. Eqs. (10) and (11). With the axion decay rateΓ a (cid:39) α π m a f a , α = g π , (15)and using Eq. (14), we find for the secondary reheatingtemperature or axion decay temperature T dec (cid:39) × GeV (cid:18) m a GeV (cid:19) / (cid:18) GeV f a (cid:19) . (16) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) m a (cid:60) (cid:71) (cid:106) m a (cid:62) (cid:71) (cid:106) (cid:68) a (cid:62) (cid:68) a (cid:61) (cid:71) (cid:106) (cid:62) H infmax m a (cid:62) H i n f m a x H infmax (cid:61) (cid:180) GeV f a (cid:61) (cid:180) GeV Axion mass m a (cid:64) GeV (cid:68) I n fl a t ond eca y r a t e (cid:71) (cid:106) (cid:64) G e V (cid:68) f a (cid:61) (cid:180) GeV, H infmax (cid:62) (cid:180) GeV R e h ea ti ng t e m p e r a t u r e T r h (cid:64) G e V (cid:68) FIG. 1: Contour plot of the final baryon asymmetry η B asa function of the axion mass m a and the inflaton decay rateΓ ϕ for an axion decay constant of f a = 3 × GeV. Theblack (bent) contours represent the full numerical result, whilethe colorful (straight) contours depict our analytical estimateaccording to Eqs. (18) and (19). In the lower part of theplot ( m a > Γ ϕ ), the effect of washout is illustrated by thedifference between the dashed ( κ = 0) and solid ( κ (cid:54) = 0) lines. This temperature should be at least of O (10) MeV [21],which imposes a lower bound on m a : m a (cid:38) × GeV (cid:18) f a GeV (cid:19) / . (17)The five differential equations in Eqs. (2), (7), (10),and (11) allow one to compute the present value of thebaryon asymmetry (i.e., the baryon-to-photon ratio), η B ≡ n B n γ = c sph g ∗ ,s g ∗ η aL (cid:39) . η aL , (18)where the sphaleron factor c sph accounts for the con-version of the lepton asymmetry into baryon asymme-try by sphalerons. Here, g ∗ ,s and g ∗ denote the effec-tive numbers of degrees of freedom contributing to theentropy density in the present epoch and during reheat-ing, respectively. In the standard model c sph = 28 / g ∗ ,s = 43 /
11, and g ∗ = 427 /
4. Last but not least, η aL inEq. (18) stands for the final lepton asymmetry after thedecay of the axion around t (cid:39) Γ − a .We determine η aL by solving the five differential equa-tions in Eqs. (2), (7), (10), and (11) numerically. Wealso present approximate analytical solutions, which willbe discussed in detail in an upcoming publication. It isconvenient to parametrize η aL as follows: η aL = C ∆ − a ∆ − ϕ η max L e − κ . (19)The approximate analytical results agree with the nu-merical results, as shown in Fig. 1. In the following, weshall present analytical expressions for the individual fac-tors on the right-hand side of Eq. (19). η max L denotes theall-time maximum value of the lepton asymmetry, whichis reached around the time when the axion oscillationsset it, i.e., at t ∼ t osc (cid:39) m − a . Note that this time doesnot necessarily coincide with the time when Γ L (cid:39) H .Integrating the Boltzmann equation for the lepton asym-metry up to t ∼ t osc , one approximately finds η max L (cid:39) σ eff g / ∗ a f a m a M Pl × min (cid:110) , (Γ ϕ /m a ) / (cid:111) , (20)which is suppressed w.r.t. the would-be equilibrium lep-ton asymmetry, η eq L = n eq L /n eq γ , evaluated at the sametime by a factor n L /n eq L (cid:39) T /T L × min (cid:8) , (Γ ϕ /m a ) / (cid:9) .Remarkably enough, the maximal lepton asymmetry israther insensitive to the axion decay constant; it only de-pends on the ratio a /f a , which is expected to be O (1).Furthermore, η max L turns out to be directly proportionalto the effective cross section σ eff . For a = f a and giventhe value of σ eff in Eq. (9), η max L is, hence, typically muchlarger than the observed value, η obs B (cid:39) × − [22], η max L (cid:39) × − (cid:18) m a GeV (cid:19) p (cid:18) Γ ϕ GeV (cid:19) q , (21)where the powers p and q are given as in Eq. (20).The two ∆ factors in Eq. (19) account for the entropyproduction in inflaton and axion decays during reheatingand at late times, respectively. We approximately have∆ ϕ (cid:39) max (cid:8) , ∆ (cid:48) ϕ (cid:9) , ∆ a (cid:39) max { , ∆ (cid:48) a } , (22)where ∆ (cid:48) ϕ (cid:39) ( m a / Γ ϕ ) / and with ∆ (cid:48) a being given as∆ (cid:48) a (cid:39) π α f a a m a M × min (cid:110) , (Γ ϕ /m a ) / (cid:111) . (23)In the region of parameter space in which we are able tosuccessfully reproduce η obs B , entropy production in axiondecays begins to play a role for f a values around 3 × GeV, cf. Fig. 2. For smaller values of f a , we alwayshave ∆ a = 1 in the entire parameter region of interest.The factor e − κ in Eq. (19) accounts for the washout of η max L during reheating due to the ∆ L = 2 washout pro-cesses, cf. the term proportional to − σ eff n L in Eq. (7).In case the axion begins to oscillate before the end ofreheating, i.e., for m a (cid:38) Γ ϕ , one can estimate κ ∼ T rh T L (cid:39) (cid:18) T rh GeV (cid:19) . (24)For m a (cid:46) Γ ϕ on the other hand, washout is always neg-ligible, so that κ can be safely set to κ = 0. A morecareful treatment of the effect of washout on the finalbaryon asymmetry in our scenario is left for future work. f a min (cid:62) (cid:180) (cid:32) GeV f a max (cid:62) (cid:180) (cid:32) GeV f a (cid:163) (cid:32) G e V f a (cid:61) (cid:32) G e V f a (cid:61) (cid:32) G e V Axion mass m a (cid:64) GeV (cid:68) I n f l a t ondecay r a t e (cid:71) (cid:106) (cid:64) G e V (cid:68) Η B (cid:61) Η B obs (cid:62) (cid:180) (cid:45) , f a min (cid:163) f a (cid:163) f a max R ehea ti ng t e m pe r a t u r e T r h (cid:64) G e V (cid:68) FIG. 2: Contour lines for successful leptogenesis ( η B = η obs B )in the m a –Γ ϕ plane for different values of the axion decay con-stant f a . The dashed segments along the individual contoursmark the regions where either m a or Γ ϕ become comparableto the maximally allowed Hubble rate H maxinf , cf. Eq. (1). For f a (cid:46) GeV, entropy production in axion decays ceases toaffect η B , cf. Eq. (23), which is reflected in the contour linesbeing no longer sensitive to changes in f a . The lower boundson m a and f a then directly follow from our constraints inEqs. (1), (3), and (12). At the same time, the regime of large f a , and hence the upper bounds on m a and f a , require a care-ful numerical analysis due to the strong impact of washout. Finally, the factor C in Eq. (19) is a numerical fudgefactor, which can, in principle, be estimated analytically,but which, in practice, is best determined by fitting η aL in Eq. (19) to the outcome of our numerical analysis.Specifically, we find C (cid:39) . m a (cid:46) Γ ϕ and C (cid:39) . m a (cid:38) Γ ϕ . The fact that these values are both of O (1)confirms the accuracy of our analytical estimate.Altogether, the parameter dependence of the final lep-ton asymmetry in Eq. (19) can be summarized as follows(here, we neglect the effect of washout and set κ → η aL ∝ T L m − / a Γ / ϕ a f − a ; m a (cid:38) Γ ϕ , ∆ a = 1 m a a f − a ; m a (cid:46) Γ ϕ , ∆ a = 1 m / a Γ / ϕ M a − f − a ; m a (cid:38) Γ ϕ , ∆ a > m a M a − f − a ; m a (cid:46) Γ ϕ , ∆ a > . (25)In the various regions of parameter space, typical valuesfor m a , Γ ϕ , and f a then yield the following final baryonasymmetries (again for a = f a and σ eff as in Eq. (9)), f a [GeV] m a [GeV] Γ ϕ [GeV] η B ∆ a × × × − × × × − × × × − × × × − f a . The range of allowed valuesspans five orders of magnitude,4 × GeV (cid:46) f a (cid:46) × GeV . (26)For smaller values of f a , it is not possible to generatea sufficiently large baron asymmetry, while keeping thebaryonic isocurvature perturbations small enough. Forlarger values of f a , the dilution of the asymmetry duringthe late-time decay of the axion is too strong. Varying f a within the interval in Eq. (26), we then find that m a ,Γ ϕ , and T rh can take values within the following ranges:1 × GeV (cid:46) m a (cid:46) × GeV , (27)3 × GeV (cid:46) Γ ϕ (cid:46) × GeV , (28)9 × GeV (cid:46) T rh (cid:46) × GeV . (29)These ranges of parameters are consistent with modelsof dynamical axions, as well as string axion models.Finally, let us conclude. As we have been able to show,lepton number violation due to the exchange of heavyMajorana neutrinos, in combination with the effectivechemical potential generated by a slowly rolling axion-like scalar field, is sufficient for a successful realizationof baryogenesis via leptogenesis. In this scenario, thebaryon asymmetry does neither depend on the concreteheavy neutrino mass spectrum nor on the amount of CP violation in the light and heavy neutrino sectors. In par-ticular, it is consistent with (almost) degenerate heavyneutrino masses close to the GUT scale. Hence, whilethermal leptogenesis assumes some of the heavy neutrinoYukawa couplings to be much smaller than O (1), ourmechanism equally applies in the case of Yukawa cou-plings of O (cid:0) − · · · (cid:1) . Furthermore, as can be seenfrom Eqs. (9) and (25), the baryon asymmetry increaseswith the light neutrino masses m i . Thus, while ther-mal leptogenesis imposes an upper bound on the neutrinomass scale, ¯ m (cid:46) . f a may potentially lie in the vicinity of the decayconstant of the QCD axion, (iii) the possibility of bary-onic isocurvature perturbations at a detectable level, etc. Such constraints may easily allow for a possibility to testour mechanism in the near future and assess whether ithas indeed the potential to serve as a viable explanationfor the baryon asymmetry of the Universe.The authors wish to thank K. Harigaya, A. Kamada,K. Kamada, M. Kawasaki, L. Pearce, and F. Takahashifor helpful comments and discussion. This work has beensupported in part by the U. S. Department of EnergyGrant de-sc0009937 (A. K.), by Grants-in-Aid for Scien-tific Research from the Ministry of Education, Culture,Sports, Science, and Technology (MEXT), Japan, No.26104009 and No. 26287039 (T. T. Y.), and by the WorldPremier International Research Center Initiative (WPI),MEXT, Japan. ∗ Corresponding author. E-mail: [email protected][1] R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. , 1440(1977).[2] S. Weinberg, Phys. Rev. Lett. , 223 (1978); F. Wilczek,Phys. Rev. Lett. , 279 (1978).[3] J. E. Kim, Phys. Rev. Lett. , 103 (1979); M. A. Shif-man, A. I. Vainshtein, and V. I. Zakharov, Nucl. Phys. B , 493 (1980); M. Dine, W. Fischler, and M. Srednicki,Phys. Lett. B , 199 (1981).[4] J. E. Kim, Phys. Rept. , 1 (1987).[5] M. B. Green and J. H. Schwarz, Phys. Lett. B , 117(1984).[6] P. Svrcek and E. Witten, JHEP , 051 (2006), hep-th/0605206.[7] L. E. Ibanez and H. P. Nilles, Phys. Lett. B , 354(1986).[8] T. S. Bunch and P. C. W. Davies, Proc. Roy. Soc. Lond.A , 117 (1978); A. D. Linde, Phys. Lett. B , 335(1982); A. A. Starobinsky and J. Yokoyama, Phys. Rev.D , 6357 (1994), astro-ph/9407016.[9] A. Dolgov and K. Freese, Phys. Rev. D , 2693 (1995),hep-ph/9410346; A. Dolgov, K. Freese, R. Rangarajan,and M. Srednicki, Phys. Rev. D , 6155 (1997), hep-ph/9610405.[10] T. Chiba, F. Takahashi, and M. Yamaguchi, Phys. Rev.Lett. , 011301 (2004), Erratum-ibid. , 209901(2015), hep-ph/0304102.[11] A. Kusenko, L. Pearce, and L. Yang, Phys. Rev.Lett. , 061302 (2015), arXiv:1410.0722 [hep-ph];L. Pearce, L. Yang, A. Kusenko, and M. Peloso,arXiv:1505.02461 [hep-ph]; L. Yang, L. Pearce, andA. Kusenko, arXiv:1505.07912 [hep-ph].[12] A. G. Cohen and D. B. Kaplan, Phys. Lett. B , 251(1987); Nucl. Phys. B , 913 (1988).[13] M. Dine, P. Huet, R. L. Singleton, and L. Susskind, Phys.Lett. B , 351 (1991); A. G. Cohen, D. B. Kaplan, andA. E. Nelson, Phys. Lett. B , 86 (1991).[14] V. A. Kuzmin, M. E. Shaposhnikov, and I. I. Tkachev,Phys. Rev. D , 466 (1992).[15] G. Servant, Phys. Rev. Lett. , no. 17, 171803 (2014),arXiv:1407.0030 [hep-ph].[16] P. J. E. Peebles, Nature , 210 (1987); Astrophys.J. Lett. , 73 (1987); K. Enqvist and J. McDonald, Phys. Rev. Lett. , 2510 (1999), hep-ph/9811412; Phys.Rev. D , 043502 (2000), hep-ph/9912478. K. Harigaya,A. Kamada, M. Kawasaki, K. Mukaida, and M. Yamada,Phys. Rev. D , 043510 (2014), 1404.3138 [hep-ph].[17] M. Fukugita and T. Yanagida, Phys. Lett. B , 45(1986).[18] W. Buchmuller, P. Di Bari, and M. Plumacher, AnnalsPhys. , 305 (2005), hep-ph/0401240.[19] T. Yanagida, in “Proceedings of the Workshop on theUnified Theory and the Baryon Number in the Universe” ,Tsukuba, Japan, Feb. 13-14 (1979), edited by O. Sawadaand A. Sugamoto, KEK report KEK-79-18, Conf. Proc. C , 95 (1979); Prog. Theor. Phys. , 1103 (1980);M. Gell-Mann, P. Ramond and R. Slansky, in ”Super-gravity” , North-Holland, Amsterdam (1979), edited byD. Z. Freedom and P. van Nieuwenhuizen, CERN re-port Print-80-0576, Conf. Proc. C , 315 (1979),arXiv:1306.4669 [hep-th]; see also: P. Minkowski, Phys.Lett. B , 421 (1977). [20] K. A. Olive et al. , [Particle Data Group Collaboration],Chin. Phys. C , 090001 (2014).[21] M. Kawasaki, K. Kohri and N. Sugiyama, Phys. Rev.Lett. , 4168 (1999), astro-ph/9811437; Phys. Rev.D , 023506 (2000), astro-ph/0002127; S. Hannestad,Phys. Rev. D , 043506 (2004), astro-ph/0403291.[22] P. A. R. Ade et al. , [Planck Collaboration], Astron. As-trophys. , A16 (2014), 1303.5076 [astro-ph.CO].[23] W. Buchmuller, R. D. Peccei, and T. Yanagida, Ann.Rev. Nucl. Part. Sci. , 311 (2005), hep-ph/0502169.[24] G. Drexlin, V. Hannen, S. Mertens, and C. Wein-heimer, Adv. High Energy Phys. , 293986 (2013),arXiv:1307.0101 [physics.ins-det].[25] K. N. Abazajian, E. Calabrese, A. Cooray, F. DeBernardis, S. Dodelson, A. Friedland, G. M. Fuller,and S. Hannestad, Astropart. Phys.35