TeV scale leptogenesis via dark sector scatterings

Abstract

We propose a novel scenario of generating lepton asymmetry via annihilation and coannihilation of dark sector particles including t-channel processes. In order to realistically implement this idea, we consider the scotogenic model having three right handed neutrinos and a new scalar doublet, all of which are odd under an in-built $Z_2$ symmetry. The lightest $Z_2$ odd particle, if electromagnetically neutral, can be a dark matter candidate while annihilation and coannihilation between different $Z_2$ odd particles into standard model leptons serve as the source of lepton asymmetry. The light neutrino masses arise at one loop level with $Z_2$ odd fields going inside the loop. We show that experimental data related to light neutrinos, dark matter relic abundance and baryon asymmetry can be simultaneously satisfied in the model for two different cases: one with fermion dark matter and the other with scalar dark matter. In both scenarios, t-channel annihilation as well as coannihilation of $Z_2$ odd particles play a non-trivial role in producing the non-zero CP asymmetry. Both the scenarios remain allowed from DM direct detection while keeping the scale of leptogenesis as low as TeV or less, lower than the one for vanilla leptogenesis scenario in scotogenic model along with the additional advantage of explaining the baryon-dark matter coincidence to some extent. Due to such low scale, the model is testable through rare decay experiments looking for charged lepton flavour violation.

1 Introduction

There have been significant progress in last few decades in gathering evidences suggesting the presence of a mysterious, non-luminous form of matter, known as dark matter (DM) in the present universe, whose amount is approximately five times more than the ordinary luminous or baryonic matter density $\Omega _B \approx 5\%$ [1]. Among different beyond standard model (BSM) proposals for DM, the weakly interacting massive particle (WIMP) paradigm remains the most widely studied scenario where a DM candidate typically with electroweak (EW) scale mass and interaction rate similar to EW interactions can give rise to the correct DM relic abundance, a remarkable coincidence often referred to as the WIMP Miracle. On the other hand, out of equilibrium decay of a heavy particle leading to the generation of baryon asymmetry has been a very well known mechanism for baryogenesis [2, 3]. One interesting way to implement such a mechanism is leptogenesis [4] where a net leptonic asymmetry is generated first which gets converted into baryon asymmetry through $B+L$ violating EW sphaleron transitions. The interesting feature of this scenario is that the required lepton asymmetry can be generated within the framework of the seesaw mechanism that explains the origin of tiny neutrino masses [5], another observed phenomena which the SM fails to address.

Although these popular scenarios can explain the phenomena of DM and baryon asymmetry independently, it is nevertheless an interesting observation that DM and baryon abundance are very close to each other, within the same order of magnitudes $\Omega _{\mathrm{DM}} \approx 5 \Omega _{B}$. Discarding the possibility of any numerical coincidence, one is left with the task of constructing theories that can relate the origin of these two observed phenomena in a unified manner. There have been several proposals already which mainly fall into two broad categories. In the first one, the usual mechanism for baryogenesis is extended to apply to the dark sector which is also asymmetric [6,7,8,9]. The second one is to produce such asymmetries through annihilations [10,11,12] where one or more particles involved in the annihilations eventually go out of thermal equilibrium in order to generate a net asymmetry. The so-called WIMPy baryogenesis [13,14,15,16,17,18] belongs to this category, where a dark matter particle freezes out to generate its own relic abundance and then an asymmetry in the baryon sector is produced from DM annihilations.

While there is no evidence yet for seesaw mechanism, recently the so-called scotogenic model [19] as an alternative to canonical seesaw mechanism has been extensively studied, where Majorana light neutrino masses can be generated at one loop level with DM particle in the loop. In the scotogenic model, the required lepton asymmetry can be generated through right handed neutrino decays at a low scale $M_N \sim 10$ TeV at the cost of a strongly hierarchical neutrino Yukawa structure [20, 21], but it can not explain the coincidence of baryon asymmetry and DM abundance.

An interesting question raised is whether leptogenesis through (co-)annihilations of $Z_2$ odd particles can be realised by lowering the scale of leptogenesis further compared to vanilla leptogenesis in the scotogenic model. Giving an answer to this question is the main purpose of this work. We examine how the (co-)annihilations of $Z_2$ odd particles can produce the lepton asymmetry while keeping the correct DM abundance, and show that the DM relic abundance is correlated with baryon asymmetry in the scenario. To include all possible annihilations producing lepton asymmetry, we consider the annihilations and coannihilations of all $Z_2$ odd particles instead of restricting them to the lightest $Z_2$ odd particle which is also the DM candidate. If we consider the neutral component of the $Z_2$ odd scalar doublet as DM, then there exists s-channel coannihilation diagrams between DM and right handed neutrinos which can produce a net leptonic asymmetry. While the model satisfies correct DM abundance and lepton asymmetry, the DM sector can be probed at direct detection experiments as well as colliders due to the electroweak gauge interactions of scalar doublet DM. We then consider the fermion DM scenario where the lightest right handed neutrino plays the role of DM. While the annihilation of a pair of the fermionic DM can not produce a net lepton asymmetry in this case, the annihilation of $Z_2$ odd scalar doublets can contribute to the generation of lepton asymmetry. Due to the natural absence of typical s-channel diagrams of scalar doublet annihilations leading to lepton asymmetry, here we show how t-channel diagrams (both tree level and one loop level) can play a non-trivial role in creating the required asymmetry.

As far as we know, the contributions of (co-)annihilations of dark sector particles to lepton asymmetry in this minimal model was not considered before. In both the scenarios we address here, the criteria for “on-shell”-ness of loop particles in one loop annihilation diagrams dictate the particle spectrum and hence the nature of dark matter candidate. We show that it is possible to satisfy the requirement of baryon asymmetry, light neutrino mass and DM related constraints in both the scenarios while keeping the scale of leptogenesis as low as 5 TeV, lower than the scale of vanilla leptogenesis in the same model [21, 22]. Due to such a low scale, the model has another advantage to predict observable rates of charged lepton flavour violation accessible by the sensitivity of the future experiments.

This paper is arranged as follows. In Sect. 2, we briefly review on the minimal scotogenic model followed by detailed discussion on leptogenesis from annihilation and coannihilations in this model in Sect. 3. We finally conclude in Sect. 4.

2 Minimal scotogenic model

The minimal scotogenic model [19] is the extension of the SM by three copies of right handed singlet neutrinos $N_i, i \in 1,2,3$ and one scalar field $\eta$ transforming as a doublet under $SU(2)_L$. An additional discrete symmetry $Z_2$ is incorporated under which these new fields are odd giving rise to the possibility of the lightest $Z_2$-odd particle being a suitable DM candidate. The Lagrangian involving the newly added singlet fermions is

$$\begin{aligned} {{{\mathcal {L}}}} \supset \frac{1}{2}(M_N)_{ij} N_iN_j + \left( y_{ij} \, {\bar{L}}_i {\tilde{\eta }} N_j + \text {h.c.} \right) . \end{aligned}$$

(1)

The electroweak symmetry breaking occurs due to the non-zero vacuum expectation value (VEV) acquired by the neutral component of the SM Higgs doublet while the $Z_2$-odd doublet $\eta$ does not acquire any VEV. After the EWSB these two scalar doublets can be written in the following form in the unitary gauge,

$$\begin{aligned} H=\left( 0, \;\; \frac{ v +h }{\sqrt{2}} \right) ^T , ~~~\eta =\left( \eta ^\pm , \;\; \frac{\eta _R+i\eta _I}{\sqrt{2}}\right) ^T. \end{aligned}$$

(2)

The scalar potential of the model is

$$\begin{aligned} \begin{aligned} V&= \mu _H^2|H|^2 +\mu _{\eta }^2|\eta |^2+ \frac{\lambda _H}{2}|H|^4+ \frac{\lambda _\eta }{2}|\eta |^4 \\&\quad +\lambda _3|H|^2|\eta |^2 +\lambda _4|H^\dag \eta |^2 + \left\{ \frac{\lambda _5}{2}(H^\dag \eta )^2 + \text {h.c.}\right\} . \end{aligned} \end{aligned}$$

(3)

The masses of the physical scalars at tree level can be written as

$$\begin{aligned} m_h^2= & {} \lambda _H v^2 ,\nonumber \\ m_{\eta ^\pm }^2= & {} \mu _\eta ^2 + \frac{1}{2}\lambda _3 v^2 , \nonumber \\ m_{\eta _R}^2= & {} \mu _\eta ^2 + \frac{1}{2}(\lambda _3+\lambda _4+\lambda _5)v^2 \nonumber \\= & {} m^2_{H^\pm }+ \frac{1}{2}\left( \lambda _4+\lambda _5\right) v^2 , \nonumber \\ m_{\eta _I}^2= & {} \mu _\eta ^2 + \frac{1}{2}(\lambda _3+\lambda _4-\lambda _5)v^2 \nonumber \\= & {} m^2_{H^\pm }+ \frac{1}{2}\left( \lambda _4-\lambda _5\right) v^2. \end{aligned}$$

(4)

Here $m_h$, $m_{\eta _R}$, and $m_{\eta _I}$ are the masses of the SM like Higgs boson, the CP even and CP odd scalars from the inert doublet, respectively. $m_{\eta ^\pm }$ is the mass of the charged scalar. Without any loss of generality, we consider $\lambda _5<0, \lambda _4+\lambda _5 <0$ so that the CP even scalar is the lightest $Z_2$ odd particle and hence a stable dark matter candidate.

Denoting the squared physical masses of neutral scalar and pseudo-scalar parts of $\eta$ as $m^2_{R,I}=m^2_{\eta _R,\eta _I}$ and the mass of the right handed neutrino $N_k$ in the internal line as $M_k$, the one loop neutrino mass can be estimated as [19]

$$\begin{aligned} (M_{\nu })_{ij} \,&= \, \sum _k \frac{y_{ik}y_{jk} M_{k}}{32 \pi ^2} \nonumber \\&\quad \times \left( \frac{m^2_{\eta _R}}{m^2_{\eta _R}-M^2_k} \, \text {ln} \frac{m^2_{\eta _R}}{M^2_k}-\frac{m^2_{\eta _I}}{m^2_{\eta _I}-M^2_k}\, \text {ln} \frac{m^2_{\eta _I}}{M^2_k} \right) \nonumber \\&\equiv \ \sum _k \frac{y_{ik}y_{jk} M_{k}}{32 \pi ^2} \left[ L_k(m^2_{\eta _R})-L_k(m^2_{\eta _I})\right] \, , \nonumber \\ \Lambda _k&= \frac{ M_{k}}{32 \pi ^2} \left[ L_k(m^2_{\eta _R})-L_k(m^2_{\eta _I})\right] , \end{aligned}$$

(5)

