Cosmic Inflation in Minimal $U(1)_{B-L}$ Model: Implications for (Non) Thermal Dark Matter and Leptogenesis

We study the possibility of realising cosmic inflation, dark matter (DM), baryon asymmetry of the universe (BAU) and light neutrino masses in minimal gauged $B-L$ extension of the standard model with three right handed neutrinos. The singlet scalar field responsible for spontaneous breaking of $B-L$ gauge symmetry also plays the role of inflaton by virtue of its non-minimal coupling to gravity. While the lightest right handed neutrino is the DM candidate, being stabilised by an additional $Z_2$ symmetry, we show by performing a detailed renormalisation group evolution (RGE) improved study of inflationary dynamics that thermal DM is generally overproduced due to insufficient annihilations through gauge and scalar portals. The non-thermal DM scenario is viable, with or without $Z_2$ symmetry, although in such a case the $B-L$ gauge sector remains decoupled from the inflationary dynamics due to tiny couplings. We also show that the reheat temperature predicted by the model prefers non-thermal leptogenesis while being consistent with light neutrino data as well as non-thermal DM scenario.


I. INTRODUCTION
Precision measurements of the cosmic microwave background (CMB) anisotropies by experiments like Planck [1][2][3] reveals that our universe is homogeneous and isotropic on large scale upto a remarkable accuracy. However, the observed isotropy of the CMB lead to the horizon problem which remains unexplained in the standard cosmology where the universe remains radiation dominated throughout the early stages. In order to solve the horizon problem, the presence of a rapid accelerated expansion phase in the early universe, called inflation [4][5][6] was proposed. Originally proposed to solve the horizon, flatness and unwanted relic problems in standard cosmology, the inflationary paradigm was also subsequently supported by the adiabatic and scale invariant perturbations observed in the CMB [1,2]. Such an early accelerated phase of expansion can be generated by the presence of one or more scalar fields whose dynamics crucially decides the period of inflation. Over the years, a variety of inflationary models have been studied with different levels of success [7]. The earliest proposal of this sort is known as chaotic inflation [8,9] where simple power law potentials like m 2 φ 2 with a scalar field φ were used. However, such simple models predict very specific values of inflationary parameters like the spectral index n s ∼ 0.967, tensor-to-scalar ratio r ∼ 0.133 for number of e-folds N e = 60 and unfortunately, the latest Planck 2018 data [2] strongly disfavour this simple model due to its large prediction of r. Modified chaotic inflation where the inflation sector is extended by an additional scalar field to assists the inflaton field has also been proposed [10][11][12]. Another class of models use the Higgs as the inflaton [13,14].
These models often suffer from problems of vacuum stability [15] and nonunitarity [16] as well as being inadequate for combining inflation with other cosmological problems like DM and baryogenesis. A possible way out is to consider an extra stabilising scalar which acts as the inflaton. We consider this possibility in our work where an additional scalar with nonminimal coupling to gravity [17][18][19][20], in addition to usual quartic chaotic type couplings, can give rise to successful inflation while predicting the inflationary parameters within observed range. The same scalar field is also responsible for several other interesting phenomenology as we discuss below.
The same CMB measurements mentioned above also suggest that the present universe has a significant amount of non-luminous, non-baryonic form of matter, known as dark matter (DM) [3,21]. This is also supported by astrophysical evidences gathered over a much longer period of time [22][23][24]. The Planck 2018 data reveal that approximately 26% of the present universe is composed of DM, which is about five times more than the ordinary luminous or baryonic matter. In terms of density parameter Ω DM and h = Hubble Parameter/(100 km s −1 Mpc −1 ), the present DM abundance is conventionally reported as [3]: Ω DM h 2 = 0.120 ± 0.001 at 68% CL. Since none of the standard model (SM) particles can satisfy the criteria of a DM candidate several proposals have been put forward among which the weakly interacting massive particle (WIMP) is perhaps the most popular one. In this framework, a DM particle having mass and interactions typically around the electroweak scale can give rise to observed abundance after thermal freeze-out, a remarkable coincidence often referred to as the WIMP Miracle [25]. The same interactions responsible for thermal freeze-out of WIMP type DM should also give rise to sizeable DM-nucleon scattering. However, null results at direct detection experiments like LUX [26], PandaX-II [27,28], XENON1T [29,30] have certainly pushed several WIMP models into a tight corner, if not ruled out yet. This has also generated interests in beyond thermal WIMP paradigm as viable alternatives. One such interesting possibility is the non-thermal origin of DM [31].
For a recent review of such feebly interacting (or freeze-in) massive particle (FIMP) DM, please see [32]. In the FIMP scenario, DM candidate does not thermalise with the SM particles in the early universe due to its feeble interaction strength and the initial abundance of DM is assumed to be zero. At some later stage, DM can be produced non thermally from decay or annihilation of other particles thermally present in the universe.
Similarly, the baryonic content of the universe also gives rise to another puzzle due to the abundance of baryons over antibaryons. Quantitatively, this excess is denoted as baryon to entropy ratio [3,21] where Y B denotes comoving baryon density, n B (nB) denotes baryon (anti-baryon) number density while s is the entropy density. Since any initial asymmetry before inflation would be washed out at the end of inflation due to the exponential expansion of the universe, there has to be a dynamical mechanism to generate the asymmetry in a post-inflationary universe.
This requires certain conditions, known as the Sakharov conditions [33] to be fulfilled. They are namely, baryon number (B) violation, C and CP violation and departure from thermal equilibrium, not all of which can be fulfilled in the required amounts within the SM alone.
Generation of baryon asymmetry of the universe (BAU) from out-of-equilibrium decays of heavy particles has been a well-known mechanism for baryogenesis [34,35]. Another interesting way, which also connects the lepton sector physics, is known as leptogenesis which was proposed a few decades back [36]. In leptogenesis, instead of creating a baryon asymmetry directly from B violating interactions, an asymmetry in lepton sector is created via lepton number (L) violating processes (decay or scattering). If this lepton asymmetry is generated before the electroweak phase transition (EWPT), then the (B + L)-violating EW sphaleron transitions [37] can convert it to the required baryon asymmetry. Since the quark sector CP violation is insufficient to produce the required baryon asymmetry, the mechanism of leptogenesis can rely upon lepton sector CP violation which may be quite large as hinted by some neutrino oscillation experiments [38,39]. An interesting feature of this scenario is that the required lepton asymmetry can be generated through CP violating out-of-equilibrium decays of the same heavy fields that take part in popular seesaw mechanisms [40][41][42][43][44][45] that also explains the origin of tiny neutrino masses [21], another observed phenomenon which the SM fails to address.
Motivated by these, we study a minimal extension of the SM, by a gauged B−L symmetry with three right handed neutrinos (RHN) required to cancel the anomalies and a singlet scalar to break the additional gauge symmetry spontaneously, apart from generating RHN masses. While the singlet scalar also plays the role of inflaton, one RHN is stabilised by an additional Z 2 symmetry to become a DM candidate. The other two RHNs can give rise to light neutrino masses with vanishing lightest neutrino mass apart from producing the required lepton asymmetry which gets converted into the observed baryon asymmetry via sphalerons. Interestingly, we find that the strict limit on inflationary observables from Planck 2018 and BICEP 2 / Keck Array (BK15) data [2] as well as the stability of inflaton potential restrict the B − L gauge couplings, scalar couplings and Yukawa couplings associated with the inflaton field to be within some limits which do not favour thermal DM scenario due to insufficient annihilations. On the other hand, with very tiny gauge and Yukawa couplings one can realise the non-thermal DM scenario (with or without Z 2 symmetry) while the inflationary potential behaviour merges with the usual case of quartic Higgs inflation with non minimal coupling to gravity. We also find that the predicted values of reheat temperature makes it difficult to realise high scale thermal N 2 leptogenesis [46,47] leaving the option of non-thermal leptogenesis [48][49][50][51][52][53][54][55][56] viable.
The structure of the paper is organised as follows. In section II, we discuss the particle contents of the proposed set up and their interactions followed by briefly mentioning the existing constraints in section III. In section IV we perform a detailed study of inflation and its predictions in view of Planck 2018 bounds. We discuss different aspects of DM phenomenology in section V and then move onto discussing the possibility of non-thermal leptogenesis in section VI. Finally we conclude in section VII.

