{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ultrarelativistic right-handed neutrino production, reproducing the $2\\leftrightarrow 2$ part of [1202.1288](https://arxiv.org/pdf/1202.1288)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importing the modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#start by importing the controller and manipulate modules\n",
    "from analytical.controller import *\n",
    "from numerical.manipulate import *\n",
    "#reimport numpy (though it is pulled by numerical) for smarter syntax highlighting in vscode\n",
    "import numpy\n",
    "#import matplotlib too\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analytical part"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The relevant part of the config file is as follows\n",
    "```ini\n",
    "[Model]\n",
    "modelpath = /Users/jacopo/NextCloud/AUTOTHERM/autotherm/analytical/models/righthandedneutrino.fr\n",
    "# Symbol for the Lagrangian in the model file\n",
    "lagrangian = Ltot\n",
    "# \"Name\" of the particle whose production rate must be computed\n",
    "produced = F[6]\n",
    "# List of the particles in the thermal bath (or leave empty for SM assumption)\n",
    "inbath = \n",
    "assumptions = Element[hhh,Reals]\n",
    "replacements =\n",
    "includeSM = yes\n",
    "noneq = hhh\n",
    "#comma-separated list of particles to be treated with a flavor expansion\n",
    "flavorexpand = \n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rhn_analytical=analytical_pipeline(\"../../MyModels/RHN/RHN.cfg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inspect the output. The `hhh` coupling in the model file stands for $\\sqrt{\\mathrm{Tr}h^\\dagger h}$, with $h$ the Yukawa matrix between the right-handed neutrinos and the left-handed lepton doublet. In this models there are three generations of massless right-handed neutrinos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "- statistics: (-1, -1, -1):  $$36 \\mathrm{Tr}[h^\\dagger h] \\vert h_t\\vert^2$$\n",
       "- statistics: (1, -1, 1):  $$- \\frac{\\mathrm{Tr}[h^\\dagger h] \\left(g_{1}^{2} + 3 g_{2}^{2}\\right) \\left(s^{2} + t^{2}\\right)}{s t}$$\n",
       "- statistics: (1, 1, -1):  $$- \\frac{\\mathrm{Tr}[h^\\dagger h] \\left(g_{1}^{2} + 3 g_{2}^{2}\\right) \\left(s^{2} + u^{2}\\right)}{s u}$$\n",
       "- statistics: (-1, 1, 1):  $$\\frac{\\mathrm{Tr}[h^\\dagger h] \\left(g_{1}^{2} + 3 g_{2}^{2}\\right) \\left(t^{2} + u^{2}\\right)}{t u}$$\n"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import display, Markdown\n",
    "out=\"\"\n",
    "s = sympy.Symbol(\"s\")\n",
    "t = sympy.Symbol(\"t\")\n",
    "u = sympy.Symbol(\"u\")\n",
    "g1 = sympy.Symbol(\"g1\")\n",
    "g2 = sympy.Symbol(\"g2\")\n",
    "g3 = sympy.Symbol(\"g3\")\n",
    "ht = sympy.Symbol(\"ht\")\n",
    "hhh=sympy.Symbol(\"hhh\")\n",
    "for key, item, in rhn_analytical[0].items():\n",
    "    sitem=sympy.simplify(sympy.sympify(item).subs(2*t*u,(s*s-t*t-u*u)))\n",
    "    sitemsubst= sitem\n",
    "    if sitemsubst.as_coeff_Add()[0] ==0:\n",
    "        sitemsubstcollect = sitemsubst\n",
    "    else:\n",
    "        sitemsubstcollect=0\n",
    "        for exprtemp in sympy.factor(sitemsubst.expand().as_independent(ht)).subs(t+u,-s)\\\n",
    "            .subs(2*t*t+2*t*u,2*t*t+(s*s-t*t-u*u)).subs(t*t+2*t*u,t*t+(s*s-t*t-u*u)).subs(t*t+t*u,t*t+(s*s-t*t-u*u)/2):\n",
    "            sitemsubstcollect += exprtemp.simplify()\n",
    "        # beautify the output\n",
    "    out+=f\"- statistics: {key}:  $${sympy.latex(sitemsubstcollect).replace('ht^{2}',r'\\vert h_t\\vert^2').\\\n",
    "                                    replace('hhh^{2}',r'\\mathrm{Tr}[h^\\dagger h]')}$$\\n\"\n",
    "    # display(sympy.pprint(sitem))\n",
    "display(Markdown(out))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical part"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First, we must call NumRate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set the degeneracy to unity, since we are treating the right-handed neutrinos as massless and thus distinct from their antiparticles. The rate for the latter is then equal to that for the former that is determined here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "rate_rhn=NumRate(*rhn_analytical,1) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The leading-log term reads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{T^{2} hhh^{2} \\left(g_{1}^{2} + 3 g_{2}^{2}\\right) \\log{\\left(\\frac{64 k^{2}}{T^{2} \\left(g_{1}^{2} + 3 g_{2}^{2}\\right)} \\right)}}{256 \\pi k \\tanh{\\left(\\frac{k}{2 T} \\right)}}$"
      ],
      "text/plain": [
       "T**2*hhh**2*(g1**2 + 3*g2**2)*log(64*k**2/(T**2*(g1**2 + 3*g2**2)))/(256*pi*k*tanh(k/(2*T)))"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "leadlog=rate_rhn.get_leadlog().simplify()\n",
    "leadlog"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "reexpress this in terms of the lepton asymptotic mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\mathrm{Tr}[h^\\dagger h] m_{\\infty\\,L}^2 \\left(2 n_{B}{\\left(k \\right)} + 1\\right) \\log{\\left(\\frac{4 k^{2}}{m_{\\infty\\,L}^2} \\right)}}{16 \\pi k}$"
      ],
      "text/plain": [
       "\\mathrm{Tr}[h^\\dagger h]*m_{\\infty\\,L}^2*(2*n_B(k) + 1)*log(4*k**2/m_{\\infty\\,L}^2)/(16*pi*k)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "minf = sympy.Symbol(r\"m_{\\infty\\,L}^2\")\n",
    "htrace=sympy.Symbol(r'\\mathrm{Tr}[h^\\dagger h]')\n",
    "T = sympy.S(\"T\")\n",
    "k= sympy.S(\"k\")\n",
    "nB= sympy.Function(\"n_B\")\n",
    "leadlog.subs(g2*g2,16*(minf/T/T-g1*g1/16)/3).subs(sympy.tanh(k/(2*T)),1/(1+2* nB(k))).subs(hhh*hhh,htrace).simplify()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This agrees with Eqs. (24) and (28) of of [1202.1288](https://arxiv.org/pdf/1202.1288) and with (4.23) of [1605.07720](https://arxiv.org/abs/1605.07720)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now prepare to compare with the integration results in (29) of of [1202.1288](https://arxiv.org/pdf/1202.1288)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#prepare to compute the total rate in the strict LO scheme\n",
    "def ratevec(x,g1,g2,ht):\n",
    "    #add factor of 2 for the \"antiparticle\"\n",
    "    #factor out factor of 2/(2\\pi^2)\n",
    "    statfac=x*x/(numpy.exp(x)+1.)\n",
    "    evalf=rate_rhn.rate(x,1,tuple([g1,g2,ht]),1)[1][0]\n",
    "    return statfac*evalf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we first test $c_Q$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "### We find $c_Q$=2.52136"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.integrate import quad\n",
    "#re-instate factor of 2/(2\\pi^2)\n",
    "yukawaint=1/numpy.pi**2*quad(ratevec,0.,20.,args=(1e-7,1e-7,1),epsrel=1e-5)[0]\n",
    "display(Markdown(r'### We find $c_Q$=%.5f'%(yukawaint*1536*numpy.pi)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It checks out! Let us then test $c_V$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "### We find $c_V$=3.16694"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "g1test=1/numpy.sqrt(2)\n",
    "g2test=1/numpy.sqrt(6)\n",
    "#re-instate factor of 2/(2\\pi^2)\n",
    "gaugeint=1/numpy.pi**2*quad(ratevec,0.,20.,args=(g1test,g2test,0.),epsrel=1e-5)[0]\n",
    "display(Markdown(r'### We find $c_V$=%.5f'%(gaugeint*1536*numpy.pi/(g1test**2+3*g2test**2)+numpy.log(g1test**2+3*g2test**2))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "### We find $c_V$=3.16693"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "g1test=1\n",
    "g2test=1\n",
    "gaugeint=1/numpy.pi**2*quad(ratevec,0.,20.,args=(g1test,g2test,0.),epsrel=1e-5)[0]\n",
    "display(Markdown(r'### We find $c_V$=%.5f'%(gaugeint*1536*numpy.pi/(g1test**2+3*g2test**2)+numpy.log(g1test**2+3*g2test**2))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Checks out as well!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}