where $M_k$ is the mass eigenvalue of the right handed neutrino mass eigenstate $N_k$ in the internal line and the indices $i, j = 1,2,3$ run over the three neutrino generations. The function $L_k(m^2)$ is defined as

$$\begin{aligned} L_k(m^2) \ = \ \frac{m^2}{m^2-M^2_k} \, \text {ln} \frac{m^2}{M^2_k} \, . \end{aligned}$$

(6)

From the physical scalar masses given above, we note that $m^2_{\eta _R}-m^2_{\eta _I}=\lambda _5 v^2$. In this model for the neutrino mass to match with experimentally observed limits ($\sim 0.1$ eV), Yukawa couplings of the order $10^{-3}$ are required if $M_k$ is as low as 1 TeV and the mass difference between $\eta _R$ and $\eta _I$ is kept around 1 GeV. Such a small mass splitting between $\eta _R$ and $\eta _I$ will correspond to small quartic coupling $\lambda _5 \sim 10^{-4}$. Thus, one can suitably choose the Yukawa couplings, quartic coupling $\lambda _5$ and $M_k$ in order to arrive at sub eV light neutrino masses. To be in exact agreement with light neutrino masses, we first rewrite the neutrino mass given above in Eq. (5) in the form of type I seesaw formula

$$\begin{aligned} M_\nu \ = \ y {\widetilde{M}}^{-1} y^T \, , \end{aligned}$$

(7)

where we have introduced the diagonal matrix ${\widetilde{M}}$ with elements

$$\begin{aligned} {\widetilde{M}}_i \&= \ \frac{2\pi ^2}{\lambda _5}\zeta _i\frac{2M_i}{v^2} \, , \end{aligned}$$

(8)

$$\begin{aligned} \text {and}\quad \zeta _i&\ = \ \left( \frac{M_{i}^2}{8(m_{\eta _R}^2-m_{\eta _I}^2)}\left[ L_i(m_{\eta _R}^2)-L_i(m_{\eta _I}^2) \right] \right) ^{-1} \, . \end{aligned}$$

(9)

The light neutrino mass matrix (7) is diagonalised by the usual Pontecorvo–Maki–Nakagawa–Sakata (PMNS) mixing matrix U, which is determined from the neutrino oscillation data (up to the Majorana phases):

$$\begin{aligned} D_\nu \ = \ U^\dag M_\nu U^* \ = \ \text {diag}(m_1,m_2,m_3) \, . \end{aligned}$$

(10)

Then the Yukawa coupling matrix satisfying the neutrino data can be written as

$$\begin{aligned} y \ = \ U D_\nu ^{1/2} O {\widetilde{M}}^{1/2} \, , \end{aligned}$$

(11)

where O is an arbitrary complex orthogonal matrix. This is the equivalent of the Casas–Ibarra parametrisation [23] for scotogenic model [24].

3 Leptogenesis from annihilations

In the minimal scotogenic model discussed in the previous section, there are different types of annihilation processes which violate lepton number. They are namely,

• 1. annihilation process of scalar doublet $\eta$: $\eta \eta \rightarrow L_\alpha L_\beta$.

• 2. coannihilation process of scalar doublet and one of the singlet fermions: $\eta N \rightarrow L X$ where ($X \equiv h,\gamma ,W^{\pm },Z$).

Interestingly, if we put the additional constraints that such lepton number violating annihilations and coannihilations also generate a non-zero CP asymmetry, they lead to two different DM possibilities namely,

• 1. the lightest neutral component of inert scalar doublet $\eta$ as DM,

• 2. the lightest right handed neutrino N as DM.

The Boltzmann equations for leptonic asymmetry is given as follows:

$$\begin{aligned} \frac{dY_{\Delta L}}{dz}&= \frac{1}{zH(z)} \times \left[ \sum _i \epsilon _{N_i}(Y_{N_i} - Y^{\mathrm{eq}}_{N_i})\langle \Gamma _{N_i \rightarrow L_\alpha \eta }\rangle \right. \nonumber \\&\quad - Y_{\Delta L}r_{N_i}\langle \Gamma _{N_i \rightarrow L_\alpha \eta }\rangle -Y_{\Delta L}r_\eta s\langle \Gamma _{\eta \rightarrow N_1 L}\rangle \nonumber \\&\quad + 2\epsilon _{\eta \eta } s\langle \sigma v\rangle _{\eta \eta \rightarrow L L}\left( Y^2_{\eta } - (Y_{\eta }^{\mathrm{eq}})^2\right) \nonumber \\&\quad -Y_{\Delta L}Y^{\mathrm{eq}}_L r^2_\eta s\langle \sigma v\rangle _{\eta \eta \rightarrow LL} \nonumber \\&\quad + \sum _i \epsilon _{N_i \eta } s\langle \sigma v\rangle _{\eta N_i\rightarrow L \mathrm{SM}}\left( Y_{\eta }Y_{N_i} - Y^{\mathrm{eq}}_{\eta }Y^{\mathrm{eq}}_{N_i}\right) \nonumber \\&\quad - \frac{1}{2}Y_{\Delta L}Y^{\mathrm{eq}}_L r_{N_i}r_\eta s\langle \sigma v\rangle _{\eta N_i \rightarrow \mathrm{SM} {\overline{L}}} \nonumber \\&\quad -Y_{\Delta L}Y^{\mathrm{eq}}_\eta s\langle \sigma v \rangle ^{wo}_{\eta L \rightarrow \eta {\overline{L}}} \nonumber \\&\quad -\sum _iY_{\Delta L}Y^{\mathrm{eq}}_\eta s\langle \sigma v \rangle ^{wo}_{\eta L \rightarrow N_iX} \nonumber \\&\quad \left. -\sum _iY_{\Delta L}Y^{\mathrm{eq}}_{N_i} s\langle \sigma v \rangle ^{wo}_{N_i L \rightarrow \eta X} \right] , \nonumber \\ H(z)&= \sqrt{\frac{4\pi ^3 g_*}{45}}\frac{M^2_{\mathrm{DM}}}{z^2M_{\mathrm{PL}}}, \quad s = g_* \frac{2\pi ^2}{45}\left( \frac{M_{\mathrm{DM}}}{z}\right) ^3, \nonumber \\ r_j&= \frac{Y^{\mathrm{eq}}_j}{Y^{\mathrm{eq}}_L}, \quad \quad \langle \Gamma _{j\rightarrow X} \rangle = \frac{K_1(M_j/T)}{K_2(M_j/T)}\Gamma _{j\rightarrow X}, \end{aligned}$$

(12)

where $z=\frac{M_{\mathrm{DM}}}{T}$, $M_{\mathrm{PL}}$ is the Planck mass and $Y=n/s$ denotes the comoving number density, as the ratio of number density to entropy density. The details of the derivation of this Boltzmann equation as well as the relevant equations for DM are presented in appendix A. In the above equation, $\epsilon _{\eta \eta }$ and $\epsilon _{N_i \eta }$ will be given appropriately later and $\epsilon _{N_i}$ is taken from [21].

$$\begin{aligned} \epsilon _{N_i}&= \sum _\alpha \frac{1}{8\pi (y^\dagger y)_{ii}}\sum _{j\ne i} \nonumber \\&\quad \times \left[ f\left( \frac{M^2_{N_j}}{M^2_{N_i}}, \frac{m^2_\eta }{M^2_{N_i}}\right) \mathfrak {I}[y^*_{\alpha ,i}y_{\alpha ,j}(y^\dagger y)_{ij}] \right. \nonumber \\&\quad -\left. \frac{M^2_{N_i}}{M^2_{N_j}-M^2_{N_i}} \left( 1-\frac{m^2_\eta }{M^2_{N_i}}\right) ^2 \mathfrak {I}[y^*_{\alpha i}y^*_{\alpha j}H_{ij}]\right] \nonumber \\ H_{ij}&= (y^\dagger y)_{ij}\frac{M_{N_j}}{M_{N_i}} + (y^\dagger y)^*_{ij}; \nonumber \\ f(x,w)&= \sqrt{x}\left[ 1+\frac{1-2w+x}{(1-w)^2} \ln \left( \frac{x-w^2}{1-2w+x}\right) \right] \end{aligned}$$

(13)

$M_{DM}$ in s is the mass of $Z_2$ odd particle taken appropriately depending on the scenario mentioned above. Here $K_n$ is the nth order Modified Bessel function of second kind. The details of the Boltzmann equations are given in appendix A. We will now discuss this general framework in the context of two specific scenarios of DM mentioned above in the upcoming subsections.

The above Boltzmann equation contains the next leading order (NLO) contributions to lepton asymmetry. For usual type I seesaw leptogenesis, such NLO effects have been calculated already, for example, see [25] and references therein. It is therefore necessary to include or consider all the diagrams which can either contribute to lepton asymmetry or washout at the same order of couplings. We first only show the annihilation of dark sector particles into SM ones in Fig. 1 for the case of scalar doublet dark matter. We do not present the usual well-known two body decay diagrams (both tree and one loop) for simplicity. Based on the diagrams shown in Fig. 1, one can construct several other diagrams just by interchanging initial and final state particles. For example, it is straightforward to consider three body decay diagram of $N_i$ into $\eta , L, X, X \equiv \gamma , W^{\pm }, Z$. This diagram, which contributes to the asymmetry at order ${\mathcal {O}}(y^4g^2)$ will be suppressed due to not only phase space but also additional ${\mathcal {O}}(g^2)$ suppression.

Now, although the 3-body decay and the co-annihilation are of the same order, the contribution in the Boltzmann evolution is different as they enter the Boltzmann equation differently. The 3-body decay will enter into the equation as an addition to the 2-body decay in which the contribution coming from the 2 body decay will dominate. As will be shown later, the contributions of coannihilation to lepton asymmetry may dominant over that from decay process. This is due to fact that the imaginary part of the interference is not suppressed as it becomes two $2\rightarrow 2$ processes (referring to the diagrams in the second line in Fig. 2 leading to the interference at order ${\mathcal {O}}(y^4g^2)$) both mediated by the leptons after cutting loop diagrams, whereas in case of the decay the loop diagram becomes $1\rightarrow 2$ and $2\rightarrow 2$ in which the $2\rightarrow 2$ is a t-channel process mediated by the right-handed neutrino after cutting. We present all such 1-loop diagrams contributing to the asymmetry arising from the interference at order ${\mathcal {O}}(y^4g^2)$ and ${\mathcal {O}}(y^6)$ in Fig. 2. Similarly, there exists several washout processes that can be constructed by swapping initial and final state particles. This is true for fermion dark matter as well, which we discuss in one of the upcoming sections. To point out the significance of such processes at the same NLO, we therefore make a comparison of different washout processes (both inverse decay and scatterings) and show them in in Fig. 3 (lightest scalar as Dark Matter) and in Fig. 6 (lightest fermion as Dark Matter).

