---
author:
- 'Shao-Ping Li,'
- 'Xin-Qiang Li,'
- 'Xin-Shuai Yan'
- 'and Ya-Dong Yang'
bibliography:
- 'reference.bib'
title: 'Freeze-in Dirac neutrinogenesis: thermal leptonic CP asymmetry'
---

Introduction {#sec:intro}
============

Recent developments in particle physics and cosmology, especially those related to the neutrino mass, dark matter, as well as baryon asymmetry of the Universe (BAU), have highlighted the importance of feeble couplings. Actually, feeble couplings are already present in the Yukawa couplings of the light charged fermions within the Standard Model (SM); *e.g.*, the SM predicts an electron Yukawa coupling with $y_e\simeq 10^{-6}$. If one also accepts feeble Yukawa couplings of the Dirac neutrinos, the smallness of neutrino masses can then be simply addressed via the Higgs-like mechanism with three right-handed Dirac neutrino singlets. Feeble couplings can also play an important role in the early Universe. For example, the feebleness allows a freeze-in production for the dark matter abundance, which can be effectively kept from large annihilation [@Hall:2009bx; @Bernal:2017kxu]. Moreover, the feebleness stirs up a new leptogenesis, named Dirac neutrinogenesis (DN) [@Dick:1999je], in which the out-of-equilibrium condition for generating the lepton-number ($L$) asymmetry can be guaranteed and the baryon-number ($B$) asymmetry is generated via thermal sphaleron transitions [@Kuzmin:1985mm], even in a theory with $B-L=0$ initially.

In the typical versions of DN mechanism [@Dick:1999je; @Murayama:2002je; @Cerdeno:2006ha; @Gu:2007mc; @Bechinger:2009qk; @Narendra:2017uxl], the lepton-number asymmetry is generated by heavy particle decays with a non-thermal distribution. In addition, to discuss the loop correction for nonzero kinetic phase (or the absorptive part of the decay amplitude), one usually focuses on the new particle sector, which is also the case in seesaw-based leptogenesis [@Giudice:2003jh; @Buchmuller:2004nz; @Davidson:2008bu], while the contribution from wavefunction correction in the lepton-doublet sector has not yet been considered to the best of our knowledge.

There could be two possible reasons for having neglected the lepton-doublet wavefunction contribution. On the one hand, the charged-lepton flavors are widely assumed to populate in the diagonal basis, and thus the leptonic CP asymmetry cannot be generated from self-energy diagrams in the lepton-doublet sector. Interestingly, however, it has been pointed out earlier that the well-known tri-bimaximal (TB) mixing pattern [@Harrison:2002er], with a minimal correction from the charged-lepton or neutrino sector, can produce compatible neutrino oscillaton data while retaining its compelling prediction [@Albright:2008rp; @He:2011gb]. In this respect, a nontrivial combination of the charged-lepton and neutrino mixings is preferred to produce the oscillation observables. On the other hand, even with a non-diagonal charged-lepton Yukawa matrix, there is no on-shell cut in the self-energy loop at zero-temperature regime, and hence no CP asymmetry either. Nevertheless, it has been illustrated in ref. [@Giudice:2003jh] and later implemented in ref. [@Hambye:2016sby] that, at high-temperature regime where thermal effects come into play, the zero-temperature cutting rules should be superseded by the thermal cuts [@Das1997], allowing consequently nonzero contributions to the leptonic CP asymmetry that would otherwise vanish at vacuum regime.

Therefore, as will be exploited in this paper, when both thermal effects and nontrivial mixings in the charged-lepton and neutrino sectors are taken into account, one can expect the leptonic CP asymmetry at finite temperature to carry a nonzero imaginary piece, *i.e.*, $\text{Im}[(Y_\nu Y^\dagger_\nu)(Y_\ell Y^\dagger_\ell)]\neq0$, where $Y_{\ell}$ and $Y_{\nu}$ denote respectively the charged-lepton and neutrino Yukawa matrices that are responsible for their respective masses and mixings. This enables us to exploit a direct interplay between the BAU and the neutrino oscillation observables in a minimal setup, without tuning additional Yukawa couplings beyond $Y_{\ell, \nu}$. Furthermore, since the feeble neutrino Yukawa couplings essentially prompt an out-of-equilibrium condition (*i.e.*, the right-handed Dirac neutrinos undergo a freeze-in production in the early Universe), there is no need to invoke much heavier dynamical degrees of freedom (d.o.f), and the evolution of the lepton-number asymmetry can be much simplified as well.

The remainder of this paper is constructed as follows. We begin in section \[sec:2\] with a brief overview of the DN mechanism, and then calculate the leptonic CP asymmetry with two different thermal cuts in a model-independent way. The Boltzmann equation for the evolution of the lepton-number asymmetry in the freeze-in regime is also derived here. In section \[sec:3\], we discuss the nontrivial mixings in the charged-lepton and neutrino sectors by focusing on minimal corrections to the TB mixing pattern. In section \[sec:4\], we identify the scalars participating in the out-of-equilibrium decay and perform our detailed numerical analyses. Our conclusions are finally made in section \[sec:con\].

Thermal leptonic CP asymmetry and evolution {#sec:2}
===========================================

Dirac neutrinogenesis
---------------------

The basic idea of DN [@Dick:1999je] can be summarized as follows. In a theory without lepton-number-violating Lagrangian, due to the feeble neutrino Yukawa couplings that prevent the left- and right-handed Dirac neutrinos from equilibration (dubbed as “L-R equilibration" from now on), the leptonic CP asymmetry from a heavy particle decay in the early Universe can result in a net lepton-number asymmetry stored in the right-handed Dirac neutrinos and lepton doublets. As the sphaleron transitions act only on the left-handed particles, the net lepton-number asymmetry stored in the lepton doublets will be partially converted to the baryon-number asymmetry via rapid sphaleron processes, while the portion stored in the right-handed Dirac neutrinos keeps intact. After the sphaleron freezes out around $T\simeq \mathcal{O}(100)$ GeV, a net baryon-number (as well as lepton-number) asymmetry survives till today. At the thermal sphaleron epoch, $10^2~\text{GeV}<T<10^{12}~\text{GeV}$, all the SM species (except for the right-handed Dirac neutrinos) are kept in chemical equilibrium, and the conversion fraction between lepton- and baryon-number asymmetries is given by [@Harvey:1990qw] $$\begin{aligned}
\label{sphaleron relation}
Y_{\Delta B}=c\,Y_{\Delta (B-L)}=-c\,Y_{\Delta L},\end{aligned}$$ where $c=(8N_f+4N_H)/(22N_f+13N_H)$, with $N_{f}$ and $N_{H}$ denoting the numbers of fermion generations and Higgs doublets, respectively.

To prompt the necessary out-of-equilibrium condition so as to keep the generated lepton-number asymmetry from being washed out by the L-R equilibration, the lepton-number-violating thermal decay rate must be sufficiently smaller than the expansion rate of the Universe, which typically requires the Dirac neutrino Yukawa couplings to be $Y_\nu\lesssim \mathcal{O}(10^{-8})$. Such feeble couplings, despite of their non-aesthetic nature, are generically present, if the sub-eV Dirac neutrino masses are generated by the Higgs-like mechanism with a vacuum expectation value (VEV) around the electroweak scale. On the other hand, such a mechanism of neutrino mass generation is often criticized on account of *naturalness*, and dynamical explanations of the smallness of neutrino masses are, therefore, more biased by enlarging the Yukawa space and/or introducing sufficiently heavy particles, which have also been considered in explicit realizations of the DN [@Murayama:2002je; @Cerdeno:2006ha; @Gu:2007mc; @Bechinger:2009qk; @Narendra:2017uxl].

Nonetheless, for these dynamical explanations with overabundant Yukawa parameters, reliable phenomenological predictions rely on particular bases and values of the unknown Yukawa couplings beyond those that can be directly fixed by the lepton flavor spectrum. In particular, a simple connection between the BAU and the low-energy neutrino oscillation observables cannot be established, if the DN realization has nothing to do with the Yukawa couplings that are directly responsible for the lepton masses and mixings. Furthermore, it is difficult, if not impossible, to detect the additional sufficiently heavy particles at colliders.

In this paper, as an underlying theory for explaining the feebleness of the Yukawa couplings for both light charged fermions and neutrinos is still unknown, if it were to exist, we shall take the feeble neutrino couplings as a starting point. In this context, we provide a new DN realization in which the decaying particle that generates the lepton-number asymmetry is in thermal equilibrium. For this purpose, we consider a Higgs-like doublet which has a vacuum mass around $\mathcal{O}(10^2)$ GeV and feeble couplings to the right-handed Dirac neutrinos. The feeble Yukawa couplings ensure that the right-handed Dirac neutrinos never reach equilibrium with the thermal bath. On top of that, the leptonic CP asymmetry is induced by the self-energy correction in the lepton-doublet sector due to thermal effects. This realization allows us to establish a simple connection between the BAU and the low-energy neutrino oscillation observables, and, at the same time, renders the detection of the scalars *at least* in principle possible at colliders.

Theoretical setup
-----------------

In this subsection, we shall adopt a real-time formalism in thermal field theory to calculate the thermal leptonic CP asymmetry. To appreciate the subtlety in calculating the CP asymmetry between thermal field theory and non-equilibrium quantum field theory (QFT), we shall use two different thermal cuts and compare the corresponding consequences that arise from the different dependence on the distribution functions. The evolution of the lepton-number asymmetry will be determined by a simplified Boltzmann equation in the freeze-in regime.

### Thermal field theory: real-time formalism

There are two equivalent approaches in thermal field theory, the real-time and the imaginary-time formalism [@Landsman:1986uw; @Das1997]. Within the real-time formalism, we do not need to perform analytic continuation for the physical region, but there is a doubling of d.o.f dual to each field presented in vacuum QFT. As a result, the interaction vertices are doubled, and the thermal propagators have a $2\times 2$ structure. In the following, we shall adopt this formalism to calculate the thermal leptonic CP asymmetry.

![\[thermal FeynRule\] Circling rules in doubled interaction vertices specified by different thermal indices $\pm$. Here $\mathcal{L}_{Y}$ can be either the Yukawa Lagrangian of the SM extended by a neutrino term or of the neutrinophilic two-Higgs-doublet model (2HDM), to be discussed later.](figures/vertices.pdf){width="82.00000%"}

![\[thermal FeynRule2\] Circling rules in thermal propagators. Here $\alpha$ and $\beta$ take the thermal indices $\pm$. Note that the propagator indices in $(c)$ and $(d)$ are completely determined by the uncircled ones, with $\dot{\alpha}$ and $\dot{\beta}$ taking the opposite signs of $\alpha$ and $\beta$, respectively. ](figures/GreenFun.pdf)

Within the real-time formalism, while both the closed-time path formulation and the thermo-field dynamics can be used, we shall follow here the former [@Das1997]. In this formulation, the circling rules necessary for evaluating the absorptive part of the decay amplitude are given in figures \[thermal FeynRule\] (for interaction vertices) and \[thermal FeynRule2\] (for thermal propagators). In order to get a compact expression for the thermal propagators and a unified rule in writing the amplitude for each vertex, we adopt a convention in which the numerator factor $\slashed p \pm m$ of the fermion propagator is decomposed into a spin summation $\sum_s u^s \bar u^s (v^s \bar v^s)$, where the Dirac spinors would then be attached to each vertex. Thus, the thermal propagators, with the subscript indices $\pm$ specifying the corresponding matrix elements, can be written, explicitly, as $$\begin{aligned}
G_{++}(p)&=\frac{i}{p^2-m^2+i \epsilon}\pm 2\pi f_{B/F}(\vert p^0 \vert)\delta(p^2-m^2),\\[0.2cm]
G_{--}(p)&=\left(G_{++}(p)\right)^*,\\[0.2cm]
G_{+-}(p)&=2\pi \left[\pm f_{B/F}(\vert p^0\vert)+\theta(-p^0)\right]\delta(p^2-m^2),\\[0.2cm]
G_{-+}(p)&=2\pi \left[\pm f_{B/F}(\vert p^0\vert)+\theta(p^0)\right]\delta(p^2-m^2),\end{aligned}$$ where $f_{B/F}(E)=(e^{E/T}\mp 1)^{-1}$ are the standard distribution functions, with $B$ and $F$ referring to the bosons and fermions, respectively. $\theta(p^0)$ denotes the Heaviside step function. Note that the circled indices in figures \[thermal FeynRule2\]$(c)$ and \[thermal FeynRule2\]$(d)$ are completely determined by the uncircled ones, with $\dot{\alpha}$ and $\dot{\beta}$ taking the opposite signs of $\alpha$ and $\beta$, respectively. For example, the propagator in figure \[thermal FeynRule2\]$(c)$ with an uncircled thermal index $\alpha=+$ is given by $G_{+-}(p)$.

### Leptonic CP asymmetry: model-independent approach

As a generic model-independent discussion, let us consider the neutrino Yukawa term, $$\begin{aligned}
\label{Yukawa}
-\mathcal{L}_\nu=Y_\nu \bar{L} \tilde{\Phi}\nu_R+ {\rm h.c.},\end{aligned}$$ added to the SM Lagrangian. Here we denote the lepton doublet by $L$, and assume that the Higgs doublet $\tilde{\Phi}\equiv i \sigma_2 \Phi^\ast$, with $\sigma_2$ being the Pauli matrix, does not populate well above the electroweak scale. It could be the SM Higgs doublet or a second Higgs doublet which may or may not couple to quarks. To forbid the appearance of Majorana neutrino mass term and, at the same time, to realize the DN, the right-handed neutrinos must carry a non-zero lepton number under some global $U(1)$ symmetry. After the Higgs doublet develops a non-vanishing VEV, $\langle\Phi\rangle=\left(0,v_\Phi/\sqrt{2}\right)^{T}$, the vacuum neutrino mass is then given by $m_\nu=v_\Phi Y_\nu/\sqrt{2}$.

![\[CPasym\] Leptonic CP asymmetry generated in $\Phi\to L\bar\nu$ decay at $\mathcal{O}(Y_\nu^2 Y_\ell^2)$, where the left and the right diagram represent the tree-level and the one-loop contribution, respectively.](figures/CPasym.pdf){width="86.00000%"}

Now, let us consider the leptonic CP asymmetry generated in $\Phi \to L \bar\nu$ decay. Since $Y_\nu \ll Y_\ell$ is a generic condition for realizing the DN, we shall not consider the CP asymmetry at $\mathcal{O}(Y_\nu^4)$, which is the case in seesaw-based leptogenesis [@Buchmuller:2004nz]. Instead, we shall determine the CP asymmetry at $\mathcal{O}(Y_\nu^2 Y_\ell^2)$. At this order, the absorptive part of the decay amplitude could arise from self-energy diagrams in the lepton-doublet sector, as well as from vertex diagrams if $\Phi$ also couples to the right-handed charged leptons. Here, let us concentrate on the former. Note that the contribution from vertex diagrams is found to be of similar size as that from the self-energy diagrams in the SM Higgs case, and is even absent in the neutrinophilic 2HDM, as will be detailed in section \[sec:4\]. Then, the CP asymmetry may arise from the interference between tree and one-loop diagrams shown in figure \[CPasym\]. At zero-temperature regime, $T=0$, there is no on-shell cut for an electroweak scalar running in the loop. At high-temperature regime, however, due to thermal bath corrections, the propagators can be on shell, producing therefore a nonzero absorptive part in the amplitude [@Giudice:2003jh; @Hambye:2016sby].

The amplitude for $\Phi \to L \bar\nu$ decay can be defined as $i \mathcal{M}\equiv c_0 I_0+c_1 I_1$, where the coupling constants have been factored out into $c_{0,1}$, while all the other factors are contained in $I_{0,1}$, with the subscripts $0$ and $1$ referring respectively to the contributions from tree and one-loop diagrams shown in figure \[CPasym\]. The thermal leptonic CP asymmetry is then given by $$\begin{aligned}
\label{CP definition}
\epsilon_D\equiv \frac{\Gamma(\Phi\to L \bar\nu)-\Gamma(\bar\Phi\to \bar L \nu)}{\Gamma(\Phi\to L \bar\nu)+\Gamma(\bar\Phi\to \bar L \nu)}
\simeq -2\frac{\text{Im}(c_0^* c_1)}{\vert c_0\vert^2} \frac{\text{Im}(I_0^*I_1)}{\vert I_0\vert^2},\end{aligned}$$ where the second equation is obtained in the rest frame of $\Phi$. As all the charged-lepton (neutrino) flavors are in (out of) L-R equilibration before the sphaleron transition decouples, an implicit summation over all final lepton flavors is assumed in the decay. With the flavor indices being specified in figure \[CPasym\], we have $c_0=Y_{\nu,ij}$ and $c_1=Y_{\nu, kj} Y^*_{\ell,kl} Y_{\ell, il}$. It can be seen that diagonal $Y_\nu$ or $Y_\ell$ would lead to $\text{Im}(c_0^* c_1)=0$. Taking into account the thermal masses and neglecting the small neutrino masses, we obtain the tree-level amplitude squared as $$\begin{aligned}
\vert I_0\vert^2=M_\Phi^2(T)-m_{L_i}^2(T).\end{aligned}$$ Here we should mention that thermal corrections to fermions would modify the Dirac equation for a spinor $\psi$ to $[(1+a)\slashed p+b\slashed u] \psi=0$ [@Weldon:1982bn], where $u$ is the four-velocity of the thermal bath, and $a,b$ are temperature-dependent functions (see *e.g.*, ref. [@Giudice:2003jh]). The $a,b$ functions would also modify the fermion propagators and hence the dispersion relations, making the expressions for spin summation and propagator poles quite lengthy and involved. Nevertheless, as illustrated in ref. [@Giudice:2003jh], to a good approximation, both the dispersion relation and the modified Dirac equations can be simplified by replacing the vacuum mass with the thermal one, *i.e.*, $p^2\simeq m^2(T)$. We shall confine ourselves to adopt this approximation in the subsequent calculations.

### Thermal effects: time-ordered and retarded/advanced cuts {#sec:TO&RAcuts}

To calculate the leptonic CP asymmetry arising purely from thermal effects, one can either use non-equilibrium QFT or thermal field theory outlined above. However, it was originally noticed that there exists a discrepancy between these two approaches in calculating the CP asymmetry in seesaw-based leptogenesis, *i.e.*, in the heavy Majorana neutrino decay $N\to H L$ [@Garny:2009rv; @Garny:2009qn]. In ref. [@Garny:2010nj], on the other hand, it was demonstrated that both approaches can yield the same result for this process, if the conventional time-ordered (TO) cut [@Kobes:1986za; @Kobes:1990ua; @Gelis:1997zv; @Giudice:2003jh] is replaced by the retarded/advanced product (dubbed as retarded/advanced (RA) cut hereafter for comparison) [@Kobes:1990kr; @Kobes:1990ua]. In this paper, we shall stick to the thermal real-time formalism, but take the subtlety into account by using both the TO and RA cuts.

![\[thermalcuts\] Thermal TO cuts (circlings) for producing a non-vanishing CP asymmetry in $\Phi\to L_i \bar\nu_{j}$ decay. The external thermal indices are fixed to $+$, while the internal index is summed over $\pm$.](figures/CPasym-Cut-1.pdf "fig:"){width="45.00000%"} ![\[thermalcuts\] Thermal TO cuts (circlings) for producing a non-vanishing CP asymmetry in $\Phi\to L_i \bar\nu_{j}$ decay. The external thermal indices are fixed to $+$, while the internal index is summed over $\pm$.](figures/CPasym-Cut-2.pdf "fig:"){width="45.00000%"}

With our definitions of the amplitudes $I_{0,1}$, the imaginary (absorptive) part of the product $I_0^* I_1$ in eq.  can be written as $$\begin{aligned}
\label{imagianry amp}
\text{Im}(I_0^* I_1)=\frac{1}{2i}I_0^* \sum_{\text{circling}}I_1.\end{aligned}$$ With the TO-cutting scheme, there are two circling diagrams contributing to the CP asymmetry, as shown in figure \[thermalcuts\]. Summing the internal thermal index over $\pm$, while fixing the external indices to $+$, we can write the corresponding amplitudes as $$\begin{aligned}
I_1^{(a)} &= i\int\frac{d^4 k}{(2\pi)^4} \left(\bar{u}_{L_i} P_R\, u_{e_l} \right)\left(\bar{u}_{e_l} P_L u_{L_k} \right)\left(\bar{u}_{L_k} P_R v_{\nu_j} \right)\,\nonumber\\[0.1cm]
&\hspace{0.5cm} \times G^F_{++}(p_i)\,G^F_{+-}(k)\,G^B_{+-}(p_i-k),
\\[0.2cm]
I_1^{(b)} & = -i\int\frac{d^4 k}{(2\pi)^4} \left(\bar{u}_{L_i} P_R\, u_{e_l} \right)\left(\bar{u}_{e_l} P_L u_{L_k} \right)\left(\bar{u}_{L_k} P_R v_{\nu_j} \right)\,\nonumber\\[0.1cm]
&\hspace{0.5cm} \times G^F_{--}(p_i)\,G^F_{-+}(k)\,G^B_{-+}(p_i-k),\end{aligned}$$ where the superscripts $B$ and $F$ denote the bosonic and fermionic propagators, respectively. The absorptive part for the TO cut is then determined to be $$\begin{aligned}
\label{ImIoI1}
\text{Im}(I_0^* I_1)^{ \text{TO}} &= \frac{1}{2 (2\pi)^2}\int d\omega\; \vert\boldsymbol{k}\vert^2 d \vert\boldsymbol{k}\vert\; d\cos\theta \; d\varphi \times \text{Tr}
\nonumber \\[0.1cm]
&\times \frac{1}{p_i^2-m_{L_k}^2} \times \delta[k^2-m_{e_l}^2]\times \delta[(p_i-k)^2-M_H^2]
\nonumber \\[0.1cm]
&\times \Big\{ \left[\theta(-\omega)-f_F(\vert \omega\vert)\right]\cdot \left[\theta(-(E_i-\omega))+f_B(\vert E_i-\omega\vert)\right]
\nonumber \\[0.1cm]
&\hspace{0.4cm} + \left[\theta(\omega)-f_F(\vert \omega\vert)\right]\cdot\left[\theta(E_i-\omega)+f_B(\vert E_i-\omega\vert)\right] \Big\},\end{aligned}$$ where the four-momenta $k$ and $p_i$ are decomposed, respectively, as $k=(\omega,\boldsymbol{k})$ and $p_i=(E_i, \boldsymbol{p_i})$, while $\cos\theta\equiv \boldsymbol{p_i}\cdot \boldsymbol{k}/\vert\boldsymbol{p_i}\vert \vert\boldsymbol{k}\vert$. The trace from spin summation is given by $$\begin{aligned}
\text{Tr}=(k\cdot p_i) \left(4 q\cdot p_i-2m_{L_i}^2\right)-2m_{L_i}^2 (k\cdot q).\end{aligned}$$

To perform the integration in eq. , a convenient way is to integrate firstly over $\cos\theta$ via $\delta[(p_i-k)^2-M_H^2]$, then over $\vert\boldsymbol{k}\vert$ via $\delta[k^2-m_{e_l}^2]$, and finally over $\omega$. Due to the appearance of Heaviside step functions in eq. , however, we must determine the sign of $\omega$ before performing the integration over $\omega$. To this end, remembering that $M_{H}$ is much larger than $m_{L_i}$ and $m_{e_l}$, the presence of $\delta[(p_i-k)^2-M_H^2]$, together with $\delta[k^2-m_{e_l}^2]$, implies that $$\begin{aligned}
\label{deltamil}
\Delta m_{il}^2\equiv m_{L_i}^2+m_{e_l}^2-M_{H}^2=2p_i\cdot k=2(E_i \omega -\vert\boldsymbol{p_i}\vert \vert\boldsymbol{k}\vert \cos\theta)<0. \end{aligned}$$ With $-1\leqslant\cos\theta \leqslant 1$, it can then be found that $$\begin{aligned}
\label{oemga sign}
\omega<0,\quad  E_i-\omega>0. \end{aligned}$$ As a consequence, the overall dependence of eq.  on the distribution functions is now simplified as $$\begin{aligned}
\label{quadratic dependence}
N(\omega)\equiv f_B(\vert E_i-\omega\vert)-f_F(\vert \omega\vert)-2f_F(\vert \omega \vert)f_B(\vert E_i-\omega\vert).\end{aligned}$$ The final integration region of $\omega$ is determined by $$\begin{aligned}
-1\leqslant\frac{\Delta m_{il}^2-2E_i \omega}{-2\vert\boldsymbol{p_i}\vert \vert\boldsymbol{k}\vert}\leqslant 1,\end{aligned}$$ where $\vert\boldsymbol{k}\vert $ takes the approximate dispersion relation, $\vert\boldsymbol{k}\vert =\sqrt{\omega^2-m_{e_l}^2}$, resulting therefore in $\omega_{min}\leqslant\omega\leqslant\omega_{max}$, with $$\begin{aligned}
\omega_{min(max)}=\frac{1}{4\,M_{\Phi}\,m_{L_i}^2}\left[\Delta m_{il}^2\, (M_{\Phi}^2+m_{L_i}^2)\mp (M_\Phi^2-m_{L_i}^2)\,\sqrt{\Delta m_{il}^4-4m_{e_l}^2 m_{L_i}^2}\,\right],\end{aligned}$$ in the limit of vanishing neutrino masses. Our final expression of the CP asymmetry is then given by $$\begin{aligned}
\label{CP result}
\epsilon_D = -2\,\frac{1}{\sum\limits_i (Y_\nu Y^\dagger_\nu)_{ii}(M_\Phi^2-m_{L_i}^2)}
\sum\limits_{i\neq k}\text{Im}[(Y_\nu Y_\nu^\dagger)_{ki}(Y_{\ell,il}Y^\dagger_{\ell,lk})]\, \mathcal{F}(M_\Phi^2,m_{L_i}^2,m_{L_k}^2,m_{e_l}^2),\end{aligned}$$ where the scalar function is defined as $$\begin{aligned}
\label{scalar fun}
\mathcal{F}(M_\Phi^2,m_{L_i}^2,m_{L_k}^2,m_{e_l}^2)&=\frac{1}{8\pi}\frac{M_\Phi^2}{(M_\Phi^2-m_{L_i}^2)(m_{L_i}^2-m_{L_k}^2)}
\nonumber \\[0.1cm]
&\times \int_{\omega_{min}}^{\omega_{max}} d\omega \left(\Delta m_{il}^2 M_\Phi-2m_{L_i}^2\omega\right) N(\omega).\end{aligned}$$ It should be emphasized here that the dependence of $\mathcal{F}$ on the index $l$ comes from the charged-lepton Yukawa coupling contribution to the thermal lepton mass $m_{e_l}^2$ present in $\Delta m_{il}^2$ (see eq. ). As $\Delta m_{il}^2$ is dominated by contributions from the gauge and top-quark Yukawa couplings (the thermal masses will be discussed later in section \[sec:4\]), it can be inferred that the $l$ dependence is very weak, rendering therefore the CP asymmetry to carry an imaginary piece, $\text{Im}[(Y_\nu Y_\nu^\dagger)_{ki}(Y_{\ell}Y^\dagger_{\ell})_{ik}]$, as mentioned already in the Introduction. Furthermore, the $i$ dependence coming from $M_\Phi^2-m_{L_i}^2$ is also overwhelmed by contributions from the gauge couplings, as well as the possibly sizable scalar potential parameters and quark Yukawa couplings. However, as the dominate contributions from the gauge couplings are canceled out in $m_{L_i}^2-m_{L_k}^2$, the $i,k$ dependence coming from $m_{L_i}^2-m_{L_k}^2$ cannot be neglected. Thus, the CP asymmetry given by eq.  displays a nontrivial dependence on the indices $i,k$.

Using the TO cut, we have obtained a quadratic dependence of the CP asymmetry on the distribution functions, as can be seen from eq. . Such a quadratic dependence was also derived for the thermal Higgs decay $H\to NL$ in ref. [@Giudice:2003jh]. Within the non-equilibrium QFT framework [@Frossard:2012pc], however, the dependence was found to be linear in $f_H+f_L$ for the same decay [@Hambye:2016sby]. Following the argument made in ref. [@Garny:2010nj], we now turn to use the RA cut to determine the absorptive part, and check if such a linear dependence can be reproduced in our case. With our convention, the imaginary amplitude is given by $$\begin{aligned}
\text{Im}(I_0^*I_1)^{\text{R/A}}=\mp \frac{1}{2i}I_0^* \sum^{RA\;\text{cut}}_{\text{circling}}I_1,\end{aligned}$$ where $\mp$ correspond to the results with the retarded/advanced cut, respectively. In this context, only the circling diagram shown in figure \[thermalcuts\](a) contributes, leading to $$\begin{aligned}
\label{RAcut_ImI0I1}
\text{Im}(I_0^*I_1)^{\text{R/A}}&=\mp \frac{1}{2i}\int \frac{d^4 k}{(2\pi)^4}  \left(\bar{u}_{L_i}  P_R v_{\nu_j} \right)^* \left(\bar{u}_{L_i}  P_R u_{e_l} \right)\left(\bar{u}_{e_l}  P_L u_{L_k} \right)\left(\bar{u}_{L_k}   P_R v_{\nu_j} \right) \nonumber \\[0.1cm]
&\hspace{0.5cm} \times \left[D_{L_k}(p_i)D^+_{e_l}(k)D_{H}^+(p_i-k)-D_{L_k}(p_i)D^-_{e_l}(k)D_{H}^-(p_i-k)\right] \nonumber \\[0.2cm]
&=\pm\frac{1}{2(2\pi)^2}\int d^4k \frac{1}{m_{L_i}^2-m_{L_k}^2}\delta[k^2-m_{e_l}^2]\delta[(p_i-k)^2-M_H^2] \times \text{Tr} \nonumber \\[0.1cm]
&\hspace{0.5cm} \times \left[f_F(-\omega)+f_B(E_i-\omega)\right],\end{aligned}$$ where the thermal propagators with the RA cut are given by $$\begin{aligned}
 D(p)=G_{++}(p), \quad D^-(p)=G_{+-}(p), \quad D^+(p)=G_{-+}(p),\end{aligned}$$ and eq.  has been used. It can be clearly seen from eq.  that the RA-cutting scheme does lead to a linear dependence on the distribution function $f_B+f_F$. At the same time, the imaginary amplitude, $\text{Im}(I_0^*I_1)^{\text{TO}}$, obtained with the TO-cutting scheme can be reconciled to match the retarded amplitude, $\text{Im}(I_0^*I_1)^{\text{R}}$, via the following replacement for the distribution functions: $$\begin{aligned}
\label{distribution dependence}
f_B-f_F-2f_B f_F\to f_B+f_F. \end{aligned}$$ It should be mentioned that the linear dependence obtained in ref. [@Hambye:2016sby] is also based on a retarded self-energy cut, though within the non-equilibrium QFT framework.

In conclusion, using the real-time formalism in thermal field theory, a quadratic dependence on the distribution functions is obtained under the conventional TO-cutting scheme, while a linear dependence appears in the RA-cutting scheme. It is also found that, albeit with a different dependence on the distribution functions, the retarded amplitude can be simply obtained from the TO result with the replacement specified by eq. , which has also been observed in ref. [@Garny:2010nj].

### Simplified Boltzmann equation: freeze-in evolution

With the CP asymmetry in hand, we now proceed to determine the evolution of the lepton-number asymmetry based on the Boltzmann equation. The general Boltzmann equation for species $X$ participating in the process $A+B\rightleftarrows C+X$ reads $$\begin{aligned}
\label{Boltzmann eq}
\dot{n}_X+3H n_X & = \int d\Pi_X\, d\Pi_A\, d\Pi_B\, d\Pi_C\, (2\pi)^4\,\delta^{(4)}(p_A+p_B-p_C-p_X) \\[0.1cm]
&\hspace{-2.0cm} \times \Big[\vert\mathcal{M}_{A+B\to C+X}\vert^2\,f_A f_B (1\pm f_C)(1\pm f_X) -\vert \mathcal{M}_{C+X\to A+B}\vert^2 f_C f_X (1\pm f_A)(1\pm f_B) \Big],\nonumber\end{aligned}$$ where the phase-space factor is given by $d\Pi_i=d^3p_i/(2\pi)^32E_i$. The Dirac delta function $\delta^{(4)}(p_A+p_B-p_C-p_X)$ enforces the four-momentum conservation in collisions. The amplitude squared is obtained by summing over the initial- and final-state spins but without average. The factors $1\pm f_i$ correspond to the Bose enhancement and the Pauli blocking effect, respectively. The Hubble parameter at radiation-dominated flat Universe is given by $H=1.66\sqrt{g^\rho_\ast}\,T^2/M_{Pl}$, where $g^\rho_*$ denotes the relativistic d.o.f at temperature $T$, and $M_{Pl}=1.2\times 10^{19}$ GeV is the Planck mass.

In the early Universe, right-handed neutrinos are produced only through feeble Yukawa interactions. Thus, the production is essentially out-of-equilibrium and the particle number density is therefore negligibly small. Such a freeze-in production mechanism [@Hall:2009bx] effectively prevents large washout effects in the neutrino number density (as well as the lepton-number asymmetry generated therein) from inverse decay and annihilation scattering. In this context, the Boltzmann equation for the lepton-number asymmetry generated in $\Phi\to L \bar \nu$ decay can be much simplified by neglecting the inverse decay and annihilation scattering, because these processes are proportional to the negligible particle-number density. As a consequence, the lepton-number asymmetry can be accumulated as the right-handed neutrinos are produced and converted to the baryon-number asymmetry via rapid sphaleron transitions. Since the lepton-number asymmetry stored in the right-handed neutrinos is equal but with an opposite sign to that in the lepton doublets, we can determine it in either sector. For the lepton doublet, the Boltzmann equation can be simplified as $$\begin{aligned}
\dot{n}_L+3H n_L =\int d\Pi_{\Phi} f^{eq}_\Phi \int d\Pi_{\nu} d\Pi_L\,  (2\pi)^4\delta^{(4)}(p_\Phi-p_{\nu}-p_L)\,\vert\mathcal{M}(\Phi\to L \bar\nu)\vert^2,\end{aligned}$$ where, as an approximation, we have set the quantum statistic factors $1\pm f\approx1$. In addition, as we consider a thermal particle with vacuum mass around $\mathcal{O}(10^2)$ GeV, for a simple estimation, we shall use $f^{eq}_{\Phi}=e^{-E/T}$ for the phase-space integration.

With our definition of the CP asymmetry $\epsilon_D$ (see eq. ), the evolution of the lepton-number asymmetry $n_{\Delta L}\equiv n_L-n_{\bar L}$ can be written as $$\begin{aligned}
\label{lepton-asymmetry Boltzmann }
\dot{n}_{\Delta L}+3Hn_{\Delta L}&\approx \int d\Pi_\Phi f_{\Phi}^{eq}\, 2\,g_{\Phi} \,M_{\Phi}\,\left[\Gamma(\Phi \to L \bar\nu)-\Gamma(\bar\Phi\to \bar L \nu)\right] \nonumber\\[0.1cm]
&=\int d\Pi_{\Phi}\,f_{\Phi}^{eq}  \times 2\,g_{\Phi}\,M_{\Phi} \times 2\,\epsilon_D \times \Gamma(\Phi\to L \bar\nu),\end{aligned}$$ where $g_{\Phi}=4$ accounts for the four d.o.f of the Higgs doublet $\Phi$. In the sphaleron-active epoch, the relativistic d.o.f $g_*^s$ present in the entropy density, $s=T^3 g_*^s 2\pi^2/45$, can be treated as a constant, and thus we can use the yield definition $Y=n/s$, with $\dot{s}=-3Hs$, and $\dot{T}=-H T$, to rewrite the above equation as $$\begin{aligned}
\label{lepton-asymmetry Boltzmann2}
\frac{dY_{\Delta L}}{dT}=-\frac{1}{sH}\,\frac{1}{\pi^2}\,g_{\Phi}\,\epsilon_D\,M_{\Phi}^2\,  K_1(M_{\Phi}/T)\,\Gamma(\Phi\to L\bar\nu).\end{aligned}$$ Here $K_1$ denotes the first modified Bessel function of the second kind. In the rest frame of $\Phi$, the decay width is given by $$\begin{aligned}
\label{decay rate}
\Gamma(\Phi \to L \bar\nu)&=\sum\limits_{i=e,\mu,\tau}\frac{1}{8\pi}\,(Y_\nu Y_\nu^\dagger)_{ii}\,\frac{\vert \boldsymbol{p_i}\vert}{M_{\Phi}^2}\,(M_{\Phi}^2-m_{Li}^2),\end{aligned}$$ where $\vert\boldsymbol{p_i}\vert=(M_{\Phi}^2-m_{Li}^2)/(2M_{\Phi})$ is the momentum of the two final-state particles, and the neutrino masses have been neglected.

Up to now, we have obtained a nonzero kinetic phase contained in $\text{Im}(I_0^*I_1)$ within the framework of thermal field theory. To generate a non-vanishing lepton-number asymmetry, however, non-diagonal Yukawa couplings $Y_{\ell,\nu}$ are also required. In the next section, we turn to exploit the Yukawa structures that can generate a nonzero coupling phase contained in $\text{Im}(c_0^* c_1)$ and, at the same time, produce compatible neutrino oscillation data.

Non-diagonal Yukawa textures in the lepton sector {#sec:3}
=================================================

When the lepton-number asymmetry is generated via the freeze-in production of right-handed neutrinos, the washout effects from inverse decay and annihilation scattering are negligible. Therefore, the flavor effects encoded in these washout processes, *e.g.*, $\Phi \bar\nu_{i}\rightleftharpoons \bar \Phi \nu_{j}$, would not play a significant role in the Boltzmann evolution. However, as mentioned already below eq. , the dependence of $\mathcal{F}$ on the indices $i, k$ from the lepton-doublet propagator cannot be neglected. In fact, such a dependence is crucial to induce a non-vanishing CP asymmetry $\epsilon_D$, because otherwise the imaginary coupling sector would vanish, *i.e.*, $\text{Im}[\text{tr}(Y_\nu Y^\dagger_\nu\,Y_\ell Y^\dagger_\ell)]=0$. Furthermore, a nonzero $\epsilon_D$ also requires the Yukawa matrices $Y_{\ell,\nu}$ to be non-diagonal. Therefore, after summing over the lepton flavors, the CP asymmetry is still texture dependent, which is a generic feature of leptogenesis [@Branco:2002kt].

As a consequence, the freeze-in DN considered here essentially puts us towards the flavor puzzle in particle physics, on which no consensus has been hitherto reached. Especially, it is not known *a priori* whether a flavor basis, in which the charged-lepton (or neutrino) Yukawa matrix is diagonal while the other is not, should be used as a natural setup, even though the charged-lepton flavors are often assumed to populate in the diagonal space for most flavor studies. On the other hand, the SM is often extended by introducing sufficient dynamical fields and free parameters, which might break the freeze-in DN or result in the BAU explanation via other avenues. In order not to spoil the formulation and results obtained up to now in this paper, we shall consider the lepton Yukawa textures from a phenomenological point of view. With such a perspective in mind, we make a general assumption: the non-trivial Yukawa textures are induced by some underlying flavor mechanism, such that it forbids us to choose the otherwise arbitrary flavor basis, in particular, the basis in which either the charged-lepton or the neutrino Yukawa matrix is diagonal. This can be realized, *e.g.*, by some symmetry-based ansatz in which a preferred basis is constructed under the symmetry invariance. In addition, we shall postulate that the observed pattern of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [@Pontecorvo:1957cp; @Maki:1962mu] is induced by the TB mixing with a minimal correction from the charged-lepton or neutrino sector [@Albright:2008rp].

The well-known TB mixing pattern has a mass-independent form [@Harrison:2002er] $$\begin{aligned}
V_{TB} =\left(
\begin{array}{ccc}
 \sqrt{\frac{2}{3}}\,\,  &  \frac{1}{\sqrt{3}}\,\,   &                   0 \\[0.2cm]
-\frac{1}{\sqrt{6}}\,\,  &  \frac{1}{\sqrt{3}}\,\,   &  \frac{1}{\sqrt{2}} \\[0.2cm]
-\frac{1}{\sqrt{6}}\,\,  &  \frac{1}{\sqrt{3}}\,\,   & -\frac{1}{\sqrt{2}} \\
\end{array}
\right).\end{aligned}$$ Although the pure TB mixing is already excluded by the observed nonzero reactor angle $\theta_{13}$ (see *e.g.*, refs. [@Esteban:2016qun; @deSalas:2017kay; @Esteban:2018azc; @Capozzi:2018ubv] for recent global analyses), it has been pointed out that, with a minimal correction from the charged-lepton or neutrino sector, this mixing pattern can readily produce the compatible neutrino oscillation data, while retaining the predictability and testability of the relations among the mixing angles [@Albright:2008rp; @He:2011gb]. As demonstrated in ref. [@He:2011gb], there exist four possible minimal corrections to the TB mixing pattern that are still compatible with the current PMNS data at $3\sigma$ level. According to which column or row of the TB matrix is invariant under the minimal corrections, the modified patterns can be classified as $\text{TM}_i$ (invariance of the $i$-th column) and $\text{TM}^i$ (invariance of the $i$-th row) [@Albright:2008rp]. Thus, on account of the observations made in refs. [@Albright:2008rp; @He:2011gb], we have $$\begin{aligned}
\label{four patterns}
\text{TM}_1:\;\;U&=V_{TB}R_{23},\qquad
\text{TM}_2:\;\;U=V_{TB}R_{13},
\nonumber \\[0.2cm]
\text{TM}^2:\;\;U&=R_{13}V_{TB},\qquad
\text{TM}^3:\;\;U=R_{12}V_{TB},\end{aligned}$$ with the unitary Euler rotation matrices given, respectively, by $$\begin{aligned}
R_{12}(\theta)&=\left(
\begin{array}{ccc}
\cos\theta\; & \sin\theta e^{i\varphi}\; & 0\\[0.1cm]
-\sin\theta e^{-i\varphi }\; &  \cos\theta\; & 0\\[0.1cm]
0\; & 0\; & 1\\
\end{array}
\right),\quad
R_{13}(\theta)=\left(
\begin{array}{ccc}
\cos\theta\; & 0\;    & \sin\theta e^{i\varphi }\\[0.1cm]
0\; & 1\;   & 0\\[0.1cm]
-\sin\theta e^{-i\varphi}\; &   0\;  &  \cos\theta \\
\end{array}
\right),
\nonumber \\[0.2cm] 
R_{23}(\theta )&=\left(
\begin{array}{ccc}
1\;&                   0\;           &                 0\\[0.1cm]
0\;&  \cos\theta\;                   &  \sin\theta e^{i\varphi}\\[0.1cm]
0\;& -\sin\theta e^{-i\varphi }\;    &  \cos\theta\\
\end{array}
\right),\end{aligned}$$ where $0\leq\theta\leq\pi$ and $0\leq\varphi\leq2\pi$. With the convention $M_f=V_f^{L\dagger} \hat{M}_f V_f^R$, where $M_f$ and $\hat{M}_f$ represent the mass matrix before and after the diagonalization, the flavor (primed) and the mass (un-primed) eigenstates are transformed to each other via the relations $f_{L(R)}^\prime =V_f^{L(R)\dagger} f_{L(R)}$, and the PMNS matrix present in the $W$-mediated charged current $\bar{\ell}_L U \gamma_\mu \nu_L$ is given by $U=V_\ell^L V_\nu ^{L \dagger}$. Then, in accordance with eq. , the charged-lepton (neutrino) mixing matrix in patterns $\text{TM}_{1,2}$ would be given by $V_\ell^L=V_{TB}$ ($V_\nu^L=R_{23}^\dagger, R_{13}^\dagger$), while the matrix for neutrinos (charged leptons) turns out to be $V_\nu^L=V_{TB}^\dagger$ ($V_\ell^L=R_{13}, R_{12}$) in patterns $\text{TM}^{2,3}$. The product of Yukawa matrices can also be rewritten in terms of the mixing and physical mass matrices as $$\begin{aligned}
Y_f Y_f^\dagger =\frac{2}{v_f^2}V_f^{L \dagger} \hat{M}_f^2 V^L_f,\end{aligned}$$ where $v_f$ denote the VEVs developed by the Higgs doublets responsible for generating the charged-lepton and neutrino masses, respectively. Given that the leptonic CP asymmetry $\epsilon_D$ is approximately proportional to $\text{Im}[(Y_\nu Y_\nu^\dagger)_{ki}(Y_{\ell}Y^\dagger_{\ell})_{ik}]$ (see eq. ), it can be seen that, after fixing the kinematics, $\epsilon_D$ will depend on the two mixing parameters, $(\theta, \varphi)$, both of which are also directly responsible for producing the current neutrino oscillation data.

Besides the requirement that eq.  should produce the observed moduli of the PMNS matrix, $\vert U \vert$, a basis-independent and rephasing-invariant measure of the low-energy CP violation, defined as [@Jarlskog:1985ht] $$\begin{aligned}
\mathcal{J}\sum\limits_\gamma \epsilon_{\alpha \beta \gamma} \sum\limits_k \epsilon_{ijk}=\text{Im}[U_{\alpha i}U^*_{\alpha j}U^*_{\beta i}U_{\beta j}],\end{aligned}$$ is also crucial to exploit how successfully the freeze-in DN can be inferred, particularly, from the sign of $\mathcal{J}$. In the standard convention of the PMNS matrix [@Tanabashi:2018oca; @Xing:2019vks], the Jarlskog invariant $\mathcal{J}$ is given by $$\begin{aligned}
\mathcal{J}=\frac{1}{8}\cos\theta_{13}\sin(2\theta_{12})\sin(2\theta_{23}) \sin(2\theta_{13})\sin\delta.\end{aligned}$$ Recently, a global fit of neutrino oscillation parameters has obtained a strong preference for values of the Dirac CP phase $\delta$ in the range $[\pi,2\pi]$ [@deSalas:2017kay]. Given the observed mixing angles, this implies a negative CP measure $\mathcal{J}<0$. In particular, the best-fit value favors $\delta\simeq 3\pi/2$ [@Esteban:2016qun], leading to (at $3\sigma$ level) $$\begin{aligned}
\mathcal{J}^{max}_{CP}=-\left(0.0329^{+0.0021}_{-0.0024}\right),\end{aligned}$$ where the uncertainties come from the determination of the mixing angles.

Corresponding to the four mixing patterns specified by eq. , the Jarlskog invariant $\mathcal{J}$ is given, respectively, by $$\begin{aligned}
\label{Jarlskog}
\text{TM}_1&:~\mathcal{J}=-\frac{\sin(2\theta)\sin\varphi}{6\sqrt{6}}, \qquad &\text{TM}_2&:~\mathcal{J}=-\frac{\sin(2\theta)\sin\varphi}{6\sqrt{3}},
\nonumber \\[0.5mm]
\text{TM}^2&:~\mathcal{J}=\frac{\sin(2\theta)\sin\varphi}{12}, \qquad &\text{TM}^3&:~\mathcal{J}=-\frac{\sin(2\theta)\sin\varphi}{12}.\end{aligned}$$ If the maximal CP violation, $\mathcal{J}=\mathcal{J}_{CP}^{max}$ is assumed, and the values of the mixing angles, $(\theta,\varphi)$ are taken to produce the $3\sigma$ ranges of the PMNS matrix moduli $\vert U\vert$, we can then establish whether such a low-energy maximal CP violation can prompt a successful DN. To this end, we need firstly specify the decaying particle for generating the leptonic CP asymmetry, which will be explored at length in the next section.

Thermal scalar implementation {#sec:4}
=============================

As an alternative to most DN applications in which the lepton-number asymmetry is generated by non-thermal heavy particle decays [@Dick:1999je; @Murayama:2002je; @Cerdeno:2006ha; @Gu:2007mc; @Bechinger:2009qk; @Narendra:2017uxl], we have considered the case where the asymmetry is accumulated via the freeze-in production of right-handed neutrinos from thermal scalar decay. In this section, we shall specify the minimal Higgs doublet for implementing the freeze-in DN described in section \[sec:2\], and consider, in particular, the CP asymmetry obtained by using both the TO and RA cuts.

SM Higgs case
-------------

If the DN were realized by the SM Higgs, the neutrino mass and BAU would then be simultaneously addressed by simply adding the missing neutrino Yukawa interactions (eq. ) to the SM. The only price to pay is to accept the non-aesthetic, feeble neutrino Yukawa couplings, which are of $\mathcal{O}(10^{-14})$ for $\mathcal{O}(10^{-2})$ eV neutrino masses.

Since the sphaleron-active epoch, $10^2~\text{GeV}<T<10^{12}~\text{GeV}$, is considered, we shall use the thermal masses of the SM Higgs and leptons that are given by [@Weldon:1982bn; @Giudice:2003jh] $$\begin{aligned}
M_H^2 &\simeq \left(\frac{3}{16}g_2^2+\frac{1}{16}g_1^2+\frac{1}{4}y_t^2+\frac{1}{2}\lambda\right)T^2,  \label{thermal Higgs masses}
\\[0.2cm]
m_{L_i}^2 &= \left(\frac{3}{32}g_2^2+\frac{1}{32}g_1^2+\frac{1}{16}(Y_{\ell}Y^\dagger_{\ell})_{ii}\right)T^2,  \label{thermal L-doublet masses}
\\[0.2cm]
m_{e_l}^2 &= \left(\frac{1}{8}g_1^2+\frac{1}{8}(Y_{\ell}Y^\dagger_{\ell})_{ll}\right)T^2,
\label{thermal L-singlet masses}\end{aligned}$$ where $g_{2}$ ($g_{1}$) is the $\text{SU(2)}_L$ ($\text{U(1)}_Y$) gauge coupling, and $\lambda$ is the SM Higgs potential parameter satisfying the tadpole equation, $m_h^2=\lambda v^2$, with $h$ being the physical SM Higgs boson. We have only kept the dominant top-Yukawa contribution to $M_{H}^2$ by assuming a diagonal Yukawa matrix in the up-type quark sector. A general charged-lepton Yukawa matrix is, however, retained for the thermal masses of leptons, because a non-diagonal $Y_\ell$ is crucial for generating a non-vanishing lepton-number asymmetry. While the renormalization-group running of the coupling constants should in principle be taken into account at a scale $\mu\simeq 2\pi T$ [@Giudice:2003jh], which would prompt additional $T$-dependent sources, as a simple estimation, we shall use here the vacuum values of these coupling constants.

Based on the analysis of index dependence of the CP asymmetry made below eq. , we can neglect the small Yukawa contributions to $M_\Phi^2 \pm m_{L_i}^2$ and $\Delta m_{il}^2$, because they are basically overwhelmed by the contributions from gauge, potential parameters, and top-quark Yukawa couplings. Under this approximation, the integration of eq.  over the sphaleron-active regime induces a semi-analytic expression for the baryon-number asymmetry. It is found that, for the mixing patterns $\text{TM}_{1,2}$, the baryon-number asymmetry is estimated to be $$\begin{aligned}
\label{asymmetry in SM Higgs}
\text{TM}_1&:\;\;Y_{\Delta B}^{\text{TO}}\simeq -Y_{\Delta B}^{\text{R}} \simeq  \mathcal{O}(10^{-15}) \sin(2\theta) \sin\varphi,  \\[0.2cm]
\text{TM}_2&:\;\;Y_{\Delta B}^{\text{TO}}\simeq -Y_{\Delta B}^{\text{R}}\simeq  -\mathcal{O}(10^{-16}) \sin(2\theta)  \sin\varphi,\end{aligned}$$ where $Y_{\Delta B}^{\text{TO}}$ and $Y_{\Delta B}^{\text{R}}$ are obtained with the input of the CP asymmetry determined in the TO- and the retarded-cutting scheme, respectively. For patterns $\text{TM}^{2,3}$, the baryon-number asymmetry is found to be even smaller, with $$\begin{aligned}
\text{TM}^2&:\;\;Y_{\Delta B}^{\text{TO}}\simeq -Y_{\Delta B}^{\text{R}} \simeq  -\mathcal{O}(10^{-17}) \tan(2\theta)\sin\varphi,
\\[0.2cm]
\text{TM}^3&:\;\;Y_{\Delta B}^{\text{TO}}\simeq -Y_{\Delta B}^{\text{R}}\simeq   -\mathcal{O}(10^{-17}) \tan(2\theta)\sin\varphi. \label{asymmetry in SM Higgs2}\end{aligned}$$ To obtain the numerical factors, we have used $g_*^\rho=g_*^s=106.75$. In addition, we have adopted a normal-ordering neutrino mass hierarchy, as suggested by the recently global analyses [@Capozzi:2018ubv; @Esteban:2018azc], and neglected the lightest neutrino mass. Explicitly, the input values of neutrino masses are given by $m_1\simeq 0$, $m_2\simeq \sqrt{\Delta m_{21}^2}$, and $m_3\simeq \sqrt{\Delta m_{31}^2}$, with the mass-squared differences taken from ref. [@Esteban:2016qun]. It can be seen from eqs. – that the dependence of $Y_{\Delta B}$ on the trigonometric functions is different between $\text{TM}_{1,2}$ and $\text{TM}^{2,3}$. This is because an additional $\theta$ dependence appears in the thermal fermion masses for patterns $\text{TM}^{2,3}$. It is also found that the baryon-number asymmetries induced by the TO- and retarded-cutting CP asymmetries have basically the same size but with an opposite sign. As will be discussed in the next subsection, such a sign difference becomes important in generating a positive baryon-number asymmetry, together with a negative CP measure in neutrino oscillations.

Compared with the observed baryon-number asymmetry of the Universe at present day [@Aghanim:2018eyx], $$\begin{aligned}
\label{BAU today}
 Y_{\Delta B}=(8.75\pm 0.23)\times 10^{-11},\end{aligned}$$ the amount of asymmetry induced by the minimal SM Higgs is negligible. Although we have followed here a phenomenological perspective, $Y_{\Delta B}^{\text{TO,R}}$ given by eqs. – are primarily controlled by the neutrino Yukawa couplings $Y_\nu\simeq \mathcal{O}(10^{-14})$, and thus the orders of magnitude estimated therein are quite reasonable. This can also be justified by noting that, even though the neutrino Yukawa couplings may be canceled in the imaginary coupling sector, the decay rate $\Gamma(H\to L \bar\nu)$ involves the couplings at $\mathcal{O}(Y_\nu^2)$. In addition, as the SM Higgs also couples to the right-handed charged leptons, an additional contribution to the leptonic CP asymmetry can be induced by the vertex correction. It is, however, expected that such an amount of asymmetry would be similar to that generated by the wavefunction correction, as no quasi-degenerate mass spectrum could resonantly enhance the latter within the SM. Based on these observations, the SM Higgs implementation should be therefore dismissed, and we are driven to consider new scalars beyond the minimal SM.

Neutrinophilic two-Higgs-doublet model {#sec:4.2}
--------------------------------------

A direct enhancement of the lepton-number asymmetry can be achieved by invoking a sufficiently large neutrino Yukawa coupling, while retaining the out-of-equilibrium condition. This can be realized by introducing another Higgs doublet which develops a smaller VEV. As a minimal extension of the SM, let us focus on the neutrinophilic 2HDM [@Davidson:2009ha], in which the second Higgs doublet couples neither to quarks nor to right-handed charged leptons.

In such a neutrinophilic 2HDM, both the right-handed Dirac neutrinos and the new Higgs doublet possess an additional $Z_2$ parity. The model Lagrangian has a soft $Z_2$-breaking scalar potential [@Davidson:2009ha] $$\begin{aligned}
V&=m_1^2 H_1^{\dagger}H_1 +m_2^2H_2^{\dagger} H_2 -(m_{12}^2H_1^{\dagger} H_2 + \text{h.c.})
+\frac{\lambda_1}{2}(H_1^{\dagger} H_1)^2+ \frac{\lambda_2}{2}(H_2^{\dagger}H_2)^2
\nonumber \\
&+\lambda_3(H_1^{\dagger} H_1) (H_2^{\dagger}H_2)
+\lambda_4(H_1^{\dagger} H_2) (H_2^{\dagger} H_1) + \left[ \dfrac{\lambda_5}{2}(H_1^{\dagger} H_2)^2+ \text{h.c.} \right].
\label{Higgs potential}\end{aligned}$$ For a real and positive soft-breaking term $m_{12}^2\ll v^2$, with $v_1^2+v_2^2=v^2=(246~\text{GeV})^2$, the tadpole equations, $\partial V/\partial H_i=0$, would induce a seesaw-like relation $$\begin{aligned}
v_1\simeq v, \qquad v_2\simeq \frac{m_{12}^2v}{ \lambda_{345} v^2+m_{2}^2},\end{aligned}$$ with $\lambda_{345}\equiv(\lambda_3+\lambda_4+\lambda_5)/2$. In the conventional neutrinophilic 2HDM [@Davidson:2009ha] (see also some phenomenological studies of the model performed in refs. [@Machado:2015sha; @Bertuzzo:2015ada]), the value of $v_2$ is tuned at eV scale so as to have $\mathcal{O}(1)$ neutrino Yukawa couplings. Apparently, when $Y_\nu\simeq \mathcal{O}(1)$, neutrinos would establish the L-R equilibration in the sphaleron-active epoch, and thus no net lepton-number asymmetry would be stored. Here we assume, instead, $Y_\nu\lesssim \mathcal{O}(10^{-8})$ to guarantee the out-of-equilibrium generation of the lepton-number asymmetry. Such an assumption is justified by the requirement $v_2\gtrsim \mathcal{O}(10^{-3})$ GeV, which in turn indicates that $m_{12}\gtrsim0.5$ GeV for $m_{2}\simeq \mathcal{O}(v)$.

At a temperature well above the electroweak scale, we use the following thermal mass for the second Higgs doublet [@Cline:1995dg]: $$\begin{aligned}
M_{H_2}^2\simeq \left(\frac{3}{16}g_2^2+\frac{1}{16}g_1^2+\frac{1}{4}\lambda_2+\frac{1}{6}\lambda_3+\frac{1}{12}\lambda_4\right)T^2,\end{aligned}$$ where marginal contributions from the soft $Z_2$-breaking term and the neutrino Yukawa couplings have been neglected. The thermal masses of $H_1$ and leptons are the same as that given by eqs. -. For $M_{H_2}$ present in eq. , we shall also include the positive mass parameter $m_2$, *i.e.*, $M_{H_2}^2=m_2^2+M_{H_2}^2(T)$. For our numerical analyses, we shall fix $m_2=500$ GeV, and assume that the effects from $\lambda_{4,5}$ are negligible. Furthermore, we shall work in the alignment limit where $H_1$ contains the SM Higgs boson. Within such a numerical setup, the Higgs mass spectrum is nearly degenerate, $M_{H^\pm}\simeq M_{H}\simeq M_{A}\simeq 500$ GeV, and $\lambda_1\simeq \lambda_2\simeq \lambda_3\simeq m_h^2/v^2$. For explicit expressions of the potential parameters $\lambda_{1-5}$ and the Higgs mass spectrum in the alignment limit, together with the theoretical and experimental constraints, the readers are referred to, *e.g.*, ref. [@Li:2018aov]. It can also be found in refs. [@Machado:2015sha; @Bertuzzo:2015ada] that such a numerical mass spectrum is phenomenologically viable.

As done in the last subsection, here the numerical integration of eq.  over the sphaleron-active regime also prompts a semi-analytic expression for the baryon-number asymmetry $Y_{\Delta B}$. Explicitly, for a TO-cutting CP asymmetry, we get $$\begin{aligned}
\text{TM}_1&:\; \;Y_{\Delta B}^{\text{TO}}\simeq 3.16\times 10^{-10}\;\frac{\sin(2\theta)\sin\varphi}{v_2^2},\label{baryon asym in pattern-1}
\\[0.2cm]
\text{TM}_2&: \;\; Y_{\Delta B}^{\text{TO}}\simeq -1.15\times 10^{-10}\;\frac{\sin(2\theta)\sin\varphi}{v_2^2},\label{baryon asym in pattern-2}
 \\[0.2cm]
\text{TM}^2&:\; \; Y_{\Delta B}^{\text{TO}}\simeq -1.33\times 10^{-12}\;\frac{\tan(2\theta)\sin\varphi}{v_2^2}, \label{baryon asym in pattern-3}
 \\[0.2cm]
\text{TM}^3&:\; \; Y_{\Delta B}^{\text{TO}}\simeq -1.33\times 10^{-12}\;\frac{\tan(2\theta)\sin\varphi}{v_2^2},\label{baryon asym in pattern-4}\end{aligned}$$ where $g_*^\rho=g_*^s=110.75$ has been used. Again, the different dependence of $Y_{\Delta B}^{\text{TO}}$ on the trigonometric functions between $\text{TM}_{1,2}$ and $\text{TM}^{2,3}$ is also due to an additional $\theta$ dependence of the thermal fermion masses in patterns $\text{TM}^{2,3}$.

![\[TO-TM-Angle\] Allowed regions for the two mixing parameters, $(\theta, \varphi)$, under the individual constraint from a positive baryon-number asymmetry (yellow bands), a negative CP measure $\mathcal{J}=\mathcal{J}_{CP}^{max}<0$ (red regions), as well as the PMNS matrix element $\vert U_{13}\vert$ (narrow blue bands) in neutrino oscillations [@Esteban:2016qun].](figures/TO-TMLow1_Angle_re.pdf "fig:"){width="42.00000%"} ![\[TO-TM-Angle\] Allowed regions for the two mixing parameters, $(\theta, \varphi)$, under the individual constraint from a positive baryon-number asymmetry (yellow bands), a negative CP measure $\mathcal{J}=\mathcal{J}_{CP}^{max}<0$ (red regions), as well as the PMNS matrix element $\vert U_{13}\vert$ (narrow blue bands) in neutrino oscillations [@Esteban:2016qun].](figures/TO-TMLow2_Angle_re.pdf "fig:"){width="42.00000%"}\
![\[TO-TM-Angle\] Allowed regions for the two mixing parameters, $(\theta, \varphi)$, under the individual constraint from a positive baryon-number asymmetry (yellow bands), a negative CP measure $\mathcal{J}=\mathcal{J}_{CP}^{max}<0$ (red regions), as well as the PMNS matrix element $\vert U_{13}\vert$ (narrow blue bands) in neutrino oscillations [@Esteban:2016qun].](figures/TO-TMUp2_Angle_re.pdf "fig:"){width="42.00000%"} ![\[TO-TM-Angle\] Allowed regions for the two mixing parameters, $(\theta, \varphi)$, under the individual constraint from a positive baryon-number asymmetry (yellow bands), a negative CP measure $\mathcal{J}=\mathcal{J}_{CP}^{max}<0$ (red regions), as well as the PMNS matrix element $\vert U_{13}\vert$ (narrow blue bands) in neutrino oscillations [@Esteban:2016qun].](figures/TO-TMUp3_Angle_re.pdf "fig:"){width="42.00000%"}

It can be seen from eqs.  and – that, for patterns $\text{TM}_{1,2}$, the baryon-number asymmetry $Y_{\Delta B}^{\text{TO}}$ has basically the same size but with an opposite sign, while the CP measure $\mathcal{J}$ has the same sign. This indicates that, to generate a positive $Y_{\Delta B}$, the product of the trigonometric functions, $\sin(2\theta)\sin\varphi$, should be positive (negative) for $\text{TM}_1$ ($\text{TM}_2$). However, if we follow the favored Dirac CP phase $\delta=[\pi,2\pi]$, which indicates a negative $\mathcal{J}$, the same factor $\sin(2\theta)\sin\varphi$ should be positive in both patterns. Therefore, for a successful DN with a TO-cutting CP asymmetry, the pattern $\text{TM}_2$ is already disfavored by the neutrino oscillation data with a Dirac CP phase in the range $\delta=[\pi,2\pi]$. For patterns $\text{TM}^{2,3}$, on the other hand, $Y_{\Delta B}^{\text{TO}}$ has basically the same value, while $\mathcal{J}$ has the opposite sign. This implies that the pattern $\text{TM}^3$ is also disfavored by the range of Dirac CP phase in realizing a successful DN. To visualize the sign significance observed above, we show in figure \[TO-TM-Angle\] the allowed regions for the two mixing parameters, $(\theta, \varphi)$, under the individual constraint from a positive baryon-number asymmetry, a negative CP measure in neutrino oscillations, as well as the PMNS matrix element $\vert U_{13}\vert$. For the two allowed patterns, $\text{TO}\text{--}\text{TM}_1$ and $\text{TO}\text{--}\text{TM}^2$, we further investigate in detail the compatibility between the freeze-in DN and the neutrino oscillation observables in a particular quadrant with $\theta=[0,\pi/2]$, which is shown in figure \[TO-TMLow1Up2\_Contour\].

![\[TO-TMLow1Up2\_Contour\] Compatibility between the freeze-in DN and the neutrino oscillation observables for patterns $\text{TM}_1$ and $\text{TM}^2$ with the TO-cutting scheme. The area enclosed by the black-dotted line represents the $3\sigma$ allowed range of $\vert U\vert$ and a maximal CP measure $\mathcal{J}=\mathcal{J}_{CP}^{max}<0$. The contours denote the variations of $v_2$ (in unit of GeV), with each contour corresponding to the best-fit point of $Y_{\Delta B}$ given by eq. .](figures/TO-TMLow1_Contour_re.pdf "fig:"){width="45.00000%"} ![\[TO-TMLow1Up2\_Contour\] Compatibility between the freeze-in DN and the neutrino oscillation observables for patterns $\text{TM}_1$ and $\text{TM}^2$ with the TO-cutting scheme. The area enclosed by the black-dotted line represents the $3\sigma$ allowed range of $\vert U\vert$ and a maximal CP measure $\mathcal{J}=\mathcal{J}_{CP}^{max}<0$. The contours denote the variations of $v_2$ (in unit of GeV), with each contour corresponding to the best-fit point of $Y_{\Delta B}$ given by eq. .](figures/TO-TMUp2_Contour_re.pdf "fig:"){width="45.00000%"}

With the input of a retarded-cutting CP asymmetry, it is found that, compared with $Y_{\Delta B}^{\text{TO}}$ given by eqs. -, $Y_{\Delta B}^{\text{R}}$ has basically the same size but with a different sign, as is observed already in the SM Higgs case. Using the same arguments as made above, the patterns $\text{TM}_1$ and $\text{TM}^2$ would be dismissed, while both $\text{TM}_2$ and $\text{TM}^3$ are favored by a Dirac CP phase in the range $\delta=[\pi,2\pi]$. The resulting $Y_{\Delta B}^{\text{R}}$ is now given by $$\begin{aligned}
\text{TM}_2&:\;\;Y_{\Delta B}^{\text{R}}\simeq 1.78\times 10^{-10}\,\frac{\sin(2\theta)  \sin\varphi}{v_2^2},
 \\[0.2cm]
\text{TM}^3&:\;\;Y_{\Delta B}^{\text{R}}\simeq 2.05 \times 10^{-12}\,\frac{\tan(2\theta) \sin\varphi}{v_2^2}.\end{aligned}$$ For the two allowed patterns $\text{R}\text{--}\text{TM}_2$ and $\text{R}\text{--}\text{TM}^3$, the trigonometric functions for prompting a positive baryon-number asymmetry and a negative CP measure $\mathcal{J}$ have the same behavior as observed for the pattern $\text{TO}\text{--}\text{TM}_1$ shown in figure \[TO-TM-Angle\], because these three patterns share a common sign in $Y_{\Delta B}$ and $\mathcal{J}$. The compatibility between the freeze-in DN and the neutrino oscillation observables are shown in figure \[R-TMLow2Up3\_Contour\], where a particular quadrant with $\theta=[0,\pi/2]$ is selected.

![\[R-TMLow2Up3\_Contour\] Same as in figure \[TO-TMLow1Up2\_Contour\] but for patterns $\text{TM}_2$ and $\text{TM}^3$ with the retarded-cutting scheme.](figures/R-TMLow2_Contour_re.pdf "fig:"){width="45.00000%"} ![\[R-TMLow2Up3\_Contour\] Same as in figure \[TO-TMLow1Up2\_Contour\] but for patterns $\text{TM}_2$ and $\text{TM}^3$ with the retarded-cutting scheme.](figures/R-TMUp3_Contour_re.pdf "fig:"){width="45.00000%"}

As shown in figure \[TO-TMLow1Up2\_Contour\], both $\text{TM}_1$ and $\text{TM}^2$ can produce the $3\sigma$-allowed range of $\vert U\vert$ as well as a maximal CP measure $\mathcal{J}=\mathcal{J}_{CP}^{max}<0$ (the area enclosed by the black-dotted line), pointing out a mixing angle in the range $0.2<\theta<0.3$. It is also observed that the maximal leptonic CP asymmetry at low energy favors a maximal CP phase $\varphi$ that is necessary for a lepton-number asymmetry at high-temperature regime: $\varphi\simeq \pi/2$ for $\text{TO}\text{--}\text{TM}_1$ and $\varphi\simeq 3\pi/2$ for $\text{TO}\text{--}\text{TM}^2$. For the two allowed patterns $\text{R}\text{--}\text{TM}_2$ and $\text{R}\text{--}\text{TM}^3$ shown in figure \[R-TMLow2Up3\_Contour\], on the other hand, the mixing angle is found at $\theta\simeq 0.2$, and a CP-violating phase $\varphi\simeq \pi/2$ is favored for both cases.

Finally, as can be seen from figures \[TO-TMLow1Up2\_Contour\] and \[R-TMLow2Up3\_Contour\], $v_2\simeq\mathcal{O}(0.1\text{--}1)$ GeV is required by the allowed range of $\vert U\vert$ and a maximal CP measure $\mathcal{J}=\mathcal{J}_{CP}^{max}<0$. With such a range of $v_2$ as input, the Dirac neutrino Yukawa couplings are then estimated to be $Y_\nu\simeq \mathcal{O}(10^{-10}\text{--}10^{-11})$ for neutrino masses at $\mathcal{O}(10^{-2})$ eV. Therefore, the feeble neutrino Yukawa couplings of $\mathcal{O}(10^{-10}\text{--}10^{-11})$ obtained in the neutrinophilic 2HDM can account for the smallness of Dirac neutrino masses in a simple while less aesthetic manner. We have further shown that, it is also the feebleness that renders the accumulation of lepton-number asymmetry to convert into the baryon-number asymmetry via rapid sphaleron transitions in the early Universe.

Conclusion {#sec:con}
==========

We have demonstrated in this paper that, when both thermal effects at high temperature and non-diagonal textures of lepton (both charged-lepton and neutrino) Yukawa matrices are considered, it is feasible to account for the matter-antimatter asymmetry of the Universe within a minimal freeze-in DN setup. While the SM Higgs cannot generate the observed baryon-number asymmetry in such a minimal setup, the second Higgs doublet of the neutrinophilic 2HDM, when being in equilibrium with the thermal bath, can realize the freeze-in DN. To establish a direct connection between the high-temperature leptonic CP asymmetry and the low-energy neutrino oscillation observables, we have considered various minimal corrections to the TB mixing pattern, and found that the patterns with a small mixing angle and a maximal CP-violating phase can produce compatible neutrino oscillation observables with a (negative) maximal CP measure and, at the same time, account for the matter-antimatter asymmetry of the Universe observed today.

Such a minimal setup realized in this paper is predictable on account of the correlation between the BAU and the neutrino oscillation observables, and might also be testable at colliders in terms of the electroweak scalars introduced to generate the neutrino masses and to implement the BAU.

This work is supported by the National Natural Science Foundation of China under Grant Nos. 11675061 and 11775092, as well as by the Fundamental Research Funds for the Central Universities under Grant Nos. CCNU18TS029 and 2019YBZZ079.