II. THE MODEL
As mentioned earlier, we study a gauged B − L extension of the SM with the minimal field content which can give rise to cancellation of triangle anomalies, spontaneous gauge symmetry breaking, light neutrino masses, dark matter, leptogenesis and cosmic inflation.
While gauged B − L symmetric extension of the SM was proposed long ago [57][58][59][60][61][62], realising a stable DM candidate in the model requires non-minimal field content or additional discrete symmetries. Also, a gauged B − L model with just SM fermion content, is not anomaly free due to the non-vanishing triangle anomalies for both both U (1) 3 B−L and the mixed U (1) B−L − (gravity) 2 . These triangle anomalies for the SM fermion content are given as Remarkably, if three right handed neutrinos are added to the model, they contribute The gauge invariant Lagrangian of the model is where L SM denotes the SM Lagrangian involving quarks, gluons, charged leptons, left handed neutrinos and electroweak gauge bosons while the second term is the kinetic term of B − L 1 For other exotic solutions to anomaly cancellation conditions, see [63][64][65][66][67][68][69].  gauge boson (Z BL ) expressed in terms of field strength tensor B αβ = ∂ α Z β BL − ∂ β Z α BL . The gauge invariant scalar Lagrangian of the model is as follows where The covariant derivatives of scalar fields are with g 1 and g 2 being the gauge couplings of SU (2) L and U (1) Y respectively and W a µ (a = 1, 3) and B µ are the corresponding gauge fields. On the other hand Z BL , g BL are the gauge boson and gauge coupling respectively for U (1) B−L gauge group.
The gauge invariant fermionic Lagrangian of the model is as follows The covariant derivative is defined as with Q R = −1 being the B − L charge of right handed neutrinos. Due to the presence of Z 2 symmetry, N R 1 has no mixing with N R 2,3 and also does not interact with SM leptons thereby qualifying for a stable DM candidate.
After breaking of both B − L symmetry and electroweak symmetry by the vacuum expectation value (VEV)s of H and Φ, the doublet and singlet scalar fields are given by where v and v BL are VEVs of H and Φ respectively. The right handed neutrinos and Z BL get masses after the U (1) B−L breaking as Here we consider diagonal Yukawa Y N in (N R 2 , N R 3 ) basis. Using equation (11) and equation (12), it is possible to relate M Z BL and M N i by Also after the breaking of SU (2) L × U (1) Y × U (1) B−L , the scalar fields h and φ can be related to the physical mass eigenstates H 1 and H 2 by a rotation matrix as where the scalar mixing angle θ is represented by .
The mass squared eigenvalues of physical scalars are given by Here M H 1 is identified as the SM Higgs mass whereas M H 2 is the singlet scalar mass.
One of the strong motivations of the minimal U (1) B−L model is the presence of heavy RHNs which can yield correct light neutrino mass via type I seesaw mechanism. The analytical expression for the light neutrino mass matrix is We consider the right handed neutrino mass matrix M N to be diagonal. Since in our case N R 1 does not interact with SM leptons, the lightest active neutrino would be massless. The Dirac neutrino Yukawa matrix Y D can be formulated through the Casas-Ibarra parametrisation [70] as where m d ν , M R are the diagonal light and heavy neutrino mass matrices respectively and U PMNS is the usual Pontecorvo-Maki-Nakagawa-Sakata (PMNS) leptonic mixing matrix. In the diagonal charged lepton basis, the PMNS mixing matrix is also the diagonalising matrix of light neutrino mass matrix In the above Casas-Ibarra parametrisation R represents a complex orthogonal matrix (RR T = I). In case of only two right handed neutrinos, the R matrix is a function of only one complex rotation parameter z = z R + iz I , z R ∈ [0, 2π], z I ∈ R [71]. For three right handed neutrinos taking part in seesaw mechanism R can depend upon three complex rotation parameters. Assuming one of them (rotation in 1-2 sector) to be vanishing, it can be represented as 2 Therefore with suitable choices of γ and γ , the Yukawa matrix can take different forms.
Furthermore, we shall use the best fit values of all three mixing angles and the mass squared differences of active neutrinos assuming a normal ordering [21].