Fig. 1

Feynman diagrams contributing to $\langle \sigma v\rangle _{\mathrm{DM DM} \rightarrow X L}$ and the asymmetry $\epsilon$ at leading order. Here $X \equiv \gamma ,W^{\pm },Z$

Fig. 2

1-loop diagrams contributing to the asymmetry arising from the interference at ${\mathcal {O}}(y^4g^2)$ and ${\mathcal {O}}(y^6)$

3.1 Scalar doublet $\eta$ as dark matter

In this scenario the only way to get asymmetry is through co-annihilations. Pure scalar annihilations give rise to vanishing leptonic asymmetry if $\eta$ is the lightest $Z_2$ odd particle, as required for it to be the DM candidate. This is particularly due to the fact the “on-shell” criteria of loop particles can not be realised in such case, resulting in a vanishing CP asymmetry, as we discuss below.

The relevant coannihilation processes, both tree level as well as one loop level, are shown in Figs. 1 and in 2. One may notice that the one loop self-energy diagrams, arising from the lepton propagator, do not contribute to the lepton asymmetry because the processes occur before electroweak symmetry breaking and thus lepton mass should be zero giving rise to vanishing self-energy loop contribution. Thus, only interference between tree and one loop vertex correction can give rise to CP asymmetry. For this scenario the Boltzmann equations for the $Z_2$ odd particles take the following form:

$$\begin{aligned} \frac{dY_{N_k}}{dz}&= -\frac{1}{zH(z)} \times \left[ (Y_{N_k}-Y^{\mathrm{eq}}_{N_k})\langle \Gamma _{N_k \rightarrow L_\alpha \eta }\rangle \right. \nonumber \\&\quad \left. + (Y_{N_k}Y_{\eta }-Y^{\mathrm{eq}}_{N_k}Y^{\mathrm{eq}}_{\eta })s\langle \sigma v \rangle _{\eta N_k\rightarrow L\mathrm{SM}} \right. \nonumber \\&\quad + \left. \sum _{l=1}^3(Y_{N_k}Y_{N_l}-Y^{\mathrm{eq}}_{N_k}Y^{\mathrm{eq}}_{N_l})s\langle \sigma v \rangle _{N_l N_k\rightarrow \mathrm{SM SM}}\right] , \nonumber \\ \frac{dY_{\eta }}{dz}&= \frac{1}{zH(z)} \times \left[ (Y_{N_k}-Y^{\mathrm{eq}}_{N_k})\langle \Gamma _{N_i \rightarrow L_\alpha \eta }\rangle \right. \nonumber \\&\quad - 2(Y^2_{\eta }-(Y^{\mathrm{eq}}_{\eta })^2)s\langle \sigma v\rangle _{\eta \eta \rightarrow \mathrm{SM SM}} \nonumber \\&\quad - \left. \sum ^3_{m=1}(Y_{N_m}Y_{\eta }-Y^{\mathrm{eq}}_{N_m}Y^{\mathrm{eq}}_{\eta })s\langle \sigma v \rangle _{\eta N_m\rightarrow L\mathrm{SM}} \right] . \end{aligned}$$

(14)

The CP asymmetry arising from the interference between tree and 1-loop diagrams in Figs. 1 and 2 can be estimated as

$$\begin{aligned} \epsilon _{N_i \eta }&= \frac{2}{(yy^\dagger )_{ii}} \sum _{j}\mathfrak {I}[(yy^\dagger )^2_{ij}]{\widetilde{\epsilon }}_{ij}, \end{aligned}$$

(15)

$$\begin{aligned} {\widetilde{\epsilon }}_{ij}&= \frac{\sqrt{x_j}}{6 x_i \left( -x_i^{3/2}+x_i (x_j-2)+\sqrt{x_i} x_j+1\right) ^2(\sqrt{x_i}-3)} \nonumber \\&\quad \times \left( x_i^{7/2} (3 x_j+1) +\sqrt{x_i} (3 x_j+5)+1 \right. \nonumber \\&\quad - 3 x_i^{5/2} \left( x_j \left( D+(x_j-3) x_j+4\right) - 3 D -2\right) \nonumber \\&\quad -3 x_i^{3/2} \left( 2 \left( D{+}3\right) {+} x_j \left( x_j \left( D{+}x_j{+}1\right) -D-4\right) \right) \nonumber \\&\quad - x_i^4+f^3 \left( 3 D+3 x_j^2+11\right) \nonumber \\&\quad - 3 x_i^2 \left( x_j \left( D+2 (x_j-1) x_j+2\right) - D +6\right) \nonumber \\&\quad \left. +x_i \left( 1-3 x_j \left( D+x_j-4\right) \right) \right) \nonumber \\&\quad + \frac{\sqrt{x_j}}{4 x_i}\left( \sqrt{x_i}-1+ \frac{\sqrt{x_j}}{(1+\sqrt{x_i})^2}(\sqrt{x_i}-1+r_j)\right. \nonumber \\&\quad \times \left( \log \left( \frac{1+\sqrt{x_i}x_j}{x_i(1+\sqrt{x_i})} \right) \right. \nonumber \\&\quad -\log \left( \frac{1+x_i+x^{3/2}_i+\sqrt{x_i}x_j}{x_i(1+\sqrt{x_i})} \right) \nonumber \\&\quad + \left. \left. \log \left( 1+\frac{1+\sqrt{x_i}}{\sqrt{x_i} (\sqrt{x_i}-1+x_i+x_j)}\right) \right) \right) \nonumber \\&\quad + \frac{\sqrt{x_i} x_j{\tilde{\Gamma }}_j}{\pi ((x_i-x_j)^2 + x_j{\tilde{\Gamma }}^2_j)} \nonumber \\ D&= \sqrt{(x_i-x_j) \left( x_i+4 \sqrt{x_i}-x_j+4\right) } \nonumber \\ x_l&= \frac{M^2_{N_l}}{m^2_\eta } \qquad {\tilde{\Gamma }}_j = \frac{\Gamma _j}{m_\eta }. \end{aligned}$$

(16)

where the details of the asymmetry is shown in appendix B. Although the above expression is an $s-$ wave approximation for actual expression shown in appendix B, we have used the actual expression for our analysis. It should be noted that in the above expression always ($1\le x_j\le x_i$) where j stands for $N_j$ inside the loop while i stands for $N_i$ as one of the initial state particles, shown in Figs. 1 and 2. This is simply to realise the “on-shell” -ness of the loop particles in order to generate the required CP asymmetry.

There are several wash-out processes in this scenario, categorised as follows:

• $\mathbf \Delta L = 2$: $L \eta \rightarrow {\overline{L}} \eta , \eta \eta \rightarrow L L$ are purely wash-out processes.

• $\mathbf \Delta L = 1$: there are two main sources of such wash-out, namely

  • 1. inverse decay of $\Gamma (N_k \rightarrow L \eta )$,

  • 2. inverse process of co-annihilation $N_k \eta \rightarrow L,X(=\gamma , W, Z, h)$.

We have taken them into account in our numerical calculations.

Adopting the Casas–Ibarra parametrisation given in Eq. (11) , we see that CP phases in U do not contribute to $\epsilon _{N_i\eta }$, but complex variables in the orthogonal matrix O can lead to non-vanishing value of $\epsilon _{N_i\eta }$. This is similar to leptogenesis from pure decay in this model [21] where, in the absence of flavour effects, the orthogonal matrix O played a crucial role. In general, this $3\times 3$ orthogonal matrix O can be parametrised by three complex parameters of type $\theta _{\alpha \beta } = \theta ^R_{\alpha \beta } + i\theta ^I_{\alpha \beta }, \theta ^R_{\alpha \beta } \in [0, 2\pi ], \theta ^I_{\alpha \beta } \in {\mathbb {R}}$ [26].¹ In general, the orthogonal matrix O for n flavours can be product of $^nC_2$ number of rotation matrices of type

$$\begin{aligned} O_{\alpha \beta }&= \begin{pmatrix} \cos {(\theta ^R_{\alpha \beta } + i\theta ^I_{\alpha \beta })} &{} \cdots &{} \sin {(\theta ^R_{\alpha \beta } + i\theta ^I_{\alpha \beta })} \\ \vdots &{} \ddots &{} \vdots \\ - \sin {(\theta ^R_{\alpha \beta } + i\theta ^I_{\alpha \beta })} &{} \cdots &{} \cos {(\theta ^R_{\alpha \beta } + i\theta ^I_{\alpha \beta })} \end{pmatrix}, \end{aligned}$$

(17)

with rotation in the $\alpha -\beta$ plane and dots stand for zero. For example, taking $\alpha =1, \beta =2$ we have

$$\begin{aligned} O_{12}&= \begin{pmatrix} \cos {(\theta ^R_{12} + i\theta ^I_{12})} &{} \sin {(\theta ^R_{12} + i\theta ^I_{12})} &{} 0 \\ -\sin {(\theta ^R_{12} + i\theta ^I_{12})} &{} \cos {(\theta ^R_{12} + i\theta ^I_{12})} &{} 0 \\ 0 &{} 0 &{} 1\end{pmatrix}. \end{aligned}$$

(18)

The above asymmetry along with this rotation (one at a time) takes the following form

$$\begin{aligned} \epsilon _{N_i\eta }&= \sum _j \times \frac{(m^4_j-m^4_i)\sin {(2\theta ^R_{ij})} \sinh {(2\theta ^I_{ij})}}{\Lambda _j^2((m^2_i-m^2_j) \cos (2\theta ^R_{ij})+(m^2_i+m^2_j) \cosh (2\theta ^I_{ij}))} \nonumber \\&\quad \times {\widetilde{\epsilon }}_{ij}, \end{aligned}$$

(19)

where $m_i$’s are the light neutrino masses and $\Lambda _i$’s are defined above in Eq. (5). On the right hand side of the above equation, a summation over index j is implicit. As an example, we have taken the benchmark values shown in Table 1 to compute the baryon asymmetry as well as scalar DM relic numerically. Here we consider $\eta _I$ as the DM candidate ($\eta _I \equiv \mathrm DM$, corresponding to positive value of quartic coupling $\lambda _5$) which is similar to the inert doublet model discussed extensively in the literature [19, 28, 29]. Typically there exists two distinct mass regions, $M_{\mathrm{DM}} \le 80$ GeV and $M_{\mathrm{DM}} \ge 500$ GeV, where correct relic abundance criteria can be satisfied. In both regions, depending on the mass differences $m_{\eta ^\pm }-m_{\eta _I}, m_{\eta _R}-m_{\eta _I}$, the coannihilations of $\eta _I, \eta ^\pm$ and $\eta _R, \eta _I$ can also contribute to the DM relic abundance [30, 31]. As for the mixing angles in the PMNS matrix U we took the best fit values obtained from the recent global fit analysis [32] shown in Table 2.

