Implementation of Casas Ibarra parametrization of arXiv:1504.07892

t13a.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## T13A model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See [arXiv:1504.07892](http://arxiv.org/abs/1504.07892)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#%%writefile Modules/LFV.py\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "def f(m1,m2):\n",
    "    return (m1**2*np.log(m1**2)-m2**2*np.log(m2**2))/(m1**2-m2**2)\n",
    "\n",
    "def M_ND(ip):\n",
    "    Mp=ip['Mp']\n",
    "    mu=ip['mu']\n",
    "    beta=ip['beta']\n",
    "    y=ip['y']\n",
    "    mz=y*ip['vev']/np.sqrt(2.)\n",
    "    return np.array([[Mp        ,-mz*np.cos(beta),mz*np.sin(beta)],\n",
    "             [-mz*np.cos(beta),   0            ,-mu ],\n",
    "             [mz*np.sin(beta) ,  -mu          ,0 ]])\n",
    "\n",
    "def eigensystem_numerical(ip):\n",
    "    MND=np.matrix(M_ND(ip))\n",
    "\n",
    "    (mchi,N)=np.linalg.eig(MND)\n",
    "    return MND,mchi,N\n",
    "\n",
    "def Lambda(ip,mS):\n",
    "    import decimal as dc\n",
    "    MND,mchi,N=eigensystem_numerical(ip)\n",
    "    #Be carefull with cancellations! use dc.Decimal\n",
    "    return float(  dc.Decimal( N[2,0]**2*mchi[0]*f(mS,mchi[0]) )\\\n",
    "                  +dc.Decimal( N[2,1]**2*mchi[1]*f(mS,mchi[1]) )\\\n",
    "                  +dc.Decimal( N[2,2]**2*mchi[2]*f(mS,mchi[2]) ) )/(16.*np.pi**2)\n",
    "\n",
    "# Function in order to compute the Ma Expression matrix    \n",
    "def neutrino_data():\n",
    "    '''From arxiv:1405.7540 (table I)\n",
    "    and asumming a Normal Hierarchy:\n",
    "    Output:\n",
    "    mnu1in: laightest neutrino mass\n",
    "    Dms2: \\Delta m^2_{12}\n",
    "    Dma2: \\Delta m^2_{13}\n",
    "    ThetSol,ThetAtm,ThetRec: in radians\n",
    "    '''\n",
    "    Dms2=np.array([7.11e-5, 7.60e-5, 8.18e-5])*1e-18 # In GeV\n",
    "    Dma2=np.array([2.30e-3, 2.48e-3, 2.65e-3])*1e-18 # In GeV\n",
    "    #input real values:\n",
    "    #\n",
    "    ThetSol = np.array([0.278,  0.323,  0.375]) \n",
    "    ThetAtm = np.array([0.392,  0.567,  0.643])\n",
    "    ThetRec = np.array([0.0177, 0.0234, 0.0294])\n",
    "    mnu1in=1E-5*1E-9\n",
    "\n",
    "    return mnu1in,Dms2,Dma2,ThetSol,ThetAtm,ThetRec\n",
    "\n",
    "\n",
    "def CasasIbarra(ip,ranMnu=False,bestfit=False):\n",
    "    #import numpy as np\n",
    "    if ranMnu==True:\n",
    "        bestfit=True\n",
    "        \n",
    "    mnu1in,Dms2,Dma2,ThetSol,ThetAtm,ThetRec=neutrino_data()\n",
    "    \n",
    "    M1 = 1./( Lambda(ip,ip['m_S1']) + 0*1j)\n",
    "    M2 = 1./(Lambda(ip,ip['m_S2']) + 0*1j) \n",
    "    M3 = 1./(Lambda(ip,ip['m_S3'])+ 0*1j)\n",
    "    DMR = np.diag([np.sqrt(M1),np.sqrt(M2),np.sqrt(M3)])\n",
    "    \n",
    "    # Phases of the PMNS matrix\n",
    "    phases1=np.random.uniform(0.,0.0*np.pi,3)       # cero value for all phases         \n",
    "    \n",
    "    delta=1.*(0 if ranMnu else phases1[0])\n",
    "    eta1 =1.*(0 if ranMnu else phases1[1])\n",
    "    eta2 =1.*(0 if ranMnu else phases1[2])\n",
    "\n",
    "    mnu1=1.*(mnu1in      if bestfit else  10**((np.log10(2.5e-3)-np.log10(1e-9))*np.random.uniform(0,1)+np.log10(1e-9))*1e-9) \n",
    "    mnu2=1.*(np.sqrt(Dms2[1]+mnu1in**2)     if bestfit else  np.sqrt(np.random.uniform(Dms2[0],Dms2[2]) + mnu1**2)  ) \n",
    "    mnu3=1.*(np.sqrt(Dma2[1]+mnu1in**2)     if bestfit else  np.sqrt(np.random.uniform(Dma2[0],Dma2[2]) + mnu1**2)  )\n",
    "    \n",
    "    # Square root of left-handed neutrino mass matrix \n",
    "    DMnu = np.diag( [np.sqrt(mnu1),np.sqrt(mnu2),np.sqrt(mnu3)] )  \n",
    "\n",
    "    #print \"NEUTRINOS\",mnu1,mnu2,mnu3\n",
    "\n",
    "    # Mixing angles (up 3 sigma range)\n",
    "    t12 = 1.*( np.arcsin(np.sqrt(ThetSol[1])) if bestfit else np.arcsin(np.sqrt(np.random.uniform(ThetSol[0],ThetSol[2]))))\n",
    "    t23 = 1.*( np.arcsin(np.sqrt(ThetAtm[1])) if bestfit else np.arcsin(np.sqrt(np.random.uniform(ThetAtm[0],ThetAtm[2]))))\n",
    "    t13 = 1.*( np.arcsin(np.sqrt(ThetRec[1])) if bestfit else np.arcsin(np.sqrt(np.random.uniform(ThetRec[0],ThetRec[2]))))         \n",
    "\n",
    "    # Building PMNS matrix\n",
    "    U12 = np.array([ [np.cos(t12), np.sin(t12),0], [-np.sin(t12), np.cos(t12),0], [0,0,1.0] ])\n",
    "    U13 = np.array([ [np.cos(t13),0, np.sin(t13)* np.exp(-delta*1j)], [0,1.0,0], [-np.sin(t13)*np.exp(delta*1j),0, np.cos(t13)] ])\n",
    "    U23 = np.array([ [1.0,0,0], [0, np.cos(t23), np.sin(t23)], [0, -np.sin(t23), np.cos(t23)] ])\n",
    "    Uphases = np.array([ [np.exp(eta1*1j),0,0], [0, np.exp(eta2*1j),0], [0,0,1.0] ])\n",
    "    U = np.dot(U23,np.dot(U13,np.dot(U12,Uphases)))\n",
    "\n",
    "    # Building R matrix of the Casas-Ibarra parametrization\n",
    "\n",
    "    # Real angles.\n",
    "    thetas = np.random.uniform(0.0,2.0*np.pi,3)#real part of mixing angles of R\n",
    "    b12 = 1.*(0 if ranMnu else thetas[0])\n",
    "    b13 = 1.*(0 if ranMnu else thetas[1])\n",
    "    b23 = 1.*(0 if ranMnu else thetas[2])\n",
    "    \n",
    "    #print \"angulos\", b12, b13, b23    \n",
    "\n",
    "    R12 = np.array([ [np.cos(b12),np.sin(b12),0], [-np.sin(b12),np.cos(b12),0], [0,0,1.0] ])\n",
    "    R13 = np.array([ [np.cos(b13),0,np.sin(b13)], [0,1.0,0], [-np.sin(b13),0,np.cos(b13)] ])\n",
    "    R23 = np.array([ [1.0,0,0], [0,np.cos(b23),np.sin(b23)], [0,-np.sin(b23),np.cos(b23)] ])\n",
    "\n",
    "    R=np.dot(R23,np.dot(R13,R12))  \n",
    "\n",
    "    # Yukawa matrix of the Casas-Ibarra parametrization\n",
    "    \n",
    "    yuk2 = np.dot( DMR, np.dot( R,np.dot( DMnu,U.transpose().conjugate() ) ) )\n",
    "\n",
    "    # Yukawa matrix in the T13A model\n",
    "    yuk = np.transpose(yuk2)\n",
    "\n",
    "    return yuk\n",
    "\n",
    "def test_CasasIbarra():\n",
    "    vev=246.\n",
    "    ip=pd.Series({'Mp':100.,'mu':50.,'y':1.*np.sqrt(2.)/vev,'beta':np.pi/6.,'vev':vev,\\\n",
    "                  'm_S1':80.,'m_S2':800.,'m_S3':1500.})\n",
    "\n",
    "    mnu1in,Dms2,Dma2,ThetSol,ThetAtm,ThetRec=neutrino_data()\n",
    "    Mnuin=np.array([mnu1in,np.sqrt(Dms2[1]+mnu1in**2),np.sqrt(Dma2[1]+mnu1in**2)])\n",
    "    h=CasasIbarra(ip,ranMnu=True)\n",
    "                  \n",
    "    #Diagonal matrix\n",
    "    Lamb = np.diag( [Lambda(ip,ip['m_S1']),Lambda(ip,ip['m_S2']),Lambda(ip,ip['m_S3'])] )\n",
    "\n",
    "    #              3x3        3x3   3x3\n",
    "    Mint = np.dot( h,np.dot(Lamb,h.transpose()) )\n",
    "    Mnu,U=np.linalg.eig(Mint) \n",
    "    lo=np.argsort(np.abs(Mnu))\n",
    "    Mnu=np.array([Mnu[lo[0]],Mnu[lo[1]],Mnu[lo[2]]])\n",
    "    U=np.matrix(U)\n",
    "    U=np.asarray(np.hstack((U[:,lo[0]],U[:,lo[1]],U[:,lo[2]])))\n",
    "    Upmns=np.array([[ 0.81311635+0.j,  0.56164206+0.j,  0.15297059+0.j],\\\n",
    "                    [-0.46875227+0.j,  0.47596125+0.j,  0.74413184+0.j],\\\n",
    "                    [ 0.34512767+0.j, -0.67677108+0.j,  0.65028286+0.j]])\n",
    "\n",
    "    #print(np.abs(Mnuin)*1e9,np.abs(Mnu)*1e9)\n",
    "    #print(np.abs(Upmns),np.abs(U))\n",
    "    np.testing.assert_array_almost_equal(np.abs(Mnuin)*1e9,np.abs(Mnu)*1e9)\n",
    "    np.testing.assert_array_almost_equal(np.abs(Upmns),np.abs(U))\n",
    "###Punto problematico del profe\n",
    "    \n",
    "# Analytical \n",
    "def Fme(x,xmin=0.996,xmax=1.005,xfit=1.001):\n",
    "    \"\"\"Fixing near to one values\n",
    "     xmin: close to 1 from below\n",
    "     xmax: close to 1 from above\n",
    "     xfit:  optimized 1 limit\n",
    "    \"\"\"\n",
    "    x=np.asarray(x)\n",
    "    if x.shape:\n",
    "        x[np.logical_and(x>xmin,x<xmax)]=xfit\n",
    "    else:\n",
    "        if x>xmin and x<xmax:\n",
    "            x=xfit\n",
    "        \n",
    "    return (2.+3.*x-6.*x**2+x**3+6*x*np.log(x))/(6.*(1.-x)**4)\n",
    "\n",
    "def muegamma(ip,h):\n",
    "    #import pandas as pd\n",
    "    #import numpy as np\n",
    "    ip=pd.Series(ip)\n",
    "    h=np.matrix(h)\n",
    "    alpha=1./137.\n",
    "    GF=1.166364E-5\n",
    "    FM=np.matrix(np.diag( [Fme(ip.mu**2/ip.m_S1**2)/ip.m_S1**2,\\\n",
    "               Fme(ip.mu**2/ip.m_S2**2)/ip.m_S2**2,\\\n",
    "               Fme(ip.mu**2/ip.m_S3**2)/ip.m_S3**2 ] ))\n",
    "    \n",
    "    return (3*alpha/(4.*16.*np.pi*GF**2)*np.abs(h[0,:]*FM*h[1,:].H)**2)[0,0]\n",
    "\n",
    "if __name__=='__main__':\n",
    "    test_CasasIbarra()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check Casas-Ibarra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the Yukawa couplings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vev=246.\n",
    "ip=pd.Series({'Mp':100.,'mu':50.,'y':1.*np.sqrt(2.)/vev,'beta':np.pi/6.,'vev':vev,\\\n",
    "                  'm_S1':80.,'m_S2':800.,'m_S3':1500.})\n",
    "h=CasasIbarra(ip)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use $h$ to build the neutrino mass matrix and get --in particular--  the mixing matrix $U$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Diagonal matrix\n",
    "Lamb = np.diag( [Lambda(ip,ip['m_S1']),Lambda(ip,ip['m_S2']),Lambda(ip,ip['m_S3'])] )\n",
    "Mint = np.dot( h,np.dot(Lamb,h.transpose()) )\n",
    "Mnu,U=np.linalg.eig(Mint) \n",
    "lo=np.argsort(np.abs(Mnu))\n",
    "Mnu=np.array([Mnu[lo[0]],Mnu[lo[1]],Mnu[lo[2]]])\n",
    "U=np.matrix(U)\n",
    "U=np.asarray(np.hstack((U[:,lo[0]],U[:,lo[1]],U[:,lo[2]])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mixing matrix with the usual PDG conventions is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.80565357+0.j, -0.56931478+0.j,  0.16371626+0.j],\n",
       "       [-0.46229000+0.j, -0.43142069+0.j,  0.77470261+0.j],\n",
       "       [ 0.37041906-0.j,  0.69982632+0.j,  0.61076415+0.j]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which is within the 3$\\sigma$ ranges:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|UPMNS_min|:\n",
      "[[ 0.77886135  0.51944855  0.13304135]\n",
      " [ 0.40568185  0.38816445  0.61682672]\n",
      " [ 0.21651102  0.55850112  0.58864607]]\n",
      "|UPMNS_max|:\n",
      "[[ 0.84215236  0.60692874  0.17146428]\n",
      " [ 0.56236395  0.61863368  0.79474455]\n",
      " [ 0.42820125  0.73537296  0.77281201]]\n"
     ]
    }
   ],
   "source": [
    "UP=[]\n",
    "UPmin=np.zeros((3,3))\n",
    "UPmax=np.zeros((3,3))\n",
    "mnu1in,Dms2,Dma2,ThetSol,ThetAtm,ThetRec=neutrino_data()\n",
    "for i in [0,2]:\n",
    "    for j in [0,2]:\n",
    "        for k in [0,2]:\n",
    "            t12 = np.arcsin(np.sqrt(ThetSol[i]))\n",
    "            t23 = np.arcsin(np.sqrt(ThetAtm[j]))\n",
    "            t13 = np.arcsin(np.sqrt(ThetRec[k]))\n",
    "            R12 = np.array([ [np.cos(t12), np.sin(t12),0], [-np.sin(t12), np.cos(t12),0], [0,0,1.0] ])\n",
    "            R13 = np.array([ [np.cos(t13),0, np.sin(t13)], [0,1.0,0], [-np.sin(t13),0, np.cos(t13)] ])\n",
    "            R23 = np.array([ [1.0,0,0], [0, np.cos(t23), np.sin(t23)], [0, -np.sin(t23), np.cos(t23)] ])\n",
    "            UP.append(np.dot(R23,np.dot(R13,R12)))\n",
    "\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        UPmin[i,j]=np.abs(np.array([ UP[l][i,j] for l in range(8)])).min()\n",
    "        UPmax[i,j]=np.abs(np.array([ UP[l][i,j] for l in range(8)])).max()\n",
    "        \n",
    "print '|UPMNS_min|:'\n",
    "print UPmin\n",
    "print '|UPMNS_max|:'\n",
    "print UPmax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "of the experimental\n",
    "\\begin{align}\n",
    "U_{\\text{PMNS}}=\\begin{pmatrix}\n",
    " 0.81311635  &  0.56164206 & 0.15297059\\\\\n",
    " -0.46875227 &  0.47596125 & 0.74413184\\\\\n",
    " 0.34512767  & -0.67677108 &  0.65028286\\\\\n",
    "\\end{pmatrix}              \n",
    "\\end{align}"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}