III. CONSTRAINT ON THE MODEL PARAMETERS
In this section, we briefly discuss the theoretical and experimental constraints on different parameters in the model.
To begin with, we consider the bounded from below criteria of the scalar potential. This gives rise to following conditions to be satisfied by the quartic couplings On the other hand, to avoid perturbative breakdown of the model, all dimensionless couplings must obey the following limits at any energy scale: The non-observation of the extra neutral gauge boson in the LEP experiment [72,73] invokes following constraint on the ratio of M Z BL and g BL : The corresponding bounds from the large hadron collider (LHC) experiment have become stronger than this by now as both ATLAS and CMS collaborations have performed dedicated searches for dilepton resonances in proton proton collisions. The latest bounds from the ATLAS experiment [74,75] and the CMS experiment [76] at the LHC rule out such gauge boson masses below 4-5 TeV from analysis of 13 TeV centre of mass energy data. However, such limits are derived by considering the corresponding gauge coupling g BL to be similar to the ones in electroweak theory and hence the bounds become less stringent for weaker gauge couplings [74]. Additionally, if such Abelian gauge bosons couple only to the third generation leptons, then the collider bounds get even weaker, as explored recently in a singlet-doublet fermion DM scenario by the authors of [77].
Additionally, the singlet scalar of the model is also constrained [78,79]