Table 1 The numerical values of the parameters chosen for generating correct scalar DM relic and baryon asymmetry. We denote it as benchmark point 1 (BP1). The fraction $\omega$ is to ensure to get the correct baryonic asymmetry. For Normal (Inverted) Hierarchy $\omega = 0.7(0.56)$

Table 2 The numerical values of light neutrino parameters used in the calculations

To perform the numerical analysis, we implement the model in SARAH 4 [33] and extract the thermally averaged annihilation rates from micrOMEGAs 4.3 [34] to use while solving the Boltzmann equations above. In Fig. 3, we plot the comoving number densities of all $Z_2$ odd particles, taking part in generating the lepton asymmetry along with the generated asymmetry $\Delta L$, as functions of temperature. The left panel corresponds to normal hierachy (NH) of neutrino mass spectrum and the right panel to the case of inverted hierarchy (IH). The horizontal solid black line labelled as “$\Delta L$ observed” correspond to the value of $\Delta L$ that is partially converted into the observed baryon asymmetry via the electroweak sphaleron processes with the conversion factor $C_s = \frac{8N_f + 4N_H}{22N_f+13N_H}$ where $N_f=3, N_H=2$ are the number of fermion generations and Higgs doublets respectively [35]. While sphalerons violate $B+L$, they conserve $B-L$ symmetry. The sphaleron processes are effective with a thermal rate $\Gamma _{\mathrm{sph}}\sim (\alpha _2 T)^4$ with SU(2) gauge coupling constant $\alpha _2$ at high temperature until EW phase transition, while they are exponentially suppressed due to finite gauge boson masses after EW gauge symmetry breaking. The shaded regions in the panels correspond to the temperature below which the sphaleron processes become inoperative($T\lesssim 200$ GeV [36]). The benchmark parameters are chosen in such a way that the generated lepton asymmetry $\Delta L$ by the epoch of sphaleron freeze-out is sufficient enough to produce the observed baryon asymmetry. The dashed horizontal black line corresponds to the observed DM relic abundance in the present universe [1]. As can be seen from this plot, the lepton asymmetry grows as the temperature cools down due to the contributions from the co-annihilation diagrams. While the lepton asymmetry gets converted into the baryon asymmetry at EW phase transition temperature, it takes a while for DM to freeze-out. Since non-zero CP asymmetry arises from coannihilations between $\eta$ and heavier right handed neutrinos $N_{2,3}$, the lepton number generating processes get frozen out much earlier compared to DM self annihilations. This makes sure that a net lepton asymmetry is created with the right amount without being washed out entirely. The limits on the parameter $\lambda _5$, keeping all the other parameters same as shown on Table 1, are set by two constraints. One is the mass difference between the neutral scalar and pseudo-scalar ($\Delta = m_{\eta _R}-m_{\eta _I}$) which needs to be more than approximately 100 keV, in order to avoid Z mediated inelastic direct detection scattering of DM off nucleons, as we discuss below. The other constraint is coming from the required leptonic asymmetry with maximal CP asymmetry. The range is shown in Table 3.

Table 3 The range of $\lambda _5$ allowed by phenomenological requirements of satisfying direct detection bounds and generating the required lepton asymmetry

Fig. 3

Top panel: comoving number densities of $Z_2$ odd particles and lepton asymmetry as a function of temperature for scalar DM scenario. The left (right) panel corresponds to the normal (inverted) hierarchy of neutrino mass spectrum. Bottom panel: rates of individual wash-out processes along with the Hubble rate (left panel); individual (co)annihilation contributions of each channel to the asymmetry. The solid (dashed) black line corresponds to the baryon asymmetry (DM abundance) observed at present epoch ($T \sim 0$). Shaded regions represent the epochs after sphaleron freeze out

Since the DM in this scenario has electroweak gauge interactions, we include DM direct detection constraints arising from tree level Z boson mediated processes $\eta _R n \rightarrow \eta _I n$, n being a nucleon (shown in Fig. 4). The cross section for Z mediated inelastic process is [38]

$$\begin{aligned} \sigma ^0_n&= \frac{G^2_F}{2\pi }\mu ^2_n \simeq 7.44\times 10^{-39} \text {cm}^2. \end{aligned}$$

(20)

We note that one can forbid such scattering if $\delta = m_{\eta _I} - m_{\eta _R} > 100$ keV. Using the expressions for physical masses above, it leads to a lower bound on the dimensionless quartic coupling $\lambda _5$ as

$$\begin{aligned} \lambda _5 \approx 1.65 \times 10^{-7} \left( \frac{\delta }{100 \; \mathrm keV} \right) \left( \frac{M_{\mathrm{DM}}}{100 \; \mathrm GeV} \right) . \end{aligned}$$

This lower limit on $\lambda _5$ becomes weaker for heavier DM masses. Now, if we avoid this bound i.e the mass difference between the scalar and pseudo-scalar dark matter is above 100 keV then the direct detection is majorly dominated by the first process shown in Fig. 4. Our benchmark point satisfies this constraint, forbidding such inelastic DM scattering. The Yukawa couplings for all three generation of leptons for the chosen benchmark point (considering only 13 rotation in O matrix, however) in this scenario can be written in terms of the following matrix.

$$\begin{aligned} y |_{NH}&= \begin{pmatrix} 4.47 - 5.04i &{} 10.7 + 4.77 i &{} 6.35 + 7.96i \\ -0.66 - 5.30i &{} 8.38 - 4.06i &{} 7.29 + 0.59i \\ -6.34 - 2.74i &{} 1.15 - 11.7i &{} 6.20 - 8.2i \end{pmatrix} \nonumber \\&\quad \times 10^{-3}, \end{aligned}$$

(21)

$$\begin{aligned} y |_{IH}&= \begin{pmatrix} 8.52 - 13.8i &{} 5.78 + 0.25i &{} -5.66 - 0.20i\\ -4.89 - 12.8i &{} 4.54 - 3.93i &{} -4.42 + 4.03i \\ -17.18 - 3.51i &{} -2.6 - 6.78i &{} 2.71 + 6.60i \end{pmatrix} \nonumber \\&\quad \times 10^{-3}. \end{aligned}$$

(22)

Now, we move onto the discussion of fermion singlet DM where the direct detection constraints are less severe due to its gauge singlet nature. This is the topic of our next subsection.

Fig. 4

Feynman diagram for inelastic scattering of scalar DM off nucleon mediated by Z boson

3.2 Right handed neutrino as dark matter

If the lightest $Z_2$ odd particle is the lightest of the right handed neutrinos (and hence the DM candidate), then the annihilation processes responsible for creating a non-zero lepton asymmetry are shown in Fig. 5. Once again, pure self annihilation of DM can not provide the asymmetry due to the absence of “on-shell” condition for loop particles. If the scalar doublet $\eta$ is the next to the lightest $Z_2$ odd particle, then the annihilation processes shown in Fig. 5 can produce the required lepton asymmetry. For this scenario the Boltzmann equations for the $Z_2$ odd particles take the following form:

$$\begin{aligned} \frac{dY_{N_k}}{dz}&= -\frac{1}{zH(z)} \times \left[ (Y_{N_k}-Y^{\mathrm{eq}}_{N_k})\langle \Gamma _{N_k \rightarrow L_\alpha \eta }\rangle \right. \nonumber \\&\quad + (Y_{N_k}Y_{\eta }-Y^{\mathrm{eq}}_{N_k}Y^{\mathrm{eq}}_{\eta })s\langle \sigma v \rangle _{\eta N_k\rightarrow L\mathrm{SM}} \nonumber \\&\quad \left. + \sum _{l=1}^3(Y_{N_k}Y_{N_l}-Y^{\mathrm{eq}}_{N_k}Y^{\mathrm{eq}}_{N_l})s\langle \sigma v \rangle _{N_l N_k\rightarrow \mathrm{SM SM}}\right] , \nonumber \\&\quad \quad \text { for}~ k = 2, 3 \nonumber \\ \frac{dY_{\eta }}{dz}&= \frac{1}{zH(z)} \times \left[ \sum ^3_{k=2}(Y_{N_k}-Y^{\mathrm{eq}}_{N_k})\langle \Gamma _{N_k \rightarrow L_\alpha \eta }\rangle \right. \nonumber \\&\quad - (Y_{\eta }-Y^{\mathrm{eq}}_{\eta })\langle \Gamma _{\eta \rightarrow L_\alpha N_1}\rangle \nonumber \\&\quad - 2(Y^2_{\eta }-(Y^{\mathrm{eq}}_{\eta })^2)s\langle \sigma v\rangle _{\eta \eta \rightarrow \mathrm{SM SM}} \nonumber \\&\quad \left. - \sum ^3_{m=1}(Y_{N_m}Y_{\eta }-Y^{\mathrm{eq}}_{N_m}Y^{\mathrm{eq}}_{\eta })s\langle \sigma v \rangle _{\eta N_m\rightarrow L\mathrm{SM}} \right] , \nonumber \\ \frac{dY_{N_1}}{dz}&= \frac{1}{zH(z)} \times \left[ (Y_{\eta }-Y^{\mathrm{eq}}_{\eta })\langle \Gamma _{\eta \rightarrow L_\alpha N_1}\rangle \right. \nonumber \\&\quad - (Y_{N_1}Y_{\eta }-Y^{\mathrm{eq}}_{N_1}Y^{\mathrm{eq}}_{\eta })s\langle \sigma v \rangle _{\eta N_1\rightarrow L\mathrm{SM}} \nonumber \\&\quad \left. - \sum ^3_{l=1}(Y_{N_1}Y_{N_l}-Y^{\mathrm{eq}}_{N_1}Y^{\mathrm{eq}}_{N_l})s\langle \sigma v \rangle _{N_l N_1\rightarrow \mathrm{SM SM}}\right] . \end{aligned}$$

(23)

Fig. 5

Feynman diagrams of the annihilation processes responsible for creating lepton asymmetry in fermion DM scenario

Now, in this scenario along with the co-annihilation channels discussed earlier, the annihilation channels shown in Fig. 5 also contribute to the asymmetry. The CP asymmetry coming from the interference of the tree (Fig. 5) and the loop (bottom two diagrams in Fig. 2) leading to ${\mathcal {O}}(y^6)$ are given as:

