# resonant leptogenesis with two sterile neutrinos. Equations from 0705.2183
import ulysses
import numpy as np
from odeintw import odeintw
from numba import jit
@jit
def fast_RHS(y0, eps2tt,eps2mm,eps2ee,eps1tt,eps1mm,eps1ee,C,d1,d2,w1,w2,n1eq,n2eq):
N1, N2, Ntt, Nmm, Nee = y0
c1t,c1m,c1e,c2t,c2m,c2e = C
c1tc = np.conjugate(c1t)
c1mc = np.conjugate(c1m)
c1ec = np.conjugate(c1e)
c2tc = np.conjugate(c2t)
c2mc = np.conjugate(c2m)
c2ec = np.conjugate(c2e)
#define the different RHSs for each equation
rhs1 = -d1*(N1-n1eq)
rhs2 = -d2*(N2-n2eq)
rhs3 = ( eps1tt * d1 * (N1-n1eq) + eps2tt * d2 * (N2-n2eq)
- w1 * ( c1t * c1tc * Ntt) - w2 * ( c2t * c2tc * Ntt))
rhs4 = ( eps1mm * d1 * (N1-n1eq) + eps2mm * d2 * (N2-n2eq)
- w1 * ( c1m * c1mc * Nmm) - w2 * ( c2m * c2mc * Nmm))
rhs5 = ( eps1ee * d1 * (N1-n1eq) + eps2ee * d2 * (N2-n2eq)
- w1 * ( c1e * c1ec * Nee) - w2 * ( c2e * c2ec * Nee))
RHStemp = [rhs1, rhs2, rhs3, rhs4, rhs5]
return RHStemp
[docs]class EtaB_2Resonant(ulysses.ULSBase):
"""
Resonant equations with two steriles and three lepton flavours. See arxiv:0705.2183.
"""
[docs] def RHS(self, y0, zzz, ETA, C, K):
k1term,k2term = K
eps2tt,eps2mm,eps2ee,eps1tt,eps1mm,eps1ee = ETA
if zzz != self._currz or zzz == self.zmin:
self._d1 = np.real(self.D1(k1term, zzz))
self._d2 = np.real(self.D2(k2term, zzz))
self._w1 = np.real(self.W1(k1term, zzz))
self._w2 = np.real(self.W2(k2term, zzz))
self._n1eq = self.N1Eq(zzz)
self._n2eq = self.N2Eq(zzz)
self._currz=zzz
return fast_RHS(y0, eps2tt, eps2mm, eps2ee, eps1tt, eps1mm, eps1ee, C,self._d1,self._d2,self._w1,self._w2,self._n1eq,self._n2eq)
@property
def EtaB(self):
_HT = [
np.real(self.hterm(2,0)),
np.real(self.hterm(1,0)),
np.real(self.hterm(0,0)),
np.real(self.hterm(2,1)),
np.real(self.hterm(1,1)),
np.real(self.hterm(0,1))
]
_K = [np.real(self.k1), np.real(self.k2)]
y0 = np.array([0+0j,0+0j,0+0j,0+0j,0+0j], dtype=np.complex128)
_ETA = [
self.epsiloniaaRES(2,1,0),
self.epsiloniaaRES(1,1,0),
self.epsiloniaaRES(0,1,0),
self.epsiloniaaRES(2,0,1),
self.epsiloniaaRES(1,0,1),
self.epsiloniaaRES(0,0,1)
]
_C = [ self.c1a(2), self.c1a(1), self.c1a(0), self.c2a(2), self.c2a(1), self.c2a(0)]
_K = [np.real(self.k1), np.real(self.k2)]
ys = odeintw(self.RHS, y0, self.zs, args = tuple([_ETA, _C, _K]), atol=1e-10)
nb = self.sphalfact*(ys[-1,2]+ys[-1,3]+ys[-1,4])
self.ys = np.real(ys[:, [2,3,4]])
return np.real(nb)