IV. INFLATION
In this section, we describe the dynamics of inflation in detail and its predictions in view of the present experimental bounds. We identify the real part of singlet scalar field Φ as the inflaton. Along with the renormalisable potential in equation (5), we also assume that Φ is non-minimally coupled to gravity. For earlier studies in this context, please see [84,85] and references therein. For works guided by the same unifying principle of inflation, dark matter and neutrino mass, one may look at [86][87][88][89] as well as references therein.
We denote the inflaton field as φ hereafter, which is same as the notation used for real part of Φ filed in earlier sections. Thus the potential responsible for inflation is given by where R stands for the Ricci scalar and ξ is a dimensionless coupling of singlet scalar to gravity. We have neglected the contribution of v BL in equation (22) by considering it to be much lower than the reduced Planck mass M P . The action for φ in Jordan frame takes the following form (apart from the couplings to the fermions and SM Higgs) where Ω(φ) 2 = 1 + ξφ 2 M 2 P , g is the spacetime metric in the (−, +, +, +) convention, D µ φ stands for the covariant derivative of φ containing couplings with the gauge bosons which just reduces to the normal derivative D µ → ∂ µ (since during inflation, there are no fields other than the inflaton).
In order to simplify the calculations we make the following conformal transformation to write the action S J in the Einstein frame [90,91]: so that it looks like a regular field theory action with no explicit couplings to gravity. In the above transformation,ĝ represents metric in the Einstein frame. To make the kinetic term of the inflaton canonical, we redefine φ by where χ is the canonical field. Using these inputs, the inflationary potential in the Einstein frame can be written as where (22). We then make another redefinition: and reach at a much simpler from of V E given by Note that for an accurate analysis, one should work with renormalisation group (RG) improved potential and in that case, λ 2 in equation (27) will be function of Φ such that, The one loop renormalisation group evolution (RGE) equations of the relevant parameters associated with the inflationary dynamics are given by, where we define s = 1 + ξφ 2 The RGE equations for rest of the couplings are provided in Appendix A.
We choose the heavy neutrino mass spectrum to be satisfying the hierarchy M N 1 M N 2 < M N 3 and a diagonal mass matrix in the (N 2 , N 3 ) basis. Note that, from this section onwards, we are denoting the RHNs as N i only without denoting the chirality explicitly.
Let us first analyse the case where the RG running of λ 2 is dominated by g BL and Y N 2,3 . Then equation (29) can be rewritten as We ignore the contributions of λ 2 and Y N 1 in the R.H.S. of equation (33) considering them to be negligible 3 . Since λ 2 is very small, β λ 2 0 or β λ 2 0 can cause sharp changes in λ 2 value from its initial magnitude during the evolution. It may also happen that λ 2 becomes negative at some energy scale. Then the inflationary potential would turn unstable along φ field direction. Therefore the most acceptable case is to make β λ 2 → 0 at least during inflation so that the inflationary potential remains stable [85]. To ensure β λ 2 0, We can further simplify the expression for ∆ by assuming λ 2 by an amount of ±10%. Fig. 1 clearly points out that indeed a small violation of the ∆ ∼ 0 criteria can cause sharp instability of the inflationary potential. In upper left panel of Fig. 2, we exhibit the behaviour of the inflationary potential V E against Φ for different values of g BL considering ξ = 0.1. The value of Σ 4 N is determined from the equality ∆ earlier defined. As it can be observed, with the increase of g BL , the potential starts to develop a local minimum near some Φ value say, Φ I . If such a local minimum exists, then the field could be trapped there and the inflaton will stop rolling. This provides an upper bound on g BL such that the local minimum of V E (Φ) does not appear. The existence of a local minimum can be further confirmed if dV E (Φ) dΦ 0 near Φ I . This condition can be converted as Fig. 2 as a function of Φ. We observe that for g BL g max BL = 0.045, the inflationary potential indeed develops a local minimum near Φ I = 4M P . Similar conclusion can be drawn for ξ = 1 as shown in lower panel of Fig. 2 .
One important point to be noted is that the value of g max BL gets enhanced with the increase of ξ. We illustrate this in Fig. 3 where g max BL is plotted against different values of ξ. Next, we move on to calculate the predictions for inflationary observables. In terms of the original field φ, the slow roll parameters ( , η) and number of e-folds (N e ) are found to respectively. The inflationary observables such as spectral index (n s ), tensor to scalar ratio (r) and scalar perturbation spectrum (P S ) can be expressed in terms of the slow roll parameters: All these quantities have to be determined at the horizon exit of the inflaton (φ t ) and we consider the number of e-foldings N e = 60 for the numerical analysis. We perform a numerical scan over g BL and ξ to estimate the inflationary observables n s and r considering ∆ ∼ 0. The initial value of λ 2 is determined to produce the correct observed value of scalar perturbation spectrum P S at horizon exit. In Fig. 4 we show the variation of λ 2 with ξ to be consistent with the observed value of P S = 2.4 × 10 −9 . It turns out that value of r does not change much with the variation of g BL for a constant value of ξ since β λ = 0 at inflationary energy scale. Contrary to this, value of n s is quite sensitive to g BL . We see from left panel of Fig. 5 that n s increases with the enhancement of g BL for different values of ξ.
The rate of increase of n s with g BL turns flatter with the rise of ξ value. In right panel of For comparison purpose we also insert the Planck 2018+BAO+BK15 1σ and 2σ bounds [2].
It is evident that the present set up is able to provide set of n s − r values, consistent with the experimental constraints. Finally, in left panel Fig. 6, we constrain the ξ − g BL plane which correctly produces the n s − r values consistent with Planck 1σ (red) and 2σ (brown) bounds. inflaton to gravity as originally studied in [92]. For completeness purpose we discuss this particular case in right panel of Fig. 6 in n s − r plane. As it is seen the n s − r contour can still satisfy the Planck 2018 1σ bounds for N e = 60. The contour of observed value of P S in ξ − λ 2 plane remains same as in Fig. 4.

A. Reheating
Once inflation ends, the thermalisation of the universe is of utmost importance, leading to a radiation dominated universe. This is the reheating epoch [93], which takes the universe from the matter-dominated phase during inflation to the radiation-domination phase. We assume that reheating of the universe occurs through the instantaneous decay of inflaton For the first case g BL is large and thus ∆ ∼ 0 is an essential condition for the stability of inflationary potential. We consider λ 3 g 2 BL so that it does not effect the evolution of ∆ significantly. This assumption was made earlier also while determining the fate of inflation. The value of ∆ as defined earlier changes at a small amount in its RG evolution (see right panel of Fig. 7). It is found that the value of λ 2 (Φ I ) switches its order at low scale, for example λ 2 (Φ = 1 TeV) becomes O(10 −6 ) from 4.34 × 10 −10 at inflationary scale where Γ H 2 is the total decay width of inflaton. Suppose we consider M Z BL = 1 TeV, M N 1 = 500 GeV along with λ 2 = 4.35 × 10 −10 for ξ = 1 at inflationary energy scale. The mass of other two heavy right handed neutrinos will be fixed by the equality related to ∆.
In Fig. 8, we show contours of different reheating temperatures in λ 3 − g BL plane. It is seen that T R is enhanced with the increase of λ 3 . This is because the mixing between inflaton and SM Higgs is directly proportional to λ 3 as given in equation (15), and hence increase in In the second case g 4 BL , Σ 4 N λ 2 2 , the inflationary potential is mainly driven by λ 2 with other couplings sufficiently small. Hence ∆ ∼ 0 is not a necessary condition for this case.
However the coupling λ 3 (we take O(10 −10 )) should be still much smaller than unity so that the stability of inflation potential remains intact. Here, due to the smallness of all relevant couplings there will not be any significant changes during their RG running unlike to the earlier case. Now, depending on the mass scale (or λ 2 ), the tree level decay of inflaton to Z BL Z BL , H 1 H 1 final states are possible. The inflaton can also decay to right handed neutrinos, if kinematically allowed. In Fig. 9  coloured region is ruled out from the requirement of predicting correct observed value of scalar perturbation spectrum P S at horizon exit discussed earlier. In the blue coloured region inflaton mass turns larger than the reheating temperature and hence remains out of equilibrium. This may have important implications for other related phenomenology as we will discuss in a while.