$$\begin{aligned} \epsilon _{\eta \eta }&= 8\sum _{ij} \left( \mathfrak {I}[(yy^\dagger )_{i1} (yy^\dagger )_{j1}(yy^\dagger )_{ij}]\right) \epsilon ^{ann}_{ij} \nonumber \\&= \sum _{j}\frac{m^4_1-m^4_j}{\Lambda ^4_1\Lambda ^2_j} \nonumber \\&\quad \times \left( \Delta m^2_{1j}\sin (4\theta ^R_{1j}) + (m^2_1 + m^2_j)\sinh (4\theta ^I_{1j})\right) \epsilon ^{ann}_{1j}, \nonumber \\ \epsilon ^{ann}_{1j}&= \frac{1}{16\pi }\left[ 1-r_1 - \frac{1}{2}(r_1-3)\ln \left[ \frac{1+r_1}{3-r_1}\right] \right] \nonumber \\&\quad \times \frac{\sqrt{r_jr_1}}{(1+r_j)(1+r_1)}\frac{1}{{\mathcal {M}}_{tree}}, \end{aligned}$$

(24)

$$\begin{aligned} {\mathcal {M}}_{tree}&= \sum _{ij} (yy^\dagger )^2_{ij}\frac{\sqrt{r_i}\sqrt{r_j}}{(1+r_i)1+r_j)}. \end{aligned}$$

(25)

At this point one may notice that the asymmetry expression given in Eq. (16) is much more complicated than the one given in Eq. (24) which is purely due to the presence of several coannihilating particles. As pointed out earlier, in this scenario we would require at least one of the right handed neutrinos to be lighter than the scalar doublet $\eta$ whose annihilations are responsible for creating the asymmetry. In our case we have considered only $N_1$ to be lighter than $\eta$ and the rest to be heavier. Alternatively, if $N_{2,3}$ are lighter than $\eta$, then their (co)annihilations can contribute more to the generation to the asymmetry. For our chosen benchmark points, as shown in Table 4, the contribution of $N_{2,3}$ annihilations to lepton asymmetry is sub-dominant compared to $\eta$ annihilations as well as $\eta -N_k\; (k=2,3)$ coannihilations. In fact in the fermion DM scenario, both $\eta -N_k\; (k=2,3)$ coannihilations (shown in Fig. 1) as well as $\eta$ annihilations (shown in Fig. 5) can contribute to lepton asymmetry. It is worthwhile to note that the $\eta$ annihilations shown in Fig. 5 can not contribute to lepton asymmetry in the scalar DM scenario due to the absence of “on-shell” condition for loop particles.

The washout effects in this scenario are categorised as follows :

• $\mathbf \Delta L = 2$: there are two processes of this type:

  • 1. $\eta \eta \rightarrow L L$ and $L \eta \rightarrow {\overline{L}} \eta$ where the former is also responsible for the source of asymmetry while the latter is purely wash-out.

• $\mathbf \Delta L = 1$: there are two main sources of such washout processes:

  • 1. the inverse decay of $\Gamma (N_k \rightarrow L \eta )$ and $\Gamma (\eta \rightarrow L N_k)$, and the second one is purely wash-out not contributing to the asymmetry.

  • 2. the inverse process of co-annihilation $N_k \eta \rightarrow L,X(=\gamma ,W,h)$

Table 4 The numerical values of the parameters chosen for generating correct Fermion DM relic and baryon asymmetry. We denote it as benchmark point 2 (BP2). The fraction $\omega$ is to ensure to get the correct baryonic asymmetry. For Normal (Inverted) Hierarchy $\omega = 0.85(0.45)$

In Fig. 6, we plot the predictions for relic abundance of all $Z_2$ odd particles, taking part in generating the lepton asymmetry along with the generated asymmetry $\Delta L$, as functions of temperature, for fermion DM scenario. The upper left panel corresponds to the NH of neutrino masses, whereas the upper right panel to the IH of neutrino masses. In the lower panel, we compare the relative contributions to the lepton asymmetry: from $\eta -\eta$ annihilations $(\Delta L_{\eta })$, from $\eta -N_k$ coannihilations $(\Delta L_{N \eta })$, from $N_2$ decay $(\Delta L_N)$, as functions of temperature. Similar to the scenario of scalar DM, here also the lepton asymmetry grows as temperature cools due to contributions from the annihilations and coannihilations and saturate around the temperature where the processes responsible for creating the asymmetry tend to go out of equilibrium. Also, the lightest right handed neutrino freezes out to give the required relic abundance for DM in the universe. The shaded regions in the panels correspond to the temperature below which the sphaleron processes become inoperative.

Table 5 The range of $\lambda _5$ allowed by phenomenological requirements of satisfying charged lepton flavour violation bounds and generating the required lepton asymmetry

The first bump in the curve for the lepton asymmetry shown in upper panel plots of Fig. 6 arises from the coannihilation diagrams while the decay contribution enters later and further increases the asymmetry at later stages. Here we notice that the mass of the lightest right handed fermion is very close to the $Z_2$ odd scalar counterparts in order to enhance the coannihilations. This is evident from the bottom panel plot of Fig. 6 which shows the co-annihilation among the lightest N and $Z_2$ odd scalar components contribute dominantly to lepton asymmetry. Along with that we have an interesting feature of $\lambda _5$. In this we have seen that for particular parameter set as shown in Table 5, the upper bound on $\lambda _5$ is set by the requirement of the required leptonic asymmetry needed to explain the observed baryonic asymmetry in the Universe and the lower bound is set by the Lepton flavor violating process $\mathrm{BR}(\mu \rightarrow e \gamma ) < 4.2 \times 10^{-13}$ [39], as we discuss below. Since the dark matter candidate is a fermion singlet in this case, the parameter $\lambda _5$ is not constrained by dark matter direct detection. While there exist some wash-out effects for the coannihilating processes, there is not much of that for the annihilation processes, as the processes are already out of equilibrium when wash-out becomes effective. In this case, large Yukawa couplings are required to achieve successful leptogenesis, which in turn leads to very small $\lambda _5$, but yet does not affect fermion DM phenomenology much. The Yukawa couplings can be represented by the following matrix for the chosen benchmark point in this scenario.

$$\begin{aligned} y\vert _{NH}&= \begin{pmatrix} 1.01 - 1.14i &{} 2.43 + 1.09i &{} 1.44 + 1.8i \\ -0.15 - 1.21i &{} 1.92 - 9.29i &{} 1.67 + 0.13i \\ -1.40 - 0.61i &{} 2.55 - 2.6i &{} 1.37 - 1.81i \\ \end{pmatrix} \nonumber \\&\quad \times 10^{-2}, \end{aligned}$$

(26)

$$\begin{aligned} y\vert _{IH}&= \begin{pmatrix} 1.88 - 3.05i &{} 1.28 + 0.55i &{} -1.25 - 0.45i \\ -1.12 - 2.92i &{} 1.04 - 0.9i &{} -1.01 + 0.92i \\ -3.87 - 0.79i &{} -0.59 - 1.53i &{} 0.61 + 1.49i \\ \end{pmatrix} \nonumber \\&\quad \times 10^{-2}. \end{aligned}$$

(27)

Fig. 6

Top panel: comoving number densities of $Z_2$ odd particles and lepton asymmetry as a function of temperature for fermion DM scenario. The left (right) panel corresponds to the normal (inverted) hierarchy of neutrino mass spectrum.Bottom panel: rates of individual wash-out processes along with the Hubble rate (left panel); individual (co)annihilation contributions of each channel to the asymmetry. The solid (dashed) black line corresponds to the baryon asymmetry (DM abundance) observed at present epoch ($T \sim 0$). Shaded regions represent the epochs after sphaleron freeze out

Fig. 7

Predictions for LFV processes for $10^2~\mathrm{GeV}< M_{\mathrm{DM}}< 10^3~\mathrm{GeV}$. The two benchmark points are highlighted with red and blue coloured points. The $\mu _\eta$ are varied from 100 GeV to 1 TeV, and the other parameters are taken as presented in Table 4

We then use the SPheno 3.1 [37] interface to check the constraints from flavour data. We particularly focus on three charged lepton flavour violating (LFV) decays namely, $\mu \rightarrow e \gamma , \mu \rightarrow 3e$ and $\mu \rightarrow e$ (Ti) conversion that have strong current limit as well as good future sensitivity [24]. The present bounds are: $\mathrm{BR}(\mu \rightarrow e \gamma ) < 4.2 \times 10^{-13}$ [39], $\mathrm{BR}(\mu \rightarrow 3e) < 1.0 \times 10^{-12}$ [40], $\mathrm{CR} (\mu , \mathrm Ti \rightarrow e, \mathrm Ti) < 4.3 \times 10^{-12}$ [41]. While the future sensitivity of the first two processes are around one order of magnitude lower than the present branching ratios, the $\mu$ to e conversion (Ti) sensitivity is supposed to increase by six order of magnitudes [24] making it a highly promising test to confirm or rule out different TeV scale BSM scenarios. It should be noted that such charged LFV process arises in the SM at one loop level and remains suppressed by the smallness of neutrino masses, much beyond the current and near future experimental sensitivities. Therefore, any experimental observation of such processes is definitely a sign of BSM physics, like the one we are studying here. We show the predictions for LFV processes in our model in Fig. 7, also highlighting the second benchmark point (BP2) mentioned above. The similar contributions for the BP1 scenario remain far more suppressed due to smallness of the corresponding Yukawa couplings. The scatter plot in Fig. 7 is obtained by only varying $\mu _\eta$ from 100 GeV to 1 TeV, and fixing the other parameters with the values presented in Table 4. It can be seen that some part of the parameter space, specially the region which generates correct DM abundance, lies close to the current experimental limits. As mentioned earlier, this bound on $\mathrm{BR}(\mu \rightarrow e \gamma )$ decides the lower bound on the parameter $\lambda _5$ shown in Table 5. If $\lambda _5$ is lower than the one chosen in BP2, the corresponding Yukawa couplings will be bigger (from Casas-Ibarra parametrisation) enhancing the decay rate. For $\mu \rightarrow e \gamma$, the latest MEG 2016 limit [39] can already rule out several points. The promising future sensitivity of the $\mu$ to e conversion (Ti) will be able to explore most part of the parameter space.

4 Conclusion

We have proposed a scenario where baryogenesis via leptogenesis can be achieved through annihilations and coannihilations of particles belonging to a $Z_2$ odd sector including t-channel processes. We have considered a popular model known as the scotogenic model to implement the idea, and addressed the the possibility of explaining the coincidence of DM abundance and baryon asymmetry in the present universe along with non-zero neutrino masses. Pointing out two different possible scenarios corresponding to scalar and fermion DM respectively, we show the non-trivial role played by t-channel annihilation as well as coannihilation processes between different $Z_2$ odd particles. For the two benchmark points chosen in our work, we could obtain successful leptogenesis along with other requirements like DM relic, DM direct detection and light neutrino masses for $M_{\mathrm{DM}} \sim 870$ GeV whereas vanilla leptogenesis in scotogenic model works for $M_1 \ge 10$ TeV. Another interesting feature is the testability of the model at DM direct detection and rare decay experiments. Even though the particle spectrum is in a few ${\mathcal {O}}(100)$ GeV regime or above, away from the reach of current collider experiments, the model can still be tested at near future run of rare decay experiments looking for charged lepton flavour violation like $\mu \rightarrow e \gamma , \mu \rightarrow 3 e$, $\mu$ to e conversion etc. We highlight our interesting results by adopting benchmark points here and leave a detailed numerical analysis of this scenario to an upcoming work.