V. DARK MATTER
In this section, we discuss the dark matter phenomenology in detail and attempt to find its consistency with the inflationary dynamics. As earlier pointed, N 1 is the DM candidate which is odd under Z 2 and hence stable. For earlier studies of DM in this model, one may refer to [95][96][97][98]. While the lightest right handed neutrino is the DM candidate, the other two RHN's take part in the usual type I seesaw mechanism, giving rise to solar and atmospheric neutrino mixing. Since DM is a singlet under SM gauge symmetry, it can interact with the visible sector particles only via gauge (Z BL ) or scalar (H 2 ) interactions. Now, depending upon the two cases namely, (i) g 4 BL , Σ 4 N λ 2 2 and (ii) g 4 BL , Σ 4 N λ 2 2 discussed in the context of inflation, DM-SM couplings can either be of order unity or very small. This will lead to completely different DM phenomenology namely, thermal or WIMP type and non-thermal or FIMP type, which we discuss separately below.

A. WIMP DM Scenario
For the first case that is, g 4 BL , Σ 4 N λ 2 2 discussed in the context of inflation earlier, it is expected that the DM stays in thermal equilibrium with the SM particles in the early dynamics, the one loop mixing can be neglected in comparison to other relevant couplings and processes. Therefore, for simplicity, we ignore such kinetic mixing for the rest of our analysis.
The evolution of comoving number density of DM (Y DM = n DM /s) is determined by the corresponding Boltzmann equation where with g and g * s being the internal degrees of freedom of the dark matter and relativistic entropy degrees of freedom respectively and z = M N 1 /T . The σv in equation (40) stands for the thermally averaged cross section of DM annihilation, given by [100] σv = 1 We have considered λ 2 = 5 × 10 −10 and λ 3 ∼ 10 −6 at inflationary energy scale.
We implement the model in FeynRules [101] and then use micrOMEGAs package [102] to estimate the relic abundance of DM numerically. The independent parameters which participate in determining the DM relic abundance are the following: In our case, we have considered the Y N matrix diagonal and the sum of fourth power of each running. Since Y N 1 is taken to be smaller than Y N 2,3 , DM mass M N 1 is smaller than M N 2,3 's.
We have already discussed the Dirac neutrino Yukawa or Y D matrix and here we use the same form as defined in equation (19) using Casas-Ibarra parametrisation. Here we work with λ 2 = 4.34 × 10 −10 (corresponding to ξ = 1, see Fig. 4), λ 3 = 10 −6 at inflationary energy scale. Since in our working range of gauge coupling 0.01 < g BL < 0.075, the T R is found to be around 10 6−7 GeV considering λ 3 = 10 −6 , it is expected that for M Z BL ∼ (1 − 10) TeV (see Fig. 8), all relevant SM and BSM fields will maintain thermal equilibrium with each other.
In Fig. 11 Fig. 11 represents the observed relic abundance, as per Planck 2018 data [3]. It is seen that the annihilation through gauge boson is the most efficient one and can satisfy correct relic in two out of three scenarios discussed. We then perform a numerical scan to find the parameter space satisfying correct DM relic.
In Fig. 12, we display the points satisfying correct DM relic (black dots) in M Z BL −g BL plane considering M Z BL 10 TeV. We use the values of relevant parameters as earlier mentioned.
We also include the LHC bound from dilepton resonance searches [74] (red curve), Planck constraints on inflation and stability bounds of the inflationary potential for comparison purpose. The shaded regions are disfavoured from the respective constraints. To conclude, we observe that with TeV scale or lower Z BL mass, it is unlikely to obtain the correct value of relic abundance for WIMP dark matter while being in agreement with LHC and inflationary observables simultaneously. We also check that direct detection limits on spin-independent DM-nucleon cross section from the XENON1T experiment [29,30] do not put any additional constraint on this parameter space as all the points shown in Fig. 13 obey these bounds.