Appendices

Appendix A: Details of the Boltzmann Equations

In this section, we present the Boltzmann Equation (BE) used in our analysis, in detail. The standard BE is written as follows:

$$\begin{aligned} \frac{d X_\alpha }{dz}&= \frac{\Delta n_\alpha }{sH(z)z} \frac{X_{i_1}X_{i_2}\cdots }{X^{\mathrm{eq}}_{i_1}X^{\mathrm{eq}}_{i_1}\cdots }\gamma ^{\mathrm{eq}} \nonumber \\&\quad \times (i_1 + i_2 \cdots \rightarrow f_1 + f_2 + \cdots ), \end{aligned}$$

(A1)

where $X_{\alpha }=n_{\alpha }/s$ is the comoving number density, $z = M_r/T$ ($M_r$ is the reference scale which is same as the mass of the decaying particle in usual leptogenesis/baryogenesis while being equal to the mass of dark matter in the case of dark matter freeze-out ) and $\Delta n_\alpha$ is the change of number of particle in the process considered. Along with that, the Hubble expansion rate, H(z), and entropy density, s, are defined respectively as

$$\begin{aligned} H(z)&= \sqrt{\frac{4\pi ^3g_*}{45}}\frac{T^2}{M_{\mathrm{PL}}} = \sqrt{\frac{4\pi ^3g_*}{45} }\frac{M_r^2}{M_{\mathrm{PL}}z^2} = \frac{H(M_r)}{z^2}, \end{aligned}$$

(A2)

$$\begin{aligned} s&= g_* \frac{2\pi ^2}{45}T^3. \end{aligned}$$

(A3)

Now, the collision term $\gamma ^{\mathrm{eq}}(i_1 + i_2 \cdots \rightarrow f_1 + f_2 + \cdots )$ are of two types: 1) decay and 2) annihilation processes. For decay process it is given as:

$$\begin{aligned}&\gamma ^{\mathrm{eq}}_D(i \rightarrow f_1 + f_2 + \cdots ) \nonumber \\&\quad = \frac{g_{i} M^2_{i}}{2\pi ^2}TK_1(M_{i}/T)\Gamma (i \rightarrow f_1 + f_2 + \cdots ), \end{aligned}$$

(A4)

and for the scattering processes it is given as

$$\begin{aligned}&\gamma ^{\mathrm{eq}}_{s}(i_1 + i_2 \rightarrow f_1 + f_2 + \cdots ) \nonumber \\&\quad = \frac{g_{i_1}g_{i_2}T}{64\pi ^4}\int ^\infty _{s_{min}} ds \frac{2\lambda (s,m^2_{i_1},m^2_{i_2})}{s}\sigma \sqrt{s}K_1(\sqrt{s}/T), \end{aligned}$$

(A5)

where $K_n$ is the nth order Modified Bessel function of second kind, $\lambda (a,b,c) = a^2 + b^2 + c^2 -2ab - 2ac - 2bc$ and $s_{min}$ = Max$\{(m_{i_1} + m_{i_2})^2,(m_{f_1} + m_{f_2} + \cdots )^2\}$.

Now, we present the explicit expressions of the BEs for $N_k$ and $\eta$ in the case of scalar dark matter.

$$\begin{aligned} \frac{dY_{N_k}}{dz}&= -\frac{1}{szH(z)} \times \left[ \left( \frac{Y_{N_k}}{Y^{\mathrm{eq}}_{N_k}}-1\right) \gamma _D(N_k \rightarrow L_\alpha \eta ) \right. \nonumber \\&\quad + \left( \frac{Y_{N_k}}{Y^{\mathrm{eq}}_{N_k}} \frac{Y_{\eta }}{Y^{\mathrm{eq}}_{\eta }}-1\right) \gamma ^{\mathrm{eq}}_{s}(\eta N_k\rightarrow L\mathrm{SM}) \nonumber \\&\quad + \left. \sum _{l=1}^{3}\left( \frac{Y_{N_k}Y_{N_l}}{Y^{\mathrm{eq}}_{N_k}Y^{\mathrm{eq}}_{N_l}}-1\right) \gamma ^{\mathrm{eq}}_{s}(N_l N_k\rightarrow \mathrm{SM SM})\right] , \nonumber \\ \frac{dY_{\eta }}{dz}&= \frac{1}{szH(z)} \times \left[ \left( \frac{Y_{N_k}}{Y^{\mathrm{eq}}_{N_k}}-1\right) \gamma _D(N_i \rightarrow L_\alpha \eta ) \right. \nonumber \\&\quad - 2\left( \frac{Y^2_{\eta }}{(Y^{\mathrm{eq}}_{\eta })^2}-1\right) \gamma ^{\mathrm{eq}}_{s}(\eta \eta \rightarrow \mathrm{SM SM}) \nonumber \\&\quad \left. - \sum ^3_{m=1}\left( \frac{Y_{N_m}Y_{\eta }}{Y^{\mathrm{eq}}_{N_m}Y^{\mathrm{eq}}_{\eta }}-1 \right) \gamma ^{\mathrm{eq}}_{s}(\eta N_m\rightarrow L\mathrm{SM}) \right] . \end{aligned}$$

(A6)

And similarly expressions of BEs for $N_k$ and $\eta$ in the case of fermion dark matter.

$$\begin{aligned} \frac{dY_{N_k}}{dz}&= -\frac{1}{szH(z)} \nonumber \\&\quad \times \left[ \left( \frac{Y_{N_k}}{Y^{\mathrm{eq}}_{N_k}}-1\right) \gamma _D(N_k \rightarrow L_\alpha \eta ) \right. \nonumber \\&\quad + \left( \frac{Y_{N_k}}{Y^{\mathrm{eq}}_{N_k}}\frac{Y_{\eta }}{Y^{\mathrm{eq}}_{\eta }}-1\right) \gamma ^{\mathrm{eq}}_{s}(\eta N_k\rightarrow L\mathrm{SM}) \nonumber \\&\quad \left. + \sum _{l=1}^{3}\left( \frac{Y_{N_k}Y_{N_l}}{Y^{\mathrm{eq}}_{N_k}Y^{\mathrm{eq}}_{N_l}}-1\right) \gamma ^{\mathrm{eq}}_{s}(N_l N_k\rightarrow \mathrm{SM SM})\right] , \nonumber \\&\quad \quad \text { for} ~k = 2,3 \nonumber \\ \frac{dY_{N_1}}{dz}&= \frac{1}{szH(z)} \nonumber \\&\quad \times \left[ \left( \frac{Y_{\eta }}{Y^{\mathrm{eq}}_{\eta }}-1\right) \gamma _D(\eta \rightarrow L_\alpha N_1) \right. \nonumber \\&\quad - \left( \frac{Y_{N_1}}{Y^{\mathrm{eq}}_{N_1}} \frac{Y_{\eta }}{Y^{\mathrm{eq}}_{\eta }}-1\right) \gamma ^{\mathrm{eq}}_{s}(\eta N_1\rightarrow L\mathrm{SM}) \nonumber \\&\quad \left. - \sum ^3_{l=1}\left( \frac{Y_{N_1}Y_{N_l}}{Y^{\mathrm{eq}}_{N_1}Y^{\mathrm{eq}}_{N_l}}-1\right) \gamma ^{\mathrm{eq}}_{s}(N_l N_1\rightarrow \mathrm{SM SM})\right] , \nonumber \\ \frac{dY_{\eta }}{dz}&= \frac{1}{szH(z)} \nonumber \\&\quad \times \left[ \sum ^3_{k=2}\left( \frac{Y_{N_k}}{Y^{\mathrm{eq}}_{N_k}}-1\right) \gamma _D(N_k \rightarrow L_\alpha \eta ) \right. \nonumber \\&\quad - \left( \frac{Y_{\eta }}{Y^{\mathrm{eq}}_{\eta }}-1\right) \gamma _D(\eta \rightarrow L_\alpha N_1) \nonumber \\&\quad - 2 \left( \frac{Y^2_{\eta }}{(Y^{\mathrm{eq}}_{\eta })^2}-1\right) \gamma ^{\mathrm{eq}}_{s}(\eta \eta \rightarrow \mathrm{SM SM}) \nonumber \\&\quad \left. - \sum ^3_{m=1} \left( \frac{Y_{N_m}}{Y^{\mathrm{eq}}_{N_m}}\frac{Y_{\eta }}{Y^{\mathrm{eq}}_{\eta }}-1\right) \gamma ^{\mathrm{eq}}_{s}(\eta N_m\rightarrow L\mathrm{SM}) \right] . \end{aligned}$$

(A7)

Finally the expressions of BEs for the lepton numbers as follows:

$$\begin{aligned} \frac{dY_L}{dz}&= \frac{1}{sH(z)z} \nonumber \\&\quad \times \left[ \sum _j\left( \left( \frac{Y_{N_j}}{Y^{eq}_{N_j}} - \frac{Y_L}{Y^{\mathrm{eq}}_L}\right) (1+\epsilon _{N_j})\gamma _D(N_j \rightarrow L_\alpha \eta )\right) \right. \nonumber \\&\quad + 2\left( \frac{Y^2_{\eta }}{(Y_{\eta }^{\mathrm{eq}})^2} - \frac{Y^2_L}{(Y_L^{\mathrm{eq}})^2}\right) (1 + \epsilon _{\eta \eta })\gamma ^{\mathrm{eq}}_{s}(\eta \eta \rightarrow L L) \nonumber \\&\quad + \sum _i \left( \frac{Y_{\eta }Y_{N_i}}{Y^{\mathrm{eq}}_{\eta }Y^{\mathrm{eq}}_{N_i}} - \frac{Y_L}{Y^{\mathrm{eq}}_L}\right) \nonumber \\&\quad \times (1 + \epsilon _{N_i \eta } )\gamma ^{\mathrm{eq}}_{s} (\eta N_i\rightarrow L\mathrm{SM}) \nonumber \\&\quad + \sum _i\left( \frac{Y^{eq}_{Ni}}{Y^{eq}_{N_i}} \frac{Y^{eq}_{X}}{Y^{eq}_{X}}-\frac{Y^{eq}_{\eta }}{Y^{eq}_{\eta }} \frac{Y_L}{Y^{eq}_L}\right) \nonumber \\&\quad \times (1+\epsilon _{\eta L})\gamma ^{eq}_s(N_iX\rightarrow \eta L)\nonumber \\&\quad + \sum _i\left( \frac{Y^{eq}_{\eta }}{Y^{eq}_{\eta }} \frac{Y^{eq}_{X}}{Y^{eq}_{X}}-\frac{Y^{eq}_{Ni}}{Y^{eq}_{N_i}} \frac{Y_L}{Y^{eq}_L}\right) \nonumber \\&\quad \times (1+\epsilon _{N_i L})\gamma ^{eq}_s(\eta X\rightarrow N_i L) \nonumber \\&\quad -\frac{Y_{\Delta L}}{Y^{\mathrm{eq}}_L}\gamma ^{\mathrm{eq}}_{s}(\eta L \rightarrow \eta {\overline{L}}) \nonumber \\&\quad \left. + \left( \frac{Y_\eta }{Y^{\mathrm{eq}}_\eta } - \frac{Y_L}{Y^{\mathrm{eq}}_L}\right) \gamma _D(\eta \rightarrow N_1 L)\right] , \end{aligned}$$

(A8)

$$\begin{aligned} \frac{dY_{{\bar{L}}}}{dz}&= \frac{1}{sH(z)z} \nonumber \\&\quad \times \left[ \sum _j\left( \left( \frac{Y_{N_j}}{Y^{\mathrm{eq}}_{N_j}} - \frac{Y_{{\bar{L}}}}{Y^{\mathrm{eq}}_L}\right) (1-\epsilon _{N_j})\gamma _D(N_j \rightarrow L_\alpha \eta )\right) \right. \nonumber \\&\quad + 2\left( \frac{Y^2_{\eta }}{(Y_{\eta }^{\mathrm{eq}})^2} - \frac{Y^2_{{\bar{L}}}}{(Y_L^{\mathrm{eq}})^2}\right) (1 - \epsilon _{\eta \eta })\gamma ^{\mathrm{eq}}_{s}(\eta \eta \rightarrow L L) \nonumber \\&\quad + \sum _i \left( \frac{Y_{\eta }Y_{N_i}}{Y^{\mathrm{eq}}_{\eta }Y^{\mathrm{eq}}_{N_i}} - \frac{Y_{{\bar{L}}}}{Y^{\mathrm{eq}}_L}\right) \nonumber \\&\quad \times (1 - \epsilon _{N_i \eta } )\gamma ^{\mathrm{eq}}_{s}(\eta N_i\rightarrow L\mathrm{SM}) \nonumber \\&\quad + \sum _i\left( \frac{Y^{eq}_{Ni}}{Y^{eq}_{N_i}} \frac{Y^{eq}_{X}}{Y^{eq}_{X}}-\frac{Y^{eq}_{\eta }}{Y^{eq}_{\eta }} \frac{Y_{{\bar{L}}}}{Y^{eq}_L}\right) \nonumber \\&\quad \times (1-\epsilon _{\eta L})\gamma ^{eq}_s(N_iX\rightarrow \eta L)\nonumber \\&\quad + \sum _i\left( \frac{Y^{eq}_{\eta }}{Y^{eq}_{\eta }} \frac{Y^{eq}_{X}}{Y^{eq}_{X}}-\frac{Y^{eq}_{Ni}}{Y^{eq}_{N_i}} \frac{Y_{{\bar{L}}}}{Y^{eq}_L}\right) \nonumber \\&\quad \times (1-\epsilon _{N_i L}) \gamma ^{eq}_s(\eta X\rightarrow N_i L) \nonumber \\&\quad +\frac{Y_{\Delta L}}{Y^{\mathrm{eq}}_L}\gamma ^{\mathrm{eq}}_{s}(\eta L \rightarrow \eta {\overline{L}}) \nonumber \\&\quad \left. + \left( \frac{Y_\eta }{Y^{\mathrm{eq}}_\eta } - \frac{Y_{{\bar{L}}}}{Y^{\mathrm{eq}}_L}\right) \gamma _D(\eta \rightarrow N_1 L)\right] . \end{aligned}$$

(A9)

where the leptonic comoving number density are defined as $Y_L = Y^{\mathrm{eq}}_L+ Y_{\Delta L}/2, Y_{{\overline{L}}} = Y^{\mathrm{eq}}_L - Y_{\Delta L}/2$. The rates $\gamma _D(N_j\rightarrow L_\alpha \eta )$ are the decay process considered in [21] and $\gamma ^{\mathrm{eq}}_s(\eta \eta \rightarrow LL)$ are shown in Fig. 5, whereas $\gamma ^{\mathrm{eq}}_{s}(\eta N_i\rightarrow L\mathrm{SM})$ and $\gamma ^{\mathrm{eq}}_s(N_iX\rightarrow \eta L)$ are shown in Fig. 1. The processes $\gamma ^{\mathrm{eq}}_s(N_iX\rightarrow \eta L)$ and $\gamma ^{\mathrm{eq}}_s(\eta X\rightarrow N_i L)$ are the same as Fig. 1 by just interchanging the one of the initial with the final particle. And finally the processes $\gamma ^{\mathrm{eq}}_{s}(\eta L \rightarrow \eta {\overline{L}})$ and $\gamma _D(\eta \rightarrow N_1 L)$ are shown in Fig. 8. There are no other further processes which contributes in the above Boltzmann equation. Now taking the difference between Eq. (A8) and Eq. (A9) and keeping the terms asymmetry to leading order i.e $Y^2_L \simeq (Y^{\mathrm{eq}}_L)^2 + Y^{\mathrm{eq}}_LY_{\Delta L}$, one would get the final Boltzmann equation for asymmetry Eq. 12)

$$\begin{aligned} \frac{dY_{\Delta L}}{dz}&= \frac{1}{sH(z)z} \nonumber \\&\quad \times \left[ \sum _j\left( \left( \frac{Y_{N_j}}{Y^{eq}_{N_j}} - 1\right) \epsilon _{N_j}\gamma _D(N_j \rightarrow L_\alpha \eta ) \right. \right. \nonumber \\&\quad \left. - \frac{Y_{\Delta L}}{Y^{\mathrm{eq}}_L}\gamma _D(N_j \rightarrow L_\alpha \eta ) \right) \nonumber \\&\quad + 2\left( \frac{Y^2_{\eta }}{(Y_{\eta }^{\mathrm{eq}})^2} - 1\right) \epsilon _{\eta \eta }\gamma ^{\mathrm{eq}}_{s}(\eta \eta \rightarrow L L) \nonumber \\&\quad - 2\frac{Y_{\Delta L}}{Y^{\mathrm{eq}}_L}\gamma ^{\mathrm{eq}}_{s}(\eta \eta \rightarrow L L) \nonumber \\&\quad + \sum _i \left( \left( \frac{Y_{\eta }Y_{N_i}}{Y^{\mathrm{eq}}_{\eta }Y^{\mathrm{eq}}_{N_i}} - 1\right) \epsilon _{N_i \eta }\gamma ^{\mathrm{eq}}_{s}(\eta N_i\rightarrow L\mathrm{SM}) \right. \nonumber \\&\quad \left. - \frac{Y_{\Delta L}}{Y^{\mathrm{eq}}_L}\gamma ^{\mathrm{eq}}_{s}(\eta N_i\rightarrow L\mathrm{SM})\right) \nonumber \\&\quad - \sum _i\frac{Y_{\Delta L}}{Y^{\mathrm{eq}}_L}\left( \gamma ^{eq}_s(N_iX\rightarrow \eta L) + \gamma ^{eq}_s(\eta X\rightarrow N_i L) \right) \nonumber \\&\quad \left. -2\frac{Y_{\Delta L}}{Y^{\mathrm{eq}}_L}\gamma ^{\mathrm{eq}}_{s}(\eta L \rightarrow \eta {\overline{L}}) - \frac{Y_{\Delta L}}{Y^{\mathrm{eq}}_L}\gamma _D(\eta \rightarrow N_1 L)\right] , \end{aligned}$$

(A10)

Now, one would notice that in Eqs. (A8),(A9) and (A10) the index j runs from $1-3$ if scalar is the Dark Matter in which case the last decay term $\gamma _D(\eta \rightarrow N_1 L)=0$. But, for the case of $N_1$ as the dark matter j runs from $2-3$ and the last decay term $\gamma _D(\eta \rightarrow N_1 L)\ne 0$. Now, one may notice that from CPT invariance and unitarity terms proportional to $\epsilon _{N_i L}$ and $\epsilon _{\eta L}$ cancel out exactly. Hence, the effects coming from the $\gamma ^{\mathrm{eq}}_s(N_iX\rightarrow \eta L)$ and $\gamma ^{\mathrm{eq}}_s(\eta X\rightarrow N_i L)$ only contribute to the wash-out but it is suppressed compared to the inverse decay processes as shown in Fig. 3 and Fig. 6.

Fig. 8

Diagrams leading to $\gamma ^{\mathrm{eq}}_{s}(\eta L \rightarrow \eta {\overline{L}})$(left) and $\gamma _D(\eta \rightarrow N_1 L)$(right), respectively

Appendix B: Details of the asymmetry

In this section we give the details of the asymmetry shown in Eq. (B9). We would first start with the the basic general expression giving the asymmetry shown as follows:

$$\begin{aligned} \delta = 4\mathfrak {I}[{\mathcal {C}}^{*}_0{\mathcal {C}}_1]\mathfrak {I}[{\mathcal {A}}_0{\mathcal {A}}_1]|{\mathcal {W}}|^2, \end{aligned}$$

(B1)

where ${\mathcal {W}}$ corresponds to the wavefunction of the incoming and outgoing particles, ${\mathcal {C}}$’s corresponds to the couplings of the tree (${\mathcal {C}}_0$) and loop (${\mathcal {C}}_1$) and ${\mathcal {A}}$’s corresponds the rest of the amplitude respectively.

Now, starting with the tree level amplitudes

$$\begin{aligned} {\mathcal {M}}_{0}&= i\frac{x^\dagger _\alpha {\bar{\sigma }}(p_1+p_2)\sigma y^{\dagger }_{N_i}}{s} + i\frac{x^\dagger _\alpha y^{\dagger }_{N_i}(2p_1-p_3)^\mu }{t - m^2_\eta }, \end{aligned}$$

(B2)

where the amplitude ${\mathcal {M}}_{0}$ corresponds to the the tree diagram shown in Fig. 1. The corresponding amplitudes for the loop correction are given as follows

$$\begin{aligned} {\mathcal {M}}_{1}&= \frac{x^\dagger _\alpha {\bar{\sigma }}^\mu (p_1 + p_2)\sigma p_1{\bar{\sigma }}x_{N_i}m_{N_j}}{s}{\mathcal {C}}_1 \nonumber \\&\quad + \frac{x^\dagger _\alpha {\bar{\sigma }}^\mu (p_1 + p_2)\sigma p_2{\bar{\sigma }}x_{N_i}m_{N_j}}{s}{\mathcal {C}}_2 \nonumber \\&\quad + \frac{x^\dagger _\alpha y^\dagger _{N_i}}{t-m^2_\eta } (2p_1-p3)^\mu {\mathcal {C}}_3 \nonumber \\&\quad + \left( \frac{x^\dagger _\alpha {\bar{\sigma }}^\mu (p_1+p_2)\sigma y^\dagger _{N_i}}{s} + \frac{x^\dagger _\alpha y^\dagger _{N_i}(2p_1-p3)^\mu }{(t-m^2_\eta )}\right) . \end{aligned}$$