B. FIMP DM Scenario
In the second case (g 4 BL , Σ 4 N λ 2 2 ), the couplings responsible for DM-SM interactions are tiny and hence it is expected that DM may never reach thermal equilibrium with the standard bath. This falls under the ballpark of FIMP dark matter, discussed earlier. For earlier work on fermion singlet as FIMP DM in U (1) B−L model, please see [103,104] and references therein. A recent study also discussed the possibility of scalar singlet responsible for breaking B − L gauge symmetry spontaneously to be a long-lived FIMP DM candidate [105]. If N 1 is a FIMP candidate, it can be produced non-thermally, due to decay or annihilation of other particles. In case Z 2 symmetry is exact, N 1 will be only pair produced as it is the only Z 2 odd particle. All scattering processes shown in Fig. 10 while discussing WIMP scenario can potentially contribute to the production of FIMP DM as well, when considered in the reverse direction. In addition, decays of H 1,2 and Z BL , if kinematically allowed, can also contribute to the relic density of N 1 . Typically, if same dimensionless couplings govern the strength of both decay and annihilation processes, the former dominates simply due to power counting. This is precisely the scenario here and FIMP is primarily produced from decays.
For our numerical calculation, we choose λ 2 ∼ 1.04 × 10 −12 at inflationary energy scale corresponding to ξ ∼ 0.01 from inflationary requirements (see Fig. 4). Then from Fig. 9, it is evident that for this choice of λ 2 , H 2 would be in thermal equilibrium with other SM particles by virtue of its coupling with Higgs as well as heavy right handed neutrinos N 2,3 which also maintain equilibrium since their masses considered here are below T R and they can interact to SM fields through Yukawa interaction. We would like to keep λ 3 ∼ 10 −10 extremely small so that it does not alter the RG running of λ 2 during inflation. Since g BL is also very small to justify FIMP nature of DM, we will investigate the possibility of production of N 1 DM from non thermal tree level decay of Z BL and H 2 (see Fig. 14). We will consider two benchmark choices of M Z BL < 10 TeV for the analysis. It is to be noted that Z BL which interacts only via gauge coupling g BL is also expected to be out of equilibrium.
Hence non-thermal production of Z BL from other bath particles and its subsequent decay into N 1 pairs play non-trivial roles. We therefore use coupled Boltzmann equations for both Z BL and N 1 to calculate the relic abundance of N 1 in this scenario.
FIG. 14. DM production channels from tree level decay of heavier particles.
The evolution of the comoving number densities for Z BL and DM are governed by the following coupled Boltzmann equations [103] where z = M H 1 /T and x represents all possible initial states. g * (z) is defined by while g * s is same as defined earlier. Here, g ρ (x) denotes the effective number of degrees of freedom related to the energy density of the universe at z. The Γ A→BC denotes the thermally averaged decay width which is given by Since initial densities of both Z BL and N 1 are almost vanishing, one can ignore Y Z BL and Y DM from first term within each bracket on right hand side of equations (44) and (45).
In left panel of Fig. 15, we show the evolution of Y Z BL against z for benchmark choices of g BL and other relevant parameters indicated in the figure. It is seen that Y Z BL starts from a vanishingly small value initially and reaches a sizeable value with the lowering of temperature very quickly. The initial increase in Z BL abundance happens primarily from H 2 decays. As expected, the production of Z BL from H 2 decay becomes efficient around suppression in the equilibrium abundance of H 2 which makes Z BL production less efficient leading to the plateau region where Y Z BL remains more or less constant. While production from H 1 decay will be mixing suppressed due to smallness of λ 3 , it is also kinematically disallowed for the chosen benchmark value of Z BL mass. For some epochs the abundance of Z BL remains constant (denoted by the plateau region) and then goes down towards zero again due to subsequent decays of Z BL into N 1 as well as other lighter particles. Similar features can be observed in right panel of Fig. 15 where the evolution of N 1 abundance is shown. The N 1 abundance begins from vanishingly small value and gets enhanced due to non-thermal production from Z BL and H 2 decays and finally gets saturated.
We use two different values of g BL as shown in the figure. For higher g BL value, the amount of production and the final abundance of the DM from Z BL decay is larger, as expected. Since the DM mass is small, the corresponding Yukawa coupling with H 2 is tiny and hence direct production of DM is primarily dominated from Z BL decay. DM production can however, increase with increase in H 2 mass. This is because H 2 decay width to Z BL pairs increases  with increase in H 2 mass which eventually increases the non-thermal production of Z BL and hence DM. Once the freeze-in abundance of DM that is Y DM saturates, one can obtain the present relic abundance using the following expression: Here Ω DM = ρ DM ρc , where ρ DM is the DM energy density and ρ c = 8πG N is the critical energy density of the universe, with G N being Newton's gravitational constant and H 0 ≡ 100 h km s −1 Mpc −1 is the present-day Hubble expansion rate. Using the above expression in Fig. 16, we have shown the relic density evolutions of the DM for two set of parameters (with lower and comparatively larger M Z BL ). Note that the final DM abundance for both sets of parameters satisfy the present constraint given by Planck 2018 data. In table III we have tabulated the two sets of parameters used in Fig. 16 as well as the left panel of Fig. 15. As mentioned earlier, for such benchmark values of parameters the contribution of 2 → 2 scattering processes to DM production in the present analysis remains sub-dominant or negligible. It should be noted that while the required FIMP DM relic abundance can be successfully generated in this model, the corresponding parameter space leads to decoupling of B −L gauge sector from inflationary dynamics leading to a usual quartic plus non-minimal inflation [92].
So far, the analysis on non thermal production of dark matter is performed by assuming H 2 in thermal equilibrium with the SM bath. This is possible when M H 2 < T R and M H 2 has sizeable couplings with other particles in the bath. However, it is also possible that M H 2 remains larger compared to the reheat temperature M H 2 > T R and hence the inflaton remains out of equilibrium afterwards (see blue coloured region of Fig. 9). In such a case, the production of Z BL and N 1 will not be possible like the way it was discussed before. Since SM Higgs mixing with H 2 is also very small, it is not possible to generate correct FIMP abundance. While interactions by virtue of gauge coupling and Yukawa coupling with H 2 are insufficient to produce correct FIMP abundance, one can turn to Yukawa couplings with ordinary leptons which are present in thermal bath for most of the epochs. However one has to get rid of the Z 2 symmetry in order to introduce such Yukawa couplings through SM Higgs. We briefly discuss this possibility in the remainder of this section.
Once the Z 2 symmetry is discarded, one can have new non-diagonal terms in right handed neutrino mass matrix. However, for simplicity we continue to choose a diagonal right handed neutrino mass matrix or the corresponding Yukawa coupling matrix Y N . The newly introduced Yukawa couplings of N 1 to SM leptons can be written as This will generate mixing of N 1 with active neutrinos once the electroweak symmetry is broken. Using Casas-Ibarra parametrisation of equation (19) and using the form of complex orthogonal matrix given in equation (20), the Yukawa coupling of N 1 with leptons can be expressed as where γ is a complex angle and m 3 the heaviest active neutrino mass with normal ordering.
In deriving this, we fix Dirac CP phase to be zero 4 and also considered the lightest active neutrino as massless. The requirement of lightest active neutrino mass to be vanishingly small arises due to tiny Yukawa couplings of N 1 to leptons for being a FIMP DM. We define the mixing of sterile N 1 with i th active neutrino by: For simplicity, we redefine δ 1 = δ and the relation between δ and δ 2,3 can be easily found using equation (50). Owing to this tiny but non-zero mixing, N 1 can now interact with SM bath directly without relying upon Z BL or H 2 mediation considered earlier in Z 2 symmetric scenario. For example, W ± boson can directly decay to N 1 through W ± → N 1 α ± , α ≡ (e, µ, τ ) if kinematically allowed. The contribution from annihilation processes continues to be sub-dominant like before. The evolution of DM comoving number density is governed by Where we have considered only the most dominant decay modes. Decay channels with more than one N 1 in final state will be suppressed due to higher powers of tiny mixing δ.
Once we obtain Y DM , it is simple to compute the relic density of the DM using equation It is to be noted that, unlike the WIMP scenario, we are not performing a complete scan of parameter space for FIMP which can be found elsewhere. We have considered two possibilities based on inflaton mass being smaller or larger compared to reheat temperature and showed that required FIMP DM abundance can be successfully produced in both the scenarios. In the case where inflaton mass is larger compared to reheat temperature so that it is not present in the thermal bath afterwards, we find that the correct FIMP abundance can be produced only when we discard the Z 2 stabilising symmetry of DM and allow for more possibilities of its production from SM bath to open up. On the other hand, such longlived dark matter can have very interesting consequences at indirect detection experiments, which have been summarised in the review article [106].

VI. LEPTOGENESIS
In this section, we briefly discuss the possibilities of generating the observed baryon asymmetry of the universe through leptogenesis. Since the lightest right handed neutrino is our DM candidate, the required lepton asymmetry can be generated only by the out of equilibrium decays of heavier right handed neutrinos N 2,3 . Usually, in such type I seesaw framework, the requirement of producing the correct lepton asymmetry pushes the scale of right handed neutrinos to a very high scale M > 10 9 GeV, known as the Davidson Also, thermal leptogenesis is not affected much by inflationary dynamics at high scale. It is of course possible to realise thermal leptogenesis and non-thermal DM in this model, but we focus mainly on non-thermal leptogenesis due to its connection to inflation as well as reheat temperature as discussed below. In fact, thermal vanilla leptogenesis is not possible in our setup as the predicted values of reheat temperature discussed earlier (see Fig. 8 and Fig.   9) falls below the Davidson-Ibarra limit on scale of such leptogenesis. This motivates us to discuss non-thermal leptogenesis in this section.
The scenario of non-thermal leptogenesis [48][49][50][51][52][53][54][55][56] arises when the reheat temperature after inflation is lower than the masses of right handed neutrinos. Thus, although the right handed neutrinos can be produced due to the decay of inflaton, they can not reach thermal equilibrium with the SM particles due to insufficient reheat temperature. The nonequilibrium abundance of right handed neutrinos will be purely decided by their couplings to inflaton which will affect the final CP asymmetry generated by subsequent decays of right handed neutrinos. Since inflaton also has to decay into other SM bath particles reproducing a radiation dominated universe, one has to solve coupled Boltzmann equations involving inflaton, right handed neutrinos and SM radiation. However, for simplicity, we assume that the decay width of N 2,3 's (Γ N 2,3 ) to be larger than that of the inflaton (Γ H 2 ) so that decays of N 2,3 to SM particles can be instantaneous [54]. This allows us to retain the same reheating description (from inflaton decay only) discussed earlier. Thus, the right handed neutrinos produced from inflaton decay turns non-relativistic and decays to SM leptons and Higgs instantaneously. The CP asymmetry generated by N i decays, following the notations of [53], can be formulated as where the first and second terms in equation (54) are the individual contributions of N 2 and N 3 respectively. The loop function G(x) containing both self-energy and vertex corrections is defined as Once the CP asymmetry parameter is calculated, the comoving lepton asymmetry (ratio of excess of leptons over antileptons and entropy) can be calculated as where Br i represents the branching ratio of the inflaton decay to N i . Finally, the baryon asymmetry generated through the standard sphaleron conversion processes is given by We have used the Casas-Ibarra parametrisation of Y D as given by equation (19). Since lepton assymmetry gets generated from N 2 and N 3 decays, the complex angle γ in equation (20) is the important parameter to tune. Note that there is not much freedom to choose γ as it appears in FIMP DM coupling discussed earlier. We consider it to be vanishingly small for leptogenesis discussions. As in the preceding analysis, here also we consider Thus it is expected that N 2 will dominantly contribute to the baryon asymmetry.
It is to noted that in the present scenario the inflaton has several other decay modes, in addition to its decay into right handed neutrinos. Thus it is difficult to generate the observed amount of baryon asymmetry where the inflaton decays to RH neutrinos are subdominant or Br φ→N 2,3 N 2,3 1. So, one needs to find the parameter space where the branching ratio of inflaton to right handed neutrinos as well as the CP asymmetry from right handed neutrino decay can be large enough to satisfy the requirement of baryon asymmetry. The decay widths of RH neutrinos N 2 and N 3 into SM leptons and Higgs depend on the strength of Yukawa couplings as defined in equation (19). Below we provide the structure of Y D i2 and In Fig. 18, we show the allowed region which satisfies the bound on Y B in M H 2 −T R plane for two different sets of complex angle γ considering M N 2 = 10 9 GeV. We vary g BL and λ 2 in specified ranges mentioned in the figure. The M N 2 < T R and M H 2 < 2M N 2 regions are shown in magenta and yellow colours respectively which are outside the regime of non-thermal leptogenesis discussed here. Similar plot is shown in Fig. 19 considering slightly higher scale of leptogenesis M N 2 = 10 10 GeV where the allowed region gets enhanced, as expected. In preparing both the figures we have taken λ 3 ∼ O(10 −15 ), such that the Br φ→N 2,3 N 2,3 does not turn very small due to other decay modes of inflaton which depend upon λ 3 or scalar mixing.
We have also confirmed that corresponding to our choices of γ, the condition Γ N 2,3 Γ H 2 is satisfied, a requirement for validating the simplistic approach adopted here.