(B3)

The cross term coming from the above two expressions will give us

$$\begin{aligned}&\mathfrak {I}[{\mathcal {A}}_0{\mathcal {A}}_1]|{\mathcal {W}}|^2 \nonumber \\&\quad = \left[ m_{N_i}m_{N_j}\frac{2m^4_{N_j} + 2m^2_\eta t + s(s+t) - m^2_{N_i}(2m^2_\eta + 3s + 2t)}{s(t-m^2_\eta )} \right. \nonumber \\&\quad + \left. 2m_{N_i}m_{N_j}\frac{(t-m^2_\eta )}{s} \right] \mathfrak {I}[{\mathcal {C}}_1] \nonumber \\&\qquad + \left[ 2m_{N_i}m_{N_j}\frac{m^2_\eta -s -t}{s} \right. \nonumber \\&\qquad +\left. m_{N_i}m_{N_j}\frac{2m^4_{N_i} + t(2m^2_\eta - s) - m^2_{N_i}(2m^2_\eta + s + 2t)}{s(m^2_\eta - t)}\right] \nonumber \\&\qquad \times \mathfrak {I}[{\mathcal {C}}_2] \nonumber \\&\qquad + \left[ 2m_{N_i}\frac{(s+t-m^2_{N_i})(m^2_\eta - t)}{(t-m^2_\eta )^2} \right. \nonumber \\&\quad \quad + \left. m_{N_i}\frac{(m^4_{N_i} + t(2m^2_\eta - s) - m^2_{N_i}(2m^2_\eta +s+2t))}{s(m^2_\eta - t)}\right] \nonumber \\&\qquad \times \mathfrak {I}[{\mathcal {C}}_3] \nonumber \\&\qquad + \left[ \frac{2(m^2_{N_i}-t)(m^2_\eta +t)}{(m^2_\eta -t)^2} -\frac{2(s+t-m^2_\eta )}{s} \right. \nonumber \\&\qquad +\left. \frac{2m^2_{N_i}+m_{N_i}(2m^2_\eta -s)t-m^3_{N_i} (2m^2_\eta + s + 2t)}{s(m^2_\eta - t)}\right] \nonumber \\&\qquad \times \frac{m_{N_i}m_{N_j}\Gamma _j}{16\pi ^2((m^2_{N_i}-m^2_{N_j})^2 + m^2_j\Gamma ^2_j)}, \end{aligned}$$

(B4)

where

$$\begin{aligned} {\mathcal {C}}_1&= \frac{(s-m^2_\eta -m^2_{N_i})}{{\mathcal {D}}} \nonumber \\&\quad \times \left[ {\mathcal {B}}_0(m^2_\eta ,0,m^2_{N_j}) - {\mathcal {B}}_0(s,m^2_{N_j},m^2_\eta ) \right. \nonumber \\&\quad + (m^2_\eta - m^2_{N_i})C_0 \nonumber \\&\quad \left. (m^2_\eta ,m^2_{N_i}, (m^2_\eta +m^2_{N_i}-s)/2,0,m^2_\eta ,m^2_{N_j})\right] , \nonumber \\&\quad - \frac{2m^2_{N_i}}{{\mathcal {D}}} \left[ {\mathcal {B}}_0(s,m^2_\eta ,m^2_{N_i}) - {\mathcal {B}}_0(m^2_{N_i},0,m^2_\eta ) \right. \nonumber \\&\qquad + (m^2_{N_i} - m^2_{N_j})C_0 \nonumber \\&\quad \times \left. (m^2_\eta ,m^2_{N_i}, (m^2_\eta +m^2_{N_i}-s)/2,0,m^2_\eta ,m^2_{N_j})\right] \nonumber \\ {\mathcal {C}}_2&= \frac{(s-m^2_\eta -m^2_{N_i})}{{\mathcal {D}}} \nonumber \\&\quad \times \left[ {\mathcal {B}}_0(s,m^2_\eta ,m^2_{N_j}) - {\mathcal {B}}_0(m^2_{N_i},0,m^2_\eta )\right. \nonumber \\&\quad + (m^2_{N_i} - m^2_{N_j})C_0 \nonumber \\&\quad \times \left. (m^2_\eta ,m^2_{N_i}, (m^2_\eta +m^2_{N_i}-s)/2,0,m^2_\eta ,m^2_{N_j})\right] , \nonumber \\&\quad - \frac{2m^2_\eta }{{\mathcal {D}}} \left[ {\mathcal {B}}_0 (m^2_\eta ,0,m^2_{N_j}) - {\mathcal {B}}_0(s,m^2_\eta ,m^2_{N_j}) \right. \nonumber \\&\quad + (m^2_\eta - m^2_{N_i})C_0 \nonumber \\&\quad \times \left. (m^2_\eta ,m^2_{N_i}, (m^2_\eta +m^2_{N_i}-s)/2,0, m^2_\eta ,m^2_{N_j})\right] , \nonumber \\ {\mathcal {C}}_3&= \frac{1}{m^2_{N_i}-t} \nonumber \\&\quad \times \left[ {\mathcal {B}}_0 (t,0,m^2_{N_j}) - {\mathcal {B}}_0(m^2_{N_j},m^2_\eta ,0) \right. \nonumber \\&\quad + (m^2_\eta - m^2_{N_j} + t - m^2_{N_i}){\mathcal {C}}_0 \nonumber \\&\quad \times \left. (m^2_\eta ,m^2_{N_i},(m^2_{N_i}-t)/2,m^2_\eta ,m^2_{N_j},0)\right] , \nonumber \\ {\mathcal {D}}&= (s-m^2_\eta - m^2_{N_i})^2 - 4m^2_\eta m^2_{N_i}. \end{aligned}$$

(B5)

where the scalar integral ${\mathcal {B}}_0$ and ${\mathcal {C}}_0$ is given by

$$\begin{aligned}&{\mathcal {B}}_0(p^2,m^2_1,m^2_3) = \mu ^{2\epsilon }\int \frac{d^d}{(2\pi )^d} \frac{1}{(l^2-m^2_1)((l+p)^2-m^2_2)}, \nonumber \\&{\mathcal {C}}_0(p^2_1,p^2_2,p1.p2,m^2_1,m^2_2,m^2_3) \nonumber \\&\quad = \mu ^{2\epsilon }\int \frac{d^d}{(2\pi )^d} \nonumber \\&\quad \times \frac{1}{(l^2-m^2_1)((l+p_1)^2-m^2_2)((l+p_2)^2-m^2_3)}. \end{aligned}$$

(B6)

Now the thermally averaged cross section is given as

$$\begin{aligned} \langle \sigma _\delta v\rangle&= \frac{T}{32\pi m^2_\eta m^2_{N_i} K_2(m_\eta /T)K_2(m_{N_i}/T)} \nonumber \\&\quad \times \int ^\infty _{(m_\eta + m_{N_i})^2}\int ^{t_+}_{t_-}ds dt |\delta |^2\frac{p_{in}p_{out}K_1(\sqrt{s}/T)}{\sqrt{s}}, \nonumber \\ p_{in}&= \frac{1}{2}\frac{\lambda ^{1/2} (s,m^2_\eta ,m^2_{N_i})}{\sqrt{s}}; \nonumber \\ p_{out}&= \frac{\sqrt{s}}{2}; \quad t_{\pm } = \frac{1}{2}(m^2_\eta + m^2_{N_i} - s \pm p_{in}p_{out}). \end{aligned}$$

(B7)

Now taking the s-wave approximation of the above expression

$$\begin{aligned} \langle \sigma _\delta v \rangle&\simeq \sum _{ij}\frac{\mathfrak {I}[(yy^\dagger )^2_{ij}]}{8\pi (1+x_{N_i})} {\widetilde{\epsilon }}_{ij}(x^2_i(x_i-3)), \end{aligned}$$

(B8)

$$\begin{aligned} {\widetilde{\epsilon }}_{ij}&= \frac{\sqrt{x_j}}{6 x_i\left( -x_i^{3/2}+x_i (x_j-2)+\sqrt{x_i} x_j+1\right) ^2(\sqrt{x_i}-3)} \nonumber \\&\quad \times \left( x_i^{7/2} (3 x_j+1) +\sqrt{x_i} (3 x_j+5)+1 \right. \nonumber \\&\quad - 3 x_i^{5/2} \left( x_j \left( D+(x_j-3) x_j+4\right) - 3 D -2\right) \nonumber \\&\quad -3 x_i^{3/2} \left( 2 \left( D+3\right) + x_j \left( x_j \left( D+x_j+1\right) -D-4\right) \right) \nonumber \\&\quad - x_i^4+f^3 \left( 3 D+3 x_j^2+11\right) \nonumber \\&\quad - 3 x_i^2 \left( x_j \left( D+2 (x_j-1) x_j+2\right) - D +6\right) \nonumber \\&\quad \left. +x_i \left( 1-3 x_j \left( D+x_j-4\right) \right) \right) \nonumber \\&\quad + \frac{\sqrt{x_j}}{4 x_i} \left( \sqrt{x_i}-1+\frac{\sqrt{x_j}}{(1+\sqrt{x_i})^2} (\sqrt{x_i}-1+r_j) \right. \nonumber \\&\quad \times \left( \log \left( \frac{1+\sqrt{x_i}x_j}{x_i(1+\sqrt{x_i})}\right) \right. \nonumber \\&\quad -\log \left( \frac{1+x_i+x^{3/2}_i+ \sqrt{x_i}x_j}{x_i(1+\sqrt{x_i})}\right) \nonumber \\&\quad + \left. \left. \log \left( 1+\frac{1+\sqrt{x_i}}{\sqrt{x_i} (\sqrt{x_i}-1+x_i+x_j)}\right) \right) \right) \nonumber \\&\quad + \frac{\sqrt{x_i} x_j{\tilde{\Gamma }}_j}{\pi ((x_i-x_j)^2 + x_j{\tilde{\Gamma }}^2_j)}, \nonumber \\ D&= \sqrt{(x_i-x_j) \left( x_i+4 \sqrt{x_i}-x_j+4\right) }, \nonumber \\ x_l&= \frac{M^2_{N_l}}{m^2_\eta }, \qquad {\tilde{\Gamma }}_j = \frac{\Gamma _j}{m_\eta }. \end{aligned}$$

(B9)