VII. CONCLUSION
To summarise, we have studied the very popular gauged B − L extension of the standard model by restricting ourselves to the minimal possible framework from the requirement of triangle anomaly cancellation, desired gauge symmetry breaking and origin of light neutrino mass. We particularly focus on the possibility of singlet scalar field responsible for breaking B − L gauge symmetry spontaneously to also drive successful inflation in agreement with Planck 2018 data and its implications for dark matter and leptogenesis. While the lightest right handed neutrino is considered to be the DM candidate, the heavier two right handed neutrinos generate light neutrino masses through type I seesaw mechanism and also generate the required lepton asymmetry via their out of equilibrium decay. We first show that the requirement of successful inflationary phase tightly constraints the scalar sector and gauge sector couplings of the model. To be more precise, the stability of the inflationary potential in fact puts an upper bound on B − L gauge coupling along with inflaton couplings to SM Higgs as well as right handed neutrinos. Since WIMP type DM in this model primarily interacts with the SM particles via B − L gauge or singlet scalar (via its mixing with SM Higgs), the bounds derived from inflation on couplings and masses involved in these portals make WIMP annihilations inefficient for most of the parameter space. The parameter space where WIMP abundance satisfies the Planck 2018 data on DM abundance along with inflationary requirements, gets ruled out by LHC data on dilepton searches. This led to our first main conclusion that thermal DM is disfavoured in such scenario. We then considered the possibility of non-thermal DM by considering two different broad scenarios related to the interplay of inflaton mass and reheat temperature. We show that in both the scenarios correct FIMP abundance can be produced. We find that for a scenario where inflaton is not part of the thermal bath after reheating, the required FIMP relic can be produced only if it is allowed to couple to SM leptons opening up several production channels from SM bath.
Such a scenario does not require any additional Z 2 symmetry needed for stabilising WIMP type DM and also have interesting consequences for indirect detection experiments due to possible decays into photons ranging from X-ray to gamma rays.
We then briefly discuss the possibility of leptogenesis by focusing primarily on nonthermal leptogenesis which is very much sensitive to the details of inflation. While resonant leptogenesis is still a viable option, thermal vanilla leptogenesis is not possible due to low reheat temperature predicted in our scenario. We find that inflationary requirements tightly constrain the scenario of non-thermal leptogenesis, precisely due to the same reason behind constraining or disfavouring WIMP type DM mentioned earlier. We show the possibility of producing observed baryon asymmetry from non-thermal leptogenesis for benchmark choices of some parameters while varying others and also show that the same parameters are also consistent with successful inflation, stability of inflaton potential, FIMP DM abundance, neutrino mass apart from other experimental limits. Since the model is very minimal, it remains very predictive, specially when the requirements of correct neutrino mass, DM abundance, baryon asymmetry along with successful inflation are to be met with. Future data from all these frontiers should be able to restrict the model parameters to even stricter ranges while ruling out some of the possibilities.