WIMP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/restrepo/WIMP/blob/master/WIMP.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "qFOeVPf9k7aL"
   },
   "source": [
    "# Dark Matter Boltzmann equation\n",
    "\n",
    "This program is made to reproduce the behavior of dark matter yield in WIMP and FIMP frameworks based on Chapter 5th, Kolb Turner (Early Universe) and https://arxiv.org/pdf/1009.3690.pdf\n",
    "\n",
    "See also\n",
    "* https://arxiv.org/pdf/1703.07364.pdf\n",
    "* [DM book](https://drive.google.com/file/d/1YVMIrldAAWWDYuMEDc6dKLS4g-13bC1d/view?usp=sharing)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 35
    },
    "colab_type": "code",
    "id": "MfyBj2Vqk7aS",
    "outputId": "aa773d2a-5642-42e1-c902-08ce97bab9b1"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "%pylab is deprecated, use %matplotlib inline and import the required libraries.\n",
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline  \n",
    "import numpy as np\n",
    "from numpy import arange\n",
    "from scipy.integrate import odeint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Expected cross sections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From: [1] https://arxiv.org/pdf/1009.3690.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\Omega h^{2} \\approx \\frac{3 \\times 10^{-27} \\mathrm{~cm}^{3} / \\mathrm{s}}{\\left\\langle\\sigma_{\\mathrm{ann}} v\\right\\rangle}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* [$1\\text{ cm}=5.068\\times 10^{13} \\hbar c\\cdot\\text{GeV}^{-1} $](https://www.wolframalpha.com/input?i=1cm+to+hbar*c/GeV)\n",
    "* [$1\\text{ s}=1.519\\times 10^{24} \\hbar \\cdot\\text{GeV}^{-1} $](https://www.wolframalpha.com/input?i=1s+to+hbar/GeV)\n",
    "* [$0.1\\text{ pb}\\approx 2.6\\times 10^{-10}\\, \\hbar c\\cdot \\text{GeV}^{-2}$](https://www.wolframalpha.com/input?i=2.6*10^{-10}GeV^{-2}(hbar*c)**2+to+picobarn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\langle\\sigma v \\rangle= 3\\times 10^{-27}\\frac{(5.068\\times 10^{13} \\hbar c\\cdot\\text{GeV}^{1})^3}{1.519\\times 10^{24}\\hbar\\cdot\\text{GeV}^{-1}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 923,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.5708296464516136e-10"
      ]
     },
     "execution_count": 923,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "3E-27*(5.068E13)**3/1.519E24"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\langle\\sigma v \\rangle=2.6\\times 10^{-10}(\\hbar c)^2\\text{GeV}^{-2}c$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 924,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.8"
      ]
     },
     "execution_count": 924,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2.6E-10*3E10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\langle\\sigma v \\rangle=7.8\\times (\\hbar c)^2\\text{GeV}^{-2}\\frac{\\text{cm}}{\\text{s}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Electroweak cross section are of order of pb\n",
    "$$\n",
    "\\sigma\\sim \\frac{\\langle\\sigma v\\rangle}{c}=2.6\\times 10^{-10}(\\hbar c)^2\\text{GeV}^{-2}\n",
    "=0.1\\text{ pb}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Expected abundances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://pdg.lbl.gov/2012/reviews/rpp2012-rev-astrophysical-constants.pdf\n",
    "* $h=0.71$\n",
    "* $\\rho_c^0=1.05375(13) \\times 10^{-5} h^{2}\\left(\\mathrm{GeV} / c^{2}\\right) \\mathrm{cm}^{-3} $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From eq. (7.11) https://arxiv.org/pdf/1009.3690.pdf\n",
    "$$ \\Omega_{\\chi} h^{2}=\\frac{\\rho_{\\chi}^{0} h^{2}}{\\rho_{c}^{0}}=\\frac{m_{\\chi} s_{0} Y_{0} h^{2}}{\\rho_{c}^{0}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From https://arxiv.org/pdf/hep-ph/0501232.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "s_0=\\left(2 \\pi^{2} / 45\\right) h_{\\mathrm{eff}}\\left(T_{\\gamma}\\right) T_{\\gamma}^{3}\n",
    "$$\n",
    "* $h_{\\mathrm{eff}}(T_\\gamma)=3.9172$\n",
    "* $T_\\gamma=2.725 K$[$=2.4\\times 10^{-13}\\text{ GeV}/k_B$](https://www.wolframalpha.com/input?i=2.725+K+to+GeV/k)$=2.4\\times 10^{-13}\\cdot 5.068\\times 10^{13}\\hbar c/(k_B\\text{cm})=\n",
    "12.2\\,\\hbar c/(k_B\\text{cm})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 927,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3120.128389914731"
      ]
     },
     "execution_count": 927,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(2*np.pi**2/45)*3.9172*(12.2)**3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$s_0\\approx 3120\\text{ cm}^3$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\frac{s_0 h^2}{\\rho_c^0}=\\frac{3120\\cdot 0.71^2}{1.05375\\times 10^{-5}}\\approx 1.5\\times 10^8$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 992,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "149256654.80427045"
      ]
     },
     "execution_count": 992,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "3120*0.71**2/1.05375E-5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So that"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ \\Omega_{\\chi} h^{2}=1.5\\times 10^8  m_{\\chi}  Y_{0} (c^2/ \\text{GeV})=1.5\\times 10^{10} Y_0 \\left(\\frac{m_\\chi}{100\\text{ GeV/}c^2}\\right) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To obtain\n",
    "$$ \\Omega_{\\chi} h^{2}=0.118$$\n",
    "we require the adimensional quantity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 943,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.866666666666666e-12"
      ]
     },
     "execution_count": 943,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y0=0.118/1.5E10\n",
    "Y0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$Y_0=8\\times 10^{-12}\\left(\\frac{100\\text{ GeV}/c^2}{m_\\chi}\\right)\\sim 10^{-11}\\left(\\frac{100\\text{ GeV}/c^2}{m_\\chi}\\right) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Expected result\n",
    "From [1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image](./Y.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use eq. (7.10) of https://arxiv.org/pdf/1009.3690.pdf\n",
    "$$\n",
    "\\frac{d Y}{d x}=-\\left(\\frac{45}{\\pi M_{P}^{2}}\\right)^{-1 / 2} \\frac{g_{*}^{1 / 2} m}{x^{2}}\\langle\\sigma v\\rangle\\left(Y^{2}-Y_{\\text{eq}}^{2}\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The equilibrium distribution of this particles is consider for the non-relativistic case, as follows (see Kolb & Turner eq 5.25):\n",
    "\n",
    "$$Y_{\\text{eq}}(x) = \\frac{45}{2\\pi^{4}}\\frac{g}{g_{*s}} x^{3/2} e^{-x} = 0.145\\frac{g}{g_{*s}} x^{3/2} e^{-x}\\,, $$\n",
    "\n",
    "where $x = M / T$ and $M= 100$ GeV is the singlet mass taken as constant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 591,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "kfDBXgSBk7aj"
   },
   "outputs": [],
   "source": [
    "# parameters\n",
    "M = 100                                #GeV Singlet Mass\n",
    "Mp = 1.22e19                            #GeV Planck Mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 592,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From fig. 1 https://arxiv.org/pdf/hep-ph/0501232.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 944,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfg=pd.DataFrame([\n",
    "    {'T':1E-1,'geff_to_half':4.2},\n",
    "    {'T':1,'geff_to_half':8.7},\n",
    "    {'T':100,'geff_to_half':10.2}\n",
    "    ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 945,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import optimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 946,
   "metadata": {},
   "outputs": [],
   "source": [
    "f=lambda T,a,b,c:a/(1+np.exp(-b*(T-c)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 974,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.9/dist-packages/scipy/optimize/minpack.py:833: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n"
     ]
    }
   ],
   "source": [
    "B,σ=optimize.curve_fit(f,dfg['T'],dfg['geff_to_half'],p0=(0,0,0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 975,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs_to_half= lambda T: f(T,a=B[0],b=B[1],c=B[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 976,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$g_*^{1/2}$')"
      ]
     },
     "execution_count": 976,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAESCAYAAAAIfCk9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhF0lEQVR4nO3dd5RV9b3+8fdnKgwgUgakzYwKEhEsOAiYYNfYicZYAraIxMQYzU31kmiSG025/mKKubkXSywgRomYYlewoEhH6UgbYGhDL0OZ8vn9cUYdJjNMn+8+5zyvtWadmb332fOstTk8s9t3m7sjIiJSnZTQAUREJLpUEiIiUiOVhIiI1EglISIiNVJJiIhIjVQSIiJSo7TQAZpS586dPS8vL3QMEZG4Mnv27C3unl3dvIQqiby8PGbNmhU6hohIXDGzgprm6XCTiIjUSCUhIiI1UkmIiEiNVBIiIlIjlYSIiNRIJSEiEs/Gj4e8PEhJib2OH9+kq0+oS2BFRJLK+PEwejQUF8d+LiiI/QwwYkST/ApLpOdJ5Ofnu+6TkCgoK3c27drPjuISdu0vYde+EnbtL2XXvhJ27y/lQGkZJWXllJQ5B8vKKSktp6SsnNJyxx3K/bPXcgd355NPauXvJclNeQv27QPg4qVTuXr+G7HpubmwenWdV2Nms909v7p52pMQaYSd+0r4cO0OVm/dy+otxbHXrXtZu62YkrKa/ytPTzXSU1M+/cpINdJSU0hLNVLMMIi92mevZhCb88n3kvQsHVqnA7AvvdVn09esabJfoZIQqYfycmfRhl28vayIt5ZuZs6aHZSVx8ogKyOV3E5t6Nu1HRf0O4qcjll0bJPOEa3SOaL1J69ptM1MIy1VpwOlCeSNjB1iqionp8l+hUpCpBbuztTlW3hh7nreXlbElj0HABjQoz3fOPNYTu/did7Zbclul4mZ/r6XFnTffYeekwDIyopNbyKRKAkzewy4FNjs7v0rpnUE/grkAauBq919e6iMknzcnfeWb+XBN5Yxu2A7R2alM6xPNmcdl80Zx2WT3S4zdERJdp+cnB4zJnaIKScnVhBNdNIaInLi2szOAPYAT1Yqid8A29z9V2b2I6CDu//wcOvRiWtpKu+v2MKDry9j5urtdGvfitvP7s3V+b3ISNNhIkk8kT9x7e7vmFlelcnDgbMqvn8CeAs4bEmINNbsgm385pWlTF+1ja5HZPLz4SdwzaBeZKalho4mEkQkSqIGXd19Q8X3G4GuIcNIYnN3Hnl3Fb98eTGd2mZy72X9uO60HFqlqxwkuUW5JD7l7m5m1R4XM7PRwGiAnCY8oy/JY39JGXc/P59Jcwu5qP9RPPCVk2iTGRcfDZFmF+VPwiYz6+buG8ysG7C5uoXcfSwwFmLnJFoyoMS/9Tv28fWnZjO/cCffPf84bj+7NykpukJJ5BNRLol/ADcCv6p4/XvYOJJoZq7exjfGzWZ/STkP35DP+f10RFOkqkiUhJlNIHaSurOZrQPuJVYOz5rZLUABcHW4hJJoJsxYwz1/X0DPDlk8M/pUendpFzqSSCRFoiTc/boaZp3bokEkKTw7cy13Pz+fM4/L5g/XnUL7imENROTfRaIkRFrK+yu28J+T5jOsT2ceuTGfdA2PIXJY+oRI0lhRtIdvjJvD0Z3b8NBXB6ogROpAnxJJCtv3HuRrj88kLcV47KZBOsQkUkc63CQJ70BpGV9/ajYbdu5nwq2D6dUxK3QkkbihPQlJaO7O3c/PZ8bqbfz3VSdyam7H0JFE4opKQhLan6Ys5/k5hdx1Xh+Gn9wjdByRuKOSkIQ1ZelmHnhtGcNP7s6d5/YJHUckLqkkJCHtPVDKjyctoHeXtvz6yyfqYUAiDaQT15KQHnhtKet37mPibUM1kqtII2hPQhLO3DXbefz91YwcnKsT1SKNpJKQhFJSVs7dz8+na7tW/ODCvqHjiMQ9HW6ShDL2nZUs2bibh2/Ip10r3TAn0ljak5CEsbJoD79/82MuHnCUhv0WaSIqCUkI5eWxm+ZapaXw08tPCB1HJGGoJCQhPDtrLdNXbeM/Lz6eLu1ahY4jkjBUEhL3Nu/az/0vLWbIMR25ZlCv0HFEEopKQuLeb15dyv7Scu6/YoBumhNpYioJiWvLN+/h+TnruGFILsdktw0dRyThqCQkrj34xjJap6fyjbOODR1FJCGpJCRuLVy/kxc/2sDXvnA0ndpmho4jkpBUEhK3fvvaMo5olcaoYceEjiKSsFQSEpfmrNnOm0s28/Uzj9WjSEWakUpC4tL/e20pndtmcNPpeaGjiCQ0lYTEnfdXbOG95Vv5xlm9aZOp4cdEmlOkS8LM7jSzBWa20MzuCp1HwnN3Hnh1Kd3at2LE4JzQcUQSXmRLwsz6A7cCpwEnAZeaWe+wqSS0KUs3M2fNDu44p48eJiTSAiJbEsDxwHR3L3b3UuBt4MrAmSSg8nLngVeXkdspi6/k9wwdRyQpRLkkFgDDzKyTmWUBFwMamCeJvbJwI4s27OKu8/qQnhrlf7oiiSOyZ/3cfbGZ/Rp4DdgLzAPKqi5nZqOB0QA5OTpGnajcnT+8+TF9urTl8pN6hI4jkjQi/eeYuz/q7qe6+xnAdmBZNcuMdfd8d8/Pzs5u+ZDSIt75eAtLNu7mtjOPJTVFg/iJtJTI7kkAmFkXd99sZjnEzkcMCZ1Jwnj4nZV0PSKTy07qHjqKSFKJdEkAfzOzTkAJcLu77wicRwJYuH4nU5dv4UcXfY6MtEjv/IoknEiXhLsPC51Bwnvk3VW0yUjlutN0zkmkpenPMom09Tv28c8P13PtaTkao0kkAJWERNP48ZCXx+NXfgsvLeXmonmhE4kkJZWERM/48TB6NLs2FPH0SRdyyZJ36fmtW2PTRaRFqSQkesaMgeJinjn5i+zJzOLWGZOguDg2XURalEpComfNGg6mpPFY/nCGFnzIgE0rPp0uIi1LJSHRk5PDi5/7AhvbdWb0jEmHTBeRlhXpS2AlOfkv7mPslB302VLAmStnxyZmZcF994UNJpKEtCchkfPe4C+yODuPW1e8S4oBubkwdiyMGBE6mkjS0Z6ERM7Yd1fSuW0mw18fB2kTQscRSWrak5BIWb55N+8sK+LGoblkpumhQiKhqSQkUp6aVkBGagrX6dGkIpGgkpDI2HOglL/NKeSSE7vRuW1m6DgigkpCImTSnHXsOVDKDUNzQ0cRkQoqCYkEd+eJaQWc2LM9J/c6MnQcEamgkpBImLZiK8s37+GGoXmY6clzIlGhkpBIeHJaAR2y0rn0xG6ho4hIJSoJCa5wxz5eW7SRawbl0Cpdl72KRIlKQoJ7enoBACN02atI5KgkJKgDpWU8M2Mt5x7flV4ds0LHEZEqVBIS1EvzN7B170FuHJoXOoqIVEMlIUE98X4Bx2S34fO9O4WOIiLVUElIMB+u3cG8tTu4YUiuLnsViSiVhATz5LQC2mSk8uVTe4aOIiI1UElIENv2HuSfH63nyoE9adcqPXQcEalB5EvCzL5jZgvNbIGZTTCzVqEzSeM9N2stB0vLuV7jNIlEWqRLwsx6AN8G8t29P5AKXBs2lTRWebkzfvoaTju6I8d1bRc6jogcRqRLokIa0NrM0oAsYH3gPNJI73xcxJptxVw/RHsRIlEX6ZJw90LgAWANsAHY6e6vhU0ljTXugwI6t83kiyccFTqKiNQi0iVhZh2A4cDRQHegjZmNrLLMaDObZWazioqKQsSUeli3vZjJSzZz7aBeZKRF+p+fiBDxkgDOA1a5e5G7lwDPA6dXXsDdx7p7vrvnZ2dnBwkpdTdhxhoAPZ5UJE5EvSTWAEPMLMtid1udCywOnEka6GBpOX+duZZzPteVHke2Dh1HROog0iXh7tOBicAcYD6xvGODhpIGe2XhRrbsOcjIIdqLEIkXaaED1Mbd7wXuDZ1DGm/cBwXkdMzijD46LCgSLyK9JyGJY9mm3cxYtY0Rg3NISdE4TSLxQiUhLWLcBwVkpKXwlfxeoaOISD2oJKTZ7T1QyvNzCrl0QDc6tskIHUdE6kElIc3uhXmF7DlQygjdYS0Sd1QS0qzcnaemFdCv2xEMzDkydBwRqSeVhDSrOWu2s2TjbkbqwUIicUklIc3qyWkFtMtMY/jJ3UNHEZEGUElIsynafYCX5m/gqvyetMmM/C05IlINlYQ0m2dmrKGkzDUkuEgcq1NJmNkVZvZtM+tbZfq3mieWxLuSsnLGT1/DsD6dOSa7beg4ItJAtZaEmf0KuBPoDbxuZndVmv21Zsolce71RZvYuGs/Nw7NCx1FRBqhLgeKLwFOcfdSM/sZ8JyZ9XD37wO6XEWq9eS01fTs0JqzP9cldBQRaYS6HG5KcfdSAHffClwI5JnZo3V8vySZpRt388HKbYwckkuqxmkSiWt1+U9+g5kN/OQHdz8IXAM40L+5gkn8enLaajLTUrhG4zSJxL26lMRNwPrKE9y93N1HAcOaI5TEr137S5g0t5DLT+pOB43TJBL36lIS5wFl1c1w9/ebNo7Eu7/NXkfxwTJu0AlrkYRQlxPX3YDnzSwFeAX4l7vPbd5YEo/Ky2PjNJ2ScyQDerYPHUdEmkCtexLu/kt3HwZcDqwGfmRmC81srJkNN7Os5g4p8eG9FVtYuWWvLnsVSSB1vjrJ3be6+1Pufg1wEjABOAOY1lzhJL488X4BndtmcNGAo0JHEZEm0qBLWN291N2nAO3c/aTK88zsjxWvv2+CfBIn1m4r5s0lm7h2UA6Zaamh44hIE2nsfQ4XVjPtlYp7KF5t5Loljoz7oIAUM746OCd0FBFpQrWeuDazaq9sIna3tVdZ9i9AK+AKIMvMrnJ3Dd2R4PYeKGXCjDVceMJRdD+ydeg4ItKE6nQzHXCUu6dW+Urh3++fuBlYAlwMLFFBJIeJs9exa38ptww7OnQUEWlidbkE9h9AX6ComnmvVDNtobtPNrMOjUomcaGs3HnsvVUMzDmSgTna5CKJpi6XwH7T3afWMO/WaqZNrHj9W2PDmVlfM5tX6WtXlVFoJbA3Fm+iYGsxo4YdEzqKiDSDSD8uzN2XAicDmFkqUAhMCplJDvXo1FX0OLI1F/TrGjqKiDSDBpWEmS0G5gEfVnx95O6FTZirOucCK9y9oJl/j9TRR+t2MGPVNn58yfGkpWpAYJFE1NA9iaeAUcRGgwVwM9tOrDDmAe8CL1WMGNtUriV2A98hzGw0MBogJ0eXX7akR6euom1mGtcM0mivIomqoX/+lQO7iN1x3QsYCvya2Anum4BHgDVm9uUmyIiZZRAbFuS5qvPcfay757t7fnZ2dlP8OqmDDTv38eJHG7h2UC/atUoPHUdEmklD9yT+A7ip0gntQmCGmY0FXgd+AeQB481sh7u/2cicFwFz3H1TI9cjTeSJ9wsod+fG0/NCRxGRZtTQkkgD/m2YT3ffaWa/BH7s7qeZWV/gx0BjS+I6qjnUJGHsPVDK09MLuKh/N3p11PiOIomsoYebJgH3mFmnauaVACdUfP8SkN/A3wGAmbUBzgeeb8x6pOno5jmR5NHQkvgeUAysMLN7zGygmfU0s7OA+4nddQ1wkNjwHQ3m7nvdvZO772zMeqRp6OY5keTS0FFgtxM7Wf0g8C1gJlAATCZ2GOq2ikXzgVWNjylR8aZunhNJKg2+ma7i8tafmdnPgROJPcFuC/Chu5dULPY+sUtiJUE8/O5K3TwnkkQafce1uzuf3VRXdd7bjV2/RMf0lVuZuXo7P72sn26eE0kS+qRLnT00ZTmd22Zw7Wm6aVEkWagkpE4+XLuDdz/ewqhhx9AqXU+eE0kWKgmpk4emLKd963RGDskNHUVEWpBKQmq1ZOMuXl+0iZs/n0fbzEgPHCwiTUwlIbX605QVtMlI5SYNwSGSdFQSclgri/bw4kfrGTk0lyOzMkLHEZEWppKQw/rzWytIT01h1Bd085xIMlJJSI3WbS9m0txCrjsth+x2maHjiEgAKgmp0dh3VmIGo8/QXoRIslJJSLU279rPMzPX8uWBPel+ZOvQcUQkEJWEVOuRqasoLSvntjOPDR1FRAJSSci/2bLnAOM+KOCyk7qT17lN6DgiEpBKQv7NQ5OXc6C0nG+f2yd0FBEJTCUhh1i7rZjx0wu4Or8nx2a3DR1HRAJTScghfvv6MlLMuPPc40JHEZEIUEnIpxZv2MUL8wq56fN5HNW+Veg4IhIBKgn51H+/upR2mWl888zeoaOISESoJASAGau2MXnJZm4761jaZ6WHjiMiEaGSENydX7+yhC7tMrn59KNDxxGRCFFJCG8u3szsgu3ceV4fWmfoqXMi8pnIl4SZHWlmE81siZktNrOhoTMlkrJy5zevLuHozm24Or9X6DgiEjHx8Jix3wOvuPtVZpYBZIUOlEhemFvIsk17eOirp5CeGvm/GUSkhUW6JMysPXAGcBOAux8EDobMlEgOlJbx29eX0b/HEVzcv1voOCISQVH/0/FooAj4i5nNNbNHzEyDCTWRR95dReGOffzwws+RkmKh44hIBEW9JNKAgcCf3f0UYC/wo8oLmNloM5tlZrOKiopCZIxLhTv28dDk5XzxhK4M65MdOo6IRFTUS2IdsM7dp1f8PJFYaXzK3ce6e76752dn6z+7uvrFvxbhOD+5tF/oKCISYZEuCXffCKw1s74Vk84FFgWMlBDe/biIlxds5PazetOzg64DEJGaRfrEdYU7gPEVVzatBG4OnCeuHSwt595/LCS3Uxa36rGkIlKLyJeEu88D8kPnSBSPTl3FyqK9/OXmQbRK141zInJ4kT7cJE1r/Y59/HHyx5zfrytn9+0SOo6IxAGVRBK578XFlJU79+hktYjUkUoiSUz9eAsvzt/A7Wf3pldHnawWkbpRSSSB2MnqBeR2ymK0TlaLSD2oJJLA2HdWsKJoL/de1k8nq0WkXlQSCW7R+l38/s2PueTEbpzzua6h44hInFFJJLCDpeX8x7PzaN86g18M7x86jojEocjfJyEN98fJH7Nk424euSGfDm0yQscRkTikPYkENW/tDv7nrRVcdWpPzuunw0wi0jAqiQS0v6SM7z47j67tMrnnMt0TISINp8NNCeiBV5eyomgv424ZzBGt0kPHEZE4pj2JBDNj1TYefW8VI4fk8IU+nUPHEZE4p5JIIHsPlPK95z6kV4cs7r7o+NBxRCQB6HBTAvmvfy1i7fZi/jp6KG0ytWlFpPG0J5Egnp21lmdmruUbZx7LaUd3DB1HRBKESiIBLFy/k5+8sIDP9+7Edy/oW/sbRETqSCUR53YWl3DbuNl0yMrg99eeQmqKhY4kIglEB67jWHm5851n57Fx537++vWhdG6bGTqSiCQY7UnEsT9NWc7kJZv5yaX9GJjTIXQcEUlAKok49c6yIn77xjK+dHJ3rh+SGzqOiCQolUQcKtyxjzufmctxXdpx/5UDMNN5CBFpHiqJOLP3QCm3PTWb0jLnzyMHkpWh00oi0nz0P0wcKSkr5/an57Bw/U4eviGfY7Lbho4kIglOJREn3J0xk+bz1tIi7r9iAOcer+G/RaT5Rb4kzGw1sBsoA0rdPT9sojAefH0Zz85ax7fP7cNXB+eEjiMiSSLyJVHhbHffEjpEKOOnF/CHycu5Jr8X3zmvT+g4IpJEdOI64l5buJGfvLCAs/tmc98V/XUlk4i0qHgoCQdeM7PZZjY6dJiWNLtgG3dMmMuAHu3504iBpKXGw+YSkUQSD4ebvuDuhWbWBXjdzJa4+zufzKwojtEAOTmJc6x+QeFObnliFt3at+KxmwbpUlcRCSLyf5q6e2HF62ZgEnBalflj3T3f3fOzs7NDRGxy89bu4KsPf0CbjDSe/NpgOmlMJhEJJNIlYWZtzKzdJ98DFwALwqZqXrNWb2PkI9Npn5XOX78+hJxOWaEjiUgSi/oxjK7ApIqTtWnA0+7+SthIzWfaiq3c8sRMuh7RiqdvHUy39q1DRxKRJBfpknD3lcBJoXO0hKkfb2HUkzPp1SGL8aMG0+WIVqEjiYhEuySSxZQlm/n6uNkc07kN40YN1nMhRCQyVBKBvTC3kO9P/JC+R7Xjqa8NpkObjNCRREQ+pZIIpLzcefCNZfxx8nIGH92RsTfk0751euhYIiKHUEkEsO9gGd99bh4vzd/INfm9+K8v9ScjLdIXmolIklJJtLBNu/Zz65OzmF+4kzEXH8+oYUdrqA0RiSyVRAtaULiTUU/MYtf+EsZen8/5/TTct4hEm0qihfzjw/X8cOJHdMhKZ+Jtp9Ov+xGhI4mI1Eol0cz2Hijlp/9YyHOz13Fqbgf+PHIgXdrpHggRiQ8qiWa0oHAn354wl1Vb93LHOb2589w+GslVROKKSqIZuDuPvbeaX7+8hA5t0hk/ajCnH9s5dCwRkXpTSTSxLXsO8P3nPmTK0iLOO74rv7nqRDrqBjkRiVMqiSbi7vxtTiH3vbiIvQfL+PnwE7h+SK4ubxWRuKaSaAIri/YwZtICpq3cyqm5Hbj/igH0Papd6FgiIo2mkmiEA6Vl/N/bK3loynIy01K474r+XDcoh5QU7T2ISGJQSTTQ9JVbGfPCApZv3sOlJ3bjnkv7aXhvEUk4Kol6Wr55N79+ZSmvL9pEjyNb85ebB3F23y6hY4mINAuVRB1t2rWfB19fxrOz1pKVkcb3LjiOW75wDK0zUkNHExFpNiqJWuzaX8L/vb2CR6euoqzcufH0PO44p48uaxWRpKCSqMH2vQd5YtpqHn9/NTuKSxh+cne+e35fcjplhY4mItJiVBLjx8OYMbBmDeTkUHjv/TzS+SSembGWfSVlnHd8F+467zj692gfOqmISItL7pIYPx5Gj4biYj7u1Iv/PeEK/r6oDaSu4vKBPbntzGM5rqvudxCR5JXcJTFmDBQXM6nfWXznsu/R+uB+rp/7IqPWz6THr+aGTiciElxyl8SaNQCcsWoOd00dz42z/0WH/btBQ2mIiACQ3ONW5+QA0GnfLu56b0KsICpNFxFJdpEvCTNLNbO5ZvavJl/5ffdBVpWrlbKyYtNFRCT6JQHcCSxuljWPGAFjx0JubuwQU25u7OcRI5rl14mIxJtIl4SZ9QQuAR5ptl8yYgSsXg3l5bFXFYSIyKciXRLA74AfAOWBc4iIJKXIloSZXQpsdvfZtSw32sxmmdmsoqKiFkonIpIcIlsSwOeBy81sNfAMcI6Zjau6kLuPdfd8d8/Pzs5u6YwiIgktsiXh7ne7e093zwOuBSa7+8jAsUREkkpkS0JERMIzdw+docmYWRFQUM2s9sDOWqZ1BrY0U7TaVJevJdZT1+VrW+5w82uaV5dtAuG2S6htUp/3NPV2qeu20mel4ctF9bOS6+7VH69394T/AsbWNg2YFaV8LbGeui5f23KHm1/TvLpsk5DbJdQ2Cbld6rqt9FlpuW1Sn23VXNslWQ43/bOO00Jpqiz1XU9dl69tucPNr2metknj39PU26U+2yoUfVbq9nuaTEIdbmoMM5vl7vmhc8ihtF2iR9skmppruyTLnkRdjA0dQKql7RI92ibR1CzbRXsSIiJSI+1JiIhIjVQSIiJSI5VEHZjZMWb2qJlNDJ0lmZlZGzN7wsweNjMN1xsR+nxEk5l9qeKz8lczu6Ch60n4kjCzx8xss5ktqDL9QjNbambLzexHh1uHu69091uaN2lyquf2uRKY6O63Ape3eNgkUp/tos9Hy6nndnmh4rNyG3BNQ39nwpcE8DhwYeUJZpYK/Am4COgHXGdm/cxsgJn9q8pXl5aPnFQep47bB+gJrK1YrKwFMyajx6n7dpGW8zj13y4/rpjfIGkNfWO8cPd3zCyvyuTTgOXuvhLAzJ4Bhrv7L4FLWzhiUqvP9gHWESuKeSTHHzjB1HO7LGrheEmrPtvFzBYDvwJedvc5Df2dyfpB68Fnf5FC7D+fHjUtbGadzOx/gVPM7O7mDic1bp/ngS+b2Z+J1l3AyaLa7aLPR3A1fV7uAM4DrjKz2xq68oTfk2gK7r6V2HE9Ccjd9wI3h84hh9LnI5rc/Q/AHxq7nmTdkygEelX6uWfFNIkGbZ9o0naJpmbdLslaEjOBPmZ2tJllEHuo0T8CZ5LPaPtEk7ZLNDXrdkn4kjCzCcA0oK+ZrTOzW9y9FPgW8CqwGHjW3ReGzJmstH2iSdslmkJsF43dJCIiNUr4PQkREWk4lYSIiNRIJSEiIjVSSYiISI1UEiIiUiOVhIiI1EglISIiNVJJiIhIjVQSkrTM7DQz81q+ys2sXZX3vVVp/l3VrPdLZvaamW01s4NmVmhmE83swqrL1pLvn2Y2/zDzHzKzHWaWaWY/rZRJT4iTJqOSkGS2HBha6WsCsLvKtEHuvrua906pmP9M5Ylm9iDwN2IDrI0iNlTzj4DWwMtmdmw98k0A+lf3YJ+KB81cBTzv7geARyryzK3H+kVqpaHCJWm5+zbgg09+NrOfAAvc/YOa3/WpbVWXM7PhwF3Aze7+eJXlnzKzy4B99Yj4d6AYuA74SZV5ZwNdiRUJ7r4OWGdmu+qxfpFaaU9C5DMDgAW1LlWzu4CZ1RQEAO7+T3dfX3mamQ0zs7fNrLji8NTDnxzeqnh+xj+p/vnE1wKbgcmNyCtSK5WECGBm7YmNyV/jOYBa3p9G7HDPa/V4z+eBN4CNxA4d3QVcDPyl0mITiA0DfWql96UDVxIb7VPP+pZmpcNNIjH9K14buifRCcjk0MdIYmYGpFaaVOafDb38K+B9d7+m0vKFwJtm1t/dFwAvAzuI7TnMrljsi0AHKg41iTQn7UmIxHxSEg3ak6ik6tj73wVKKn3dDmBmWcT2PJ41s7RPvoCpFcudCuDuB4k92/vqisKB2OGnAmLPFRBpVioJkZgBwCZ339LA928FDhB7dGRlTwGDKr4q60BsD+N/OLREDgDpHPo4yglADjDUzFoBw4FnXA+DkRagw00iMf1pxF6Eu5ea2TTgAuCeStM3AZsAPtsRAGKHkBz4KfBSNausfIJ7SsU6rgW6Ae3QoSZpISoJkZj+xP7qb4zfAS+Y2fXufth1ufteM/sA6OvuP69l2TIzexb4CtADWOzuHzYyq0idqCQk6ZlZN2Innht1PsLd/25mvwMeN7OziV2+uqVi3RdULLan0lt+QOwkdTkwkdiNfDnAJcAYd19WadkJwB3AFcC9jckpUh8qCZHY+Qho3D0SALj7d8zsHeCbwKPEDg0VETvJfLG7v1xp2almdgbwM2J7ManETki/QsUhqkrLTjOz1UAeOtQkLch07kukfszsLWInqq/h0EtagzGzFGIXorwJFLn7VYEjSYLQ1U0iDXMlsauR7gwdpMI9xPKcETqIJBbtSYjUk5n1JXYYCWCNu28OmQfAzLoD3St+3ObuK0PmkcShkhARkRrpcJOIiNRIJSEiIjVSSYiISI1UEiIiUiOVhIiI1EglISIiNVJJiIhIjf4/TIWow3Zd4c8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "T=np.logspace(np.log10(1E-1),np.log10(100))\n",
    "plt.semilogx(dfg['T'],dfg['geff_to_half'],'ro')\n",
    "plt.plot(T,gs_to_half(T) )\n",
    "plt.xlabel('$T$ [GeV]',size=15)\n",
    "plt.ylabel(r'$g_*^{1/2}$',size=15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 977,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import misc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 985,
   "metadata": {},
   "outputs": [],
   "source": [
    "m=100;Mp=1.22e19"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 986,
   "metadata": {},
   "outputs": [],
   "source": [
    "global k\n",
    "k=1\n",
    "fk=lambda σv:np.sqrt(np.pi*Mp**2/45)*m*σv\n",
    "def BEx(x,m=100):\n",
    "    T=m/x\n",
    "    return k*gs_to_half(T)/x**2\n",
    "def Yeq(x):\n",
    "    return 0.145*x**(3/2)*np.exp(-x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 987,
   "metadata": {},
   "outputs": [],
   "source": [
    "def BE(Y,x):\n",
    "    return -(Y**2-Yeq(x)**2)*BEx(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 990,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import integrate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 991,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEkCAYAAAAfLVccAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABd1ElEQVR4nO3dd3gUVffA8e9NQkghtACioReR3qQpHURAiiCCqCiCIi8o6k99EXlVFAVRRECwgaCANAEB6U1AilRRSgDphJ4CaaSf3x+7GbIh1JTdJOfzPPuQc6fs3UzIyczcOdeICEoppVRWcXN2B5RSSuUumniUUkplKU08SimlspQmHqWUUllKE49SSqkspYlHKaVUltLEo5RSKktp4lEqExljKhhjvjPG/GOMSTTGrL/BelWMMWuNMdHGmLPGmI+MMe53u94N3uNxY8wqY0yIMSbOGHPGGDPPGNP2Dj/Tb8aYvTdZPsEYc9kYk/dO9qtyD008SmWuqkB74BBwOK0VjDGFgDWAAJ2Bj4A3gQ/vZr0bvMeXwHzgDPAi0Bp4B/AGlhtjyt/BZ5oFVDPGVEnjfdyBbsACEYm9g32qXMRo5QKlMo8xxk1EkuxfzwOKiEjzVOsMAf4LlBaRcHvbf4FhQPEUbbe1Xhp96AwsBF4QkR/TWN4R2CUiZ2/zM/kCF4ExIvJeqmWtgdVAGxFZfTv7U7mPnvGoHMsY42aMiTTGvG6MGWOMuWiMCTPGvG1f3ssYc8C+zgJjjHdG9yE56dxCO2BlqsQxG9vZSLO7WC+114EdaSUdex9/S510jDFNjDEb7Jf0Qowxk4wxfvb1o4DfgB5p7O4pbElp3U36o3I5TTwqJysH+GL7xRsLPA0sBT4zxkzA9kvybWyXnB4HXki9A2PjcatXOvv5AHAwZYOInAKi7cvudL2U/fcAGgGrbrczxpiHsV3SO4/tstnr2C4XTk2x2iygojGmbort8gBdgbkikni776dyn/T+h1HKlVW3/ztGRMYDGGP+BZ4BKgOtxX6t2RjTD6iUxj6ex/EX7o2YdPSzEHA5jfYw+7I7XS8lfyAvcDplozHGACkHJSTKtevunwJbRKRHivXPAGuNMdVEZB+w3N6Xp4Bd9tUetfdj1g36ohSgiUflbDWw/XL8JkWbr/3fT8XxBqcvEJrGPn4D6mVK77JW6pu5bwKfp4hfBSYYY3ywnSG9mupMbhMQD9QF9olInDFmAdDdGPNf+/eyB3AS2JpZH0LlDJp4VE5WHdgkIvEp2moACcDG5Ab7L9sywL409hEKXMnEPoLtjKVAGu2F7MvudL2UQrBdZiyRqn06sN7+9Y5U+3IHvra/UiuZ4utZQB+gkTFmN7aRdl+nSuhKXUcTj8rJqgNzU7XVBA6mGupbHdv9zn/S2EdWXGo7SKp7NMaYkoAPjvd0bnc9i4gkGGO2Am2A91O0XwAu2PeRcpPL2M6OhgHL0thlykEIv9v38RRwL+CHXmZTt0ETj8qR7CPUKgB/p1pU4wZtUcDRNHaVFZfalgNvG2P8RCTC3tYDuApsuIv1UhsLLDTG9BKR6TfriIhEGWP+BCqJyEe3WDfRGDMXeBIIAAJFJPX3VqnraOJROVVVbGcxaSWZ8Wm07U9r6LOIhGC7XHVX7Jfx2tvDACC/MaabPV4mItHAt8AgYIExZhS20XjDsA2KSDl0+nbXS/0ZFhljxgI/GmNaYEumwdgGHrSxrxaZYpP/YhtIkATMAyKAUsBjwFARSfkg7Cxs94e6AB/czvdEKX2AVOVIxpgXgK+A/Cke4CyE7Z5NWxFZmWLdDcBhEXkpE/pRBjh+g8VlReSEfb0qwARsN/YvA5OBYamHJd/uejfoSxdgALYBAn7AJWwDASaLyPJU6zbAVhHhIWz3fE4CK4APReRKqnWPY7tHVlFEjtyqH0pp4lFKKZWl9AFSpZRSWSpH3eOx15D6GogD1ovIz07uklJKqVRc/ozHGDPFXmNrX6r2tsaYQ8aYI8aYd+zNXYF59mv1nbK8s0oppW7J5RMP8CPgMF+IvfT6RGxFE6sAPe03XUtwrTSI1opSSikX5PKJR0Q2cn0pk/rAERE5JiJx2Cr0dgaCuPaEtst/NqWUyo2y6z2eAByLHgYBDbA9nzHBGPMYtmcVrmMvBtkPwMvLq26pUqUyuasqsyQlJeHmpn9fZFd6/LKvw4cPB4tI0bvdPrsmnjTZ5wm5rrR9qnW+B74HqFSpkhw6dCgruqYywfr162nevLmzu6Hukh6/7MsYczI922fXPzfO4FissIS9TSmllIvLrolnB7ZJqMoaYzyxFSlc7OQ+KaWUug0un3iMMbOwlfWoZIwJMsb0FZEE4BVgJRCIbcbD/c7sp1JKqdvj8vd4RKTnDdqXkXbZdqWUUi7M5c94lFJK5SyaeJRSSmUpTTxKKaWylCYepZRSWUoTj1JKqSyliUcppVSW0sSjlFIqS2niUUoplaU08SillMpSmniUUkplKU08SimlspTL12pzJSHBwZw6fZoy1aoRk5TExOHDiYuMpGu7dmAM85cuxcfPjz7/+x8A8ydMoESxYrR55BHyFymCW548Tv4ESinlfLk68Rw9epS8Hh50rVyZqKgolh4/ThLQaNcuwhMS2N+gAQC+QCyQkLzh77/b/v3kEwA+HzfOYb8ftW5t++Ltt9N834DNmyng4cGBBg1wA9qWKIGvnx+LDh+mcKFCTAoMJCBvXnZMmcJzffrg5eeXkR9bKaWcKlcnnsTERBKBxfv2kc8YkuztPm5u3Ovjw2EPDxITEmhfrRoenp4s2rcPjOGT8uXxcnNjWJEiJMXH06dDBxBh8pIleHp68kmlSgjwf76+uAFd69cnLjqamdu3kydPHtoWLszlhAQOAALsOXeO8KAg4oDzwcF03LfP1pHXX+fl11+npJsb93p7syMqiqKFCjHp8GHu9/bGHD5Mpdq1nfCdU0qpu5erE09AiRJs+P13yleokPYK8fE33b7/pUsO8aeplveJjHSIp6XegYhDGHjgAH8fPEjZOnU4Hh7O8x4e+Hl6UsnfnxPBwQgQfPkynfftI8/27cQPHowb0K5wYe4tVYrVZ8/y4fvv8/zAgTftt1JKOVOuHlzg6+Nz46TjBJWrVOGprl1pkD8/T5UoQWx8PMFRUaw+dYp/o6MREQ6fP8+fderwSr58GCC/pyc7IyKYvGcPJy9epN///R/l/vyTFq+/jneePLzdrRsJZ3RWcKWU68jViSc7Kl+sGA3y52fMgAEkiRAWG8v5uDjmL11K6VKl6DBgALXz5WPLlCnEJCQwev58ipUoQWkPD3w8PPhpwgRnfwSlVC6Xqy+15SRd27en68mTVhx76RK9+vXDz9OT49u28fvevQC8tmcPo7ZvJ6FnT8zly2ybNImC7duDl5ezuq6UymU08eRQefPmZe5PP1lxbGwsfV97jbpDh7IiNJRVx44BUOmJJ3jc3Z3Ie+9l9ODB3DtwIBjjrG4rpXIBvdSWS+TNm5cZ337LGyVLsrJmTZatW8f9depQrHZtpgIzg4K479VXeevoUQ5ERcGqVXDhgrO7rZTKgTTx5FLtWrTg0K5d7N29m+MXLlCwWDHyFS/OuDNnqD5rFubRRyldvDiXW7aEadMgIsLZXVZK5RA5KvEYYx43xkwyxswxxrRxdn+yiwB/f8IuXCDi3DnONGpEPfsDsqeAUuvX0+355+lXsCD8/LNT+6mUyhlcJvEYY6YYYy4aY/alam9rjDlkjDlijHnnZvsQkYUi8hLQH+iRmf3NqYp5evLnzz8TGxvL18uXU/DRR5kPTEpKotbMmWwLD4f16+E//4HNm697FkkppW7FZRIP8CPQNmWDMcYdmAi0A6oAPY0xVYwx1Y0xS1K9iqXY9H/27dRd8vT05D9t23Jq+XLe+uADPLy9OfHOOzTcvZu8jzxClW+/hcaNoXx5eO89OHjQ2V1WSmUTLpN4RGQjEJqquT5wRESOiUgcMBvoLCJ7RaRDqtdFYzMKWC4iu7P6M+RUnw8bRnx0NEGNGjHUz4+4hAQCgfvLl2d14cIkfvIJtGoFSfaiQzExTu2vUsq1ufpw6gDgdIo4CGhwk/VfBVoDBYwxFUTk29QrGGP6Af0AihYtyvr16zOut7lAa6DIN98wYd48juzaRZujRwFoW6oUgzduhMREGj79NNGlSnHhkUcIbtyYRB+fTOlLZGSkHr9sTI9f7mXEha7RG2PKAEtEpJo97ga0FZEX7XEvoIGIvJIR71epUiU5dOhQRuwqV7oYHk758uWJDA6Grl35v1GjGOTpSanvvsPMnAknToC3Nzz+uK1SdwYXNF2/fj3NmzfP0H2qrKPHL/syxuwSkQfvdntXP+M5A5RMEZewtykXUCx/fiIuXWLNrl3M9PXly6AgxrRogXuePMRGR+O+bRvMmAFz5sDzz9s2OngQtm+H9u2hSBHnfgCllFO4zD2eG9gBVDTGlDXGeAJPAYud3CeVSuu6dZnywAP8WrAgAInx8RR9+GFW5ctH6CefwPnzkDxH0bx5tiR0zz3QpAl89BFs2gSJic77AEqpLOUyiccYMwvYClQyxgQZY/qKSALwCrASCATmish+Z/ZT3VjnWrWIjIqiy9ChhB04wKO1a+Pv78+chQvB3d220rvvwo4dMHQoREXBsGHQtu21gQlr1sC6dXDlirM+hlIqk7nMpTYR6XmD9mXAsizujrpLvj4+LPj4Y/5+6SVqly+PJCbyzsGD1I6OpkhMDIUKFcI8+CA8+KDtbCc0FAIDIXla8P/+F/76y/b1Aw9AvXq2xPT00877UEqpDOUyiUflLDVLlyYpIYEB48czu04dqmzeDI89xhNdujBnzpxrKxYuDA8/fC1eswZ27rTdB9qxA1avhujoa4mnenUoWBCqVqWEmxuEhUHNmlCuXJZ+PqXU3XOpUW1ZTUe1ZY0zsbHc/8ADRJ84gZuXF0dOnqR0kSKEh4dT0H5f6IZEbM8FeXtDQgIMHAgHDsD+/bakA/D66/Dll7b1mjSBMmXgvvtsr3vvhYYN4f77r1VZ0OrbLkFHtWVfOX1Um8oBAvLm5fLhw9xbsSIhZ85wf61adO3YkbXz57Nlyxbuv//+G29sjC3pAHh4wHff2b4WYdNvv9E4IAAKFLC1RUSAvz/s3Wurrh0ebmsfM8aWeI4cgRo1bAMb/P2vvQYMsCWsS5dsZ1z+/rYzsQIFIH9+29fJlwKVUummiUdliTx58hB84gST16zh5SefZO7335PHy4vy5csDEBYWRqFChW5/h8aQkD8/1K17ra1oUVix4locFQXnzl1LTD4+8MorcPEihIRAcDAcP277GmDfvrTvJS1cCJ0725LSSy/ZklHyq0ABW8mgypVtCW/JEtv7JL+8vaFlS9vlweBg23v7+l5b7uV1beCFUrmEJh6VpV5s3Rq/adN4qlMn4uvUYcCRI3x6771Ur16dl156iQ8++CDj3szXFypUuBYHBMDnn994/YYNbQMdgoNtgx7Cw22vmjVtywsVsp0ZXbliaz9/Hg4fhqtXbct37rSN2kvtn39siWf2bHj11euXHz1qu0c1diyMHAmenpA3r+3l6Qm//27bfupUmD//WnvyOhMm2M7IFi2y3Rvz8Lj28vSEN9+0vc+6dbb3Srncx8eWVME2qCM4+PrlyZ//2DHb/TY3t2svLy8oVcq2/MIFiI+3tbu72/719LyW+KOibP/atzUJCbZh9Jp4cx1NPCrL9ejYkW6Jibx34gQjT51i6sMP81CtWjzyyCOArZRKUlIS+fPnz9qOeXvbRtLdSN26trmJbqR3b9sZU3S04ys5+bVta0s+ye1RURAba7u0B7b37trV1hYbC3Fxtn+TL/Mln8GlXBYbC19/bVu+di18843tXlgyL69riWfqVNsDvSkVLXot8Xz4oS15pVSunC1Zge1sb906x+U1a8KePbavO3WyJb6UGjeGP/6wfV2vni2x2zUDaNcOltkHrZYtC0FBjomta1eYPt22vGRJW9I35tp9ul694KuvbF8XK2ZLZCmXv/wyfPyx7Y+D5AEoKZe//rqtqkZwMNSpc/3yIUNs+zh1ynbmmnr5sGG2Y37wIHTpcq09eZ0RI2zfl927rz1EndLYsbY6h5s22S75pvb997Y/iFautI34TG3GDNuAm19/tY0STW3+fNvnnjkTRo++fvmyZVC8OEyebPvZSe33321n9l99dfOf/TukiUc5hbubGyPKlePi4sX8EBLChrVriR4zBoAxY8Ywbtw4Dh06RJHsVN3AmGtnIWldNqxQwfEMLLW2bW2vG3nlFdvrRsaPt71EbM9FJSQ4JqGvvoJPP73WnnIZwKhRtl/CKZd7el5b/sEHtukwkpKuvZLPZsD2bNbFi7Zf/snLixe/tvytt2yXNe3Ljh09SrlWra4tHzAALl+2LU9MtH2OGjWuLe/Vy5ZAkgeJiECDFKUbn3rKtm3K5clna+7utgQg4rg8+XjkyWN7yDn18hIlbF/nzWtLAKmXFy1q+9rLy5YAktuT10n+/nh52e4zppYvn+1fH5+0fzaS6xzmy5f2yM28ea8tTz7zTMnD49ry5M+SUvLZZr58tsE4qbm5XVterNj1y++SjmrTUW1O92j37qxesQKMYeyPP9LgvvtYt24dQ4YMAWD58uU8+OCDFE3+T26no6KyNz1+2ZeOalPZ3sq5c1m3fz9tWrXita5d8StUiCv2G/6RkZE8+eST9OzZk0mTJjm5p0qpjOAyJXNU7tayalUmT7TN3Rfh7s6vwcEA5MuXjx07dvCu/ab9v//+S506ddi5c6fT+qqUSh9NPMpl9H7iCf49d45GK1fy5P79dHjtNX777TcqV65M2bJlAQgJCSFPnjzce++9AGzevJmxY8cSHR3tzK4rpe6AJh7lUioUL87qmjUp+eOPLB0/nie6dycpuYAo0LBhQ7Zt20ZAQAAAy5YtY9iwYeSxj/xaunQpv/zyi1P6rpS6PZp4lMvxdXdn1QcfYNzdiY+JYeDXXxMbG0tMGlNqf/LJJxw5csRKPBMnTuTTTz+1lo8ZM4apU6dmWd+VUremiUe5pPvLluViWBgF69fn29de476yZWnVqhXx8fHXrZtyyPXixYtZvPjalE0LFixg5cqVVty2bVs+SvG8w9atW7lw4UImfQqlVFo08SiXVcTPj8Nr1uDu40PouXOcCQuzzmxuxMPDw7oMB7Bp0yam2x9AFBHuuece/O0PbCYlJdGyZUs+++wza3m7du2s6tkiwtq1azUxKZXBNPEol1bUz49dO3aQp2hRgsaO5Y/Ll9m7dy8LFiy47X0kJytjDD/99BMDBw4EbIll2bJl9O3bF4CoqCguX77MVXsJnODgYFq3bs3cuXMBuHjxIiVLlrQSU2hoKEOGDGHfvn3W9hs3biQ0NNTav1Lqepp4lMur+cADnA0KokKBAjz69dc81qkTr732Wpr3fO6Eu7s7LVq0oEqVKoBt6PbWrVvp3bs3APnz5+f333+ns72kjIjQunVr64zq7NmzjB49muSHkA8dOkSzZs34w14iZvv27Xh5ebF69WoA/v77b9q0acM///xjrT9kyBBOnToFQFBQEHPnziXMPt3DlStXOHjwIHFxcQAkJCQQHx+vCU1le5p4VLZQxNOTlTVrcvV//+P0iRO8O2oUXl5eJCUlEWx/5iej5c2bl+bNm1PKXorknnvuYerUqTRu3BiAatWqERsby+OPPw5AhQoVWL16NQ899BAARYsW5bXXXrOGgsfGxhKePFUDcPToUcaMGWP1/88//6RHjx4EBQUBthF6lStX5vjx4wDMmDEDT09Ph7hUqVKcP38egPnz59O8eXOu2KcNX7p0Kb169bLO4FavXs2bb75Jgr1Uzh9//MGoUaOsRLZjxw6mpajHtX//fof7Y0ePHmXHjh1WfPr0aQ4ePGjFFy5c4PTp01Z85coV6+wPICYmhtjY2Nv75qucTURy7ev+++8Xlb28/tFHQp48krd8eZm7ZImMHTtWihQpIsePH3d21+5aUlKSiIiEh4fLvn375OrVqyIicuLECZk5c6aEh4eLiMiePXtk+PDhcuXKFRER+f333+WFF16w4l9++UWaNWsmkZGRIiLy3XffSdmyZa39ffrpp+Lr6ysJCQkiIjJ06FBxd3e3+vHWW2+Jt7e3FQ8aNEgKFChgxf3795dixYpZcZ8+faREiRJW/Mwzz0i5cuWs+Mknn5QHHnjAijt16iQ1a9a04gYNGkjDhg2tuH379vLII49Ycbt27aRLly5W/Oijj8rTTz9txW3atJG+fftacevWreWVV15xiN9++20rbtmypbz//vsO8YgRIxziMWPGWHGzZs1k4sSJVty4cWOZNGmSiIgkJCRI48aN5ccffxQRkatXr0qTJk3k559/FhHbsWzWrJn88ssvIiISEhIizZs3l4ULF4qIyIULF6RFixaybNkyEREJCgqSVq1ayerVq0XEduxbt24t69evFxGRf//9Vx555BHZvHmziIgEBgZKmzZtZPv27SIisnfvXmnbtq389ddfIiLy119/Sbt27WTfvn0iIrJ9+3Zp3769HDx4UEREtmzZIh06dJCjR4+KiMjGjRulY8eOcurUKRERWbdunXTq1EnOnTsnIiKrVq2Sxx9/XC5duiQiIsBOScfvXi2Zo7KVL997j/uqVeO/3brxbN++TBkzhhdeeIHSpUsDkJiYiHs2K7Nv7BWN/fz8qFq1qtVeunRp63MB1KxZk5rJRS+B5s2bO9Q669atG926dbPifv360a9fPysePHgwgwcPtuJhw4ZZFSEAhg4dyispipC++eabPPfcc1Y8aNAgevToYcUDBgzgiSeesOKXX36Zrl27WnHfvn0JSZ7rCHj++ecdzvhatGjh8Pk6duxIYmKiFT/yyCMOg0maNWuGn5+fFTdq1MhhRGOdOnWss1OAypUrO+y/XLly3HPPPVZ83333OcwBVahQIXySi3Jiu/TqmaJIqpeXFx4e135lenp6OvyseXh44OZ244tISUlJ1tmliJCQkGA9oyYixMTEWJ8/KSmJ6Oho6+w0KSmJyMhIK05MTCQ8PNyK4+PjCQ0NteK4uDiCg4OtUaBxcXFcvHjRWh4TE8PZs2etODo6mqCgIGv9yMhITp06ZS2PiIjg2LFjDscnPbRIqBYJzZYefvpptsyaRR5vb+LsVQtCQkJo2LAhX375JR06dHByD9WtaJHQ7Cu9RUJz3D0eY4yvMWanMUZ/8+Rgm2fOJH/58sSPHs2P584Btr/aypUrR8mSJQHbKLOE1KX/lVJO5zKJxxgzxRhz0RizL1V7W2PMIWPMEWPMO7exq8HA3MzppXIl8ydPptVDD9F3xw6atGuHn58fK1eutC5H/e9//6NatWrWqDCllGtwmcQD/Ag4zIJljHEHJgLtgCpAT2NMFWNMdWPMklSvYsaYR4ADwMWs7rzKeh7A3KpVke7d2bRiBV+mmkGxZcuWPPvss9Z1+nHjxrF582Yn9FQplZLLJB4R2QiEpmquDxwRkWMiEgfMBjqLyF4R6ZDqdRFoDjQEngZeMsa4zOdTmaNwnjwMGDIEjOGbZctISHHzs2PHjvzvf/8DbDdThw8fziL71M5if3g0MjLSKf1WKjdzqcEFxpgywBIRqWaPuwFtReRFe9wLaCAiN5n/F4wxvYFgEVmSxrJ+QD+AokWL1k1+Kl1lP5GRkeSzTx380cKF/D5uHOXbtKFuwYK8/PLL140wio2NJS4uDj8/P44ePcqLL77I22+/Tfv27YmIiODQoUNUrVoVb29vZ3ycXCfl8VPZS4sWLdI1uMDpz9KkfAFlgH0p4m7A5BRxL2BCRr2fPseTvf3+++/W14mJiVKkdm0BJI+XlwQFBd1029jYWFm7dq31XML8+fMFkD///FNERP755x/58ssvJSwsLLO6n+ulPH4qeyGdz/G4+qWoM0DJFHEJe5tSDtzc3NiycCHGywuvF1/Eu1ixm67v6elJy5YtredA2rRpw6pVq6hduzYAv//+O2+88Yb13MW0adNo2bIlUVFRABw7dox//vnHWq6Uun2unnh2ABWNMWWNMZ7AU8DiW2yjcqmKpUqx8+JFYrp14+XDhxk5ciSTJk26rW3z5cvHI488Yg1EGDRoEBcvXrQeMExKSsLNzc16wHD8+PFWaRyAL774gqeeesqKt2/fzoYNG6xYE5RS17hM4jHGzAK2ApWMMUHGmL4ikgC8AqwEAoG5IrLfmf1Urq2Onx/v+vszr2FDhg0fnq5RbEWLFrW+7t27N2vWrLGqDAwcOJB58+ZZcVxcnEMdss8++4z//Oc/Vvzkk0/SpEkTKx41ahQjRoyw4lWrVjkkqnPnzlk115TKaVymZI6I9LxB+zJgWRZ3R2Vj79x/Px/FxRGXkEDP118H4OrVqxk6aKBixYpUrFjRiocMGeKwfOzYsQ4FMtu3b++QSP7++2+Hh1vff/99ChQoQLNmzQDo0KED9957L0uW2MbHtGnThsqVKzNu3DgA+vfvT5UqVRg0aBBgS3SVKlWyKmn/+uuvlClTxrp0uGfPHooVK8Z9990HwKVLl8iXL58OpFBO4TJnPEplFC9PTzbu2YPx8aHb888TduUKzZo1Y+jQoVnWhxIlSlCjRg0r7tOnD2+88YYVz5w5k5QjKufNm8d3331nxanrptWuXZsKFSpY8cmTJx0mqBs3bhzLly93eL+UU343btyYMWPGWHFAQIA1E6uI4O3tbZ2BxcXFUbFiRb7//nvAVgGiWbNm/PLLLwCEh4fTpUsX6/3CwsLo27evNR1EaGgob7/9Nrt27QJspYxGjhxJYGAgYJvn6Ouvv+acveLEpUuXmDp1qlWV+9KlS8yePdv6fJcuXWLRokVW3bdLly6xatUqK5EHBwfzxx9/WEPjQ0ND2bVrl1WVOywsjH/++cc6Iw0LC+PAgQNWXbKwsDAOHTpk/SFw+fJljhw5YtUlu3z5MsePH7fqql25coXTp09bl0/Dw8M5e/as9b2NiIjg4sVrjxJGRUU5/BESHR3t8EdITEyMw7D+mJgYou1loJKPR8opQFKfXcfHxzs8JB0fH+8wU29CQoLDHzkJCQkONdcSExMd4qSkJOuzwrUBaBkqPSMTsvtLR7Vlb7caFdV/wgQBpEiZMjJo0CBZtGhR1nTMSRITE62vDx06JGfOnLHiJUuWyN69e614woQJsmXLFhGxVVoePHiwrFq1SkRslZZ79uxpfb+SKy3PnTtXREQuXbok1atXl9mzZ4uIrbJyiRIlZObMmSIicuTIEfHx8bEqNe/bt08Aa/vdu3cLIMOHDxcRkT///FMAq1Lzxo0bBZA1a9aIiMjq1asFkD/++ENERJYuXSqAbNu2TUREFixYIIBVmXn27NkCyIEDB0RE5KeffhJAjhw5IiIikyZNEkBOnz4tIiITJ04UQC5cuCAiImPGjBFALl++LCK2qt6AREdHi4jIhx9+KMANq3y//fbb11X5LliwoBWnrvL9wgsvXFflu3z58lb85JNPSuXKla24c+fODlW+27ZtK/Xr17fili1bSuPGja344YcfllatWllxvXr1pF27dlZcs2ZN6dy5sxU/8MAD0r17dysuV66cPPvss1YcEBCQ7lFtTv/l78yXJp7s7VaJJzExUTz8/ASQXzdtstoXLVpklZtXmS8pKUliYmIkPj5eRETi4uLk/PnzsnLlShERiYmJkePHj0tUVJSIiERFRcmBAwes6R3Cw8Nl165dEhERISIioaGhsmnTJmu6iAsXLsjq1autOCgoSBYvXmzFJ06ckHnz5lnbHz16VObMmWO93+HDh2XmzJnW9BEHDhyQ6dOnS2xsrIjYhtZPnTrV6v/u3btl8uTJ1nQW27dvl++//976vFu2bJHvvvvOijds2CDffvutFa9Zs8YhXrFihcP6S5Yscdjfr7/+ak3HICIyb948+eGHH6x49uzZMnXqVCueMWOGTJs2zYp//PFHmTFjhhX/8MMP1h8JIrbpM+bMmWPFEydOtKZzEBEZN26cLFiwwIq/+OKLdCcel3qANKtpders7XaqG6/asYPHN2+mYfPmrK1Zk8TERB588EEKFy7M2rVrrcEBKutpdersS6tTK3UTberVY2z37vweEkKb/v0JDw9ny5Yt/PzzzxhjCAkJ4fvvv9cq1kplIU08Ksd76d57yf/GG6z5/nteeestfHx8uPfeewHb9NH9+/e3ppNWSmU+TTwqxzPGMHncOHB3Z+Pp0w7LBg0axLZt26yh0YMHD7YKiyqlMocmHpUrPNmiBS1ee40za9YwbNo0Dh8+DNiSUr169QDbQJvUD26OHz+evXv3OqXPSuVUmnhUrrFg+HDcixXjw9696dipk8OzCmBLQtOmTWP8+PGA7fmTN998k2XLbM8vx8TEMGHCBE6nOmtSSt0ZTTwq1yjo48OLr78OxlCua9frpk1IljzSzd/fn4sXL9KvXz8A/vrrL1599VV2794NwNGjRxk8eDAnT57Mkv4rlVNo4lG5yrdDhtBh61Y2tm3LyZiYW569FCpUyCoU2qhRI06cOEHr1q0B2L9/P19++aX1lPlvv/1G3bp1OXXqFACnTp1i27ZtDk+RK6U08ahcaEKNGiRevUqjtm0pV64c69atu+1tS5cuja+vLwCdOnUiMjKSSpUqAeDh4UHRokWtqRZmzZpFw4YNrdItU6dOpV27dlYi2rdv3x29t1I5hSYeleuU9vLigcmTObdhA9UaNKB+/fp3vS9PT0/rkl27du1YsWKFNXXCs88+y9KlS8mfPz9gq5F19epV8uTJA8C3337LE088Ye3rtddeo2rVqlY8ZcoUPvjgAyvetm2bQwXr8PBwhxpdSmUXmnhUrrRy8mRMwYIcuXwZHx8fYmNjOXHiRIa+R0BAAO3bt7fil156ifXr11vxu+++61DYs2HDhnTt2tWKt2/fzsqVK6145MiRVjVqgJ49e9KoUSMr7tevHy+++KIVDx8+nNGjR1vxjBkzWLz42nRWW7ZsYf/+a7OMnDhxgkuXLllxXFwcubmyicpE6am3k91fWqste0vv1Mm97cUgB4wfL7169ZISJUpY9bxc0YkTJ2Tfvn1WPH/+fIcaXEOGDJF33nnHih9//HHp1auXFdeoUeO6YpBPPvmkFZcrV06eeeYZKy5ZsqS88MILDtsPGTLEih955BEZN26cFXft2tWhZljv3r2tGl+JiYnyxhtvyOrVq0XEVq+tb9++1lTjMTEx8vXXX8s///wjIiLR0dHy448/yqFDh0TEVr9t7ty5cvz4cRERiYyMlCVLllhTnEdERMjatWutQp+RkZGydetWCQ0NteI9e/ZY9dsiIyMd6sFFRkbKoUOHrEKgkZGRcvToUYmJibHe/+TJk1b9tujoaDlz5oxVvy06OlouXLhgFQ69evWqBAcHW4VbY2Ji5PLly1YcFxcnkZGRVr23+Ph4iYmJseLExESJj493iFMWgXU2tEioJp7cKr2JJzY+XjxLlBCMkW++/15++umnjOmYi4qKipIrV65Y8a5duxwS2YIFC2Tjxo1WPHbsWFm8eLEVv/baazJ9+nQr7tSpk3zzzTdW3LBhQysRJSUlSZkyZeSLL74QEZHY2Fjx8/OTzz77TERshT8BGT16tIjYKl4DMn78eBEROXv2rABWMc3jx48LIFOmTBERW/VtwKqA/c8//whgFbfcuXOnAFaF7c2bNwsgK1asEBHbzw4g69atExGRlStXCiCbN28WEZHFixcLIDt27BARW2FOwEqMM2fOFEAOHjwoIiJTp04VQI4dOyYitsKbgFUhfPz48QLIpUuXRERk9OjRAliJ8JNPPhHASnQffPCBAFbiGTJkiOTJk8f6Xr/55pvi6+trxa+++qoULlzYil9++WW57777rLhPnz5StmxZK3722WflgQcesOIePXpIrVq1rLhLly7SoEEDK+7QoYM0bdrUitObeFxmIjilspqnhwdDP/2UD/7zH/7y8+M7+9TV27dv5/jx4/To0cPJPcxYyfeektWpU8ch7tKli0P82muvOcRjx451iBctWuQQb9261fraGONQhsjT05Pw8HArzpcvH6tWrbKKhBYuXJhz587h5+cH2GZ/PXbsGP7+/gDcd9997N27l4CAAABKlSrF9u3bKVeuHADlypVj/fr1VKlSBYAKFSqwfPly6zNWqlSJ+fPnU6tWLQCqVKnC7NmzrXtq1atXZ8aMGVYFi1q1avHjjz9SpkwZAOrWrcsPP/xgvX+9evX47rvvuOeeewDbZdKJEyda/X3ooYcYN26cdX8veT6k5IEpTZs25bPPPiNv3rwANG/enBEjRuDhYfuV3LJlS9zd3a3vV+vWrR2OX5s2bShcuLAVP/roo9Ykf2CbeDD5ewPw2GOP8cADD1hxhw4drEkCATp27EjDhg2tuHPnzg7Hq3Pnzg5zAKWXVqfW6tTZVkZVN26/Zw9bIyM53qABBfPkoUuXLuzdu5d9+/bh5eWV/o6qNGl16uxLq1MrlU4jypfn8qlT1G/bltOnTzN37lzWrl2Ll5cX8fHx/Pbbb+TmP9CUymiaeFSuV8vPj1KzZvHvunX8MHMmefLkoXTp0gBMnz6dTp068eeffzq5l0rlHJp4lAJ++uIL8PRk9qZNDu3PPfcc8+bNs65/T5o0iR9++MEZXVQqx8hRiccY42aM+cQY85Ux5nln90dlH83vv5/qzz/PoSVLmDRnjlUGx8PDgyeeeMKq3zZ//nx+/fVXa7tNmzYRGRnplD4rlV25TOIxxkwxxlw0xuxL1d7WGHPIGHPEGPPOLXbTGSgBxANBmdVXlTNNGzYMfHzo17Onw4OXKS1fvpyZM2cCEBERQatWrazqAiLCzp07dTZTpW7BZRIP8CPQNmWDMcYdmAi0A6oAPY0xVYwx1Y0xS1K9igGVgC0i8n/Af7K4/yqbq3XffdR78UWoUIG2zz6b5jrGGGuIrI+PDytXruTll18GIDAwkHr16jF9+nQAQkND+fXXX7l8+XKW9F+p7MJlEo+IbARCUzXXB46IyDERiQNmA51FZK+IdEj1uojtLCfMvm1i1vVe5RTzP/2UPJMm8WNCAiJy07MXd3d3mjdvzv333w9AiRIlmDVrFm3b2v5++uOPP+jatSuBgYEA7Ny5k0GDBnH+/HnAVrtNR8up3MjVHyANAFLWrQ8CGtxk/QXAV8aYJsDGtFYwxvQD+oHtIbWUtbNU9hIZGZkpx68N8P3mzazr2pU6NWpY8/HcjuLFi3Po0CEOHTqEt7c3EydOJCIigvXr17NmzRomT55My5YtKViwIIsWLWLy5Mn8/PPP5M+fn3/++YeDBw/SpUsX8uTJQ1xcHHny5LHuL+U0mXX8VDaQnrIHGf0CygD7UsTdgMkp4l7AhIx6Py2Zk72lt2TOjRyKiBDuuUd8ixaViRMnZui+k5KSrDIo69evl4EDB1o1uN59913x9PS04sGDB0v+/Pmt9adOnSqDBg2y9rVnzx6rxIuIuFQtr9uRWcdPZT7SWTLHZS613cAZoGSKuIS9TalMc3++fNTv35+oS5comKIMSUYwxlhnMM2aNWPChAnWtAoff/wxly5dsuKWLVvyzjvvWOsfPnyYLVu2WPsaNWoUvXr1suKePXtaJWEA3nvvPYdq1nPmzGHOnDlWvH37dv766y8rDgkJ0RF6KmukJ2tl9Ivrz3g8gGNAWcAT+BuomlHvp2c82Vtm/sX8V1iYEBAgxR54QNauXSsdO3aUy5cvZ9r73Y2jR4/K9u3brXjmzJkO1aJff/116d27txU3btxYmjVrZsWNGjWSRx55xIrr1q0r7du3t+LmzZtLv379rLhHjx7y4YcfWvGgQYPk+++/t+KRI0daRTlFRH766Ser+rSIrfDm/v37rXjcuHFy8uRJEbGdCe7bt88qopmUlCTnzp2TqKgoK46JibGqPyvnIqdUpwZmAee4NhS6r729PXAYOAoMzcj31MSTvWX2pZoHR40SQP4zaJBUqVJFLl68mKnvl9liY2OtasgitorOydWWRUTmzJkjv/32mxV/8MEH8vXXX1txr169ZMSIEVbcuHFjGTx4sBUXL15cBgwYYMV+fn7y+uuvW3HevHnlv//9rxUbY+S9994TEdu0AIAMHz5cRGzTEgAyatQoEREJDQ0VQMaOHSsiIufPnxc3NzerevWpU6ckX758Mm3aNBEROXbsmBQrVsyqVn3o0CEpXbq0LF26VEREDhw4IJUrV5a1a9da34u6detaly7/+usvady4sezatUtERHbs2CHNmjWzvl9//vmntGrVSgIDA0VEZNOmTfLoo4/K0aNHRURkw4YN0r59ezl16pSIiKxdu1Y6deok586dExGRVatWSZcuXSQ4OFhERFasWCFPPvmk9cfN0qVLpWfPnlbiXbx4sTz77LPWtAwLFy6U3r17W5dhFyxY4PBHwi+//CIDBw604tmzZzscixkzZsjbb79txdOmTZN3333XiqdOnSoffPCBFU+ePFk+/vhjK05v4nGZwQUi0vMG7cuAZVncHaX4+uWXqT91Kid8fdmzZw958uQhKSmJoKAgSpUq5ezu3TFPT088PT2tuHr16g7Lu3fv7hAPGzbMIZ42bZpD/McffzjE586dS/4jEoCDBw86FFndtm2bNS24iDB69Gg6d+4M2C5Bzp0716oWnSdPHr755hsaNLCNJfLy8mLEiBE89NBDgG0o+7vvvmtVWPb19aVfv37WNOS+vr507drVOk4+Pj40a9bMev+8efNStWpVqxq2h4cH99xzj1Ut2hiDp6enw8COlJ8tMTGRmJgYqy0uLo7Lly+TmGgbTBsTE8OFCxesOCoqilOnTllxeHg4//77rzVqMiQkhH379lnLz58/z86dO0lKSgIgKCiIzZs3W+93/Phxh2nTDx8+zIoVK6w4MDDQYZLBvXv3snz5cr788ksA/v77b5YvX85nn30GwI4dO9i0aROffPIJAH/++Se7d++2fgY2b97MoUOHGDp0KBlBq1NrdepsKyuqG7fevZvAmBiON2yIp5sbw4cPZ8yYMfz999/ZMvm4Eq1OnX1pdWqlMtGbpUtzNi6ODxcu5OrVq/Tq1Yu33nqLkiVL3npjpVSaNPEodROPFi5M2aAgRjzxBBMmTKBMmTIMHToUYwynT5/m+eefJyhIqzMpdSc08Sh1E27GMKRlS6hXj49HjuTKlSvWss2bN7Ns2TLrOvzVq1ed1U2lshVNPErdwrP33EPBl18mPCyMMWPGWO1PPfUUp06dsu71PP/889bNcqXUjd008ZicWqtDqTvg7e7OoJYtoVkzRo8Zw6VLl64t8/a2vm7ZsiWNGze24jfeeIOVK1dmaV+Vyg5udcazwhhTMCs6opQrGxAQQJ4+fcDbmwMHDqS5Tv/+/Xn77bcB2/DYX375xSoQevXqVd5++23279+fZX1WylXdKvGUB3YbY2pkRWeUclX3eHryTP36yKxZ1Hr44Vuu7+/vz+nTpxkwYAAA+/btY/z48Zw8eRKAQ4cO8dxzz3H48GHA8RkRpXK6WyWeusABYKsx5pks6I9SLuuVgACuurvz45kzzJ0795bJIvkhRIB69eoRHh5O69atAThx4gQrV67E3d0dgF9++YVSpUpx/PhxAI4cOcLKlSuJjY3NxE+klHPcNPGIyBUR6QCMBqYZY8bZJ2dTKtep6+dHfT8/Pp85kx49ejB58uQ72j5v3rxWInr00Uc5f/485cqVA2zTKTRt2pT77EVJ586dS9u2ba0n27///nuaN29uxbt27WLRokV6pqSypdsa1SYiHwCdsE1LsM4YUzFTe6WUi/rPffdxpkEDajdtymuvvcbff/991/tKWam6adOmzJgxwyrZ8vLLL7N582Z8fX0BW0kXHx8fPDxsVa6mTJlCnz59rO1feeUVq9wM2BLVe++9Z8V//vknGzdem6LqypUrOvxbOc2dDKdeC3wHNAEOGmOuGGO2GGO+M8a8YoxpnhkdVMqV9ChWjEKenpT46CMKFSpEly5duHDhQoa/j7+/v1WXDKBPnz4sW3atZOFHH33kkEgefvhhnnzySSves2ePwyRrI0aM4PXXX7finj170qxZM4f9v/TSS1acXBoo2c8//8ySJUus+M8//7QGToCtTlvKKb71TEzd1K2qiAIVgS+xTUsdC8wEngdGAIuB40ASkJieaqXOeGl16uzNWROJvfnvv+Kxfr38tnGj+Pj4yEMPPWRVCXZVx48fd5iSYMGCBTJr1iwrfvfddx2qE3fs2FF69eplxdWqVZMuXbpY8f333y89evSw4rJlyzqsHxAQ4FAtuWbNmlYlahGRVq1ayauvvmrF3bt3lxkzZlhx3759rSkW4uPj5e2335Z169aJiK3K9ogRI2Tbtm0iInL16lX59ttvZd++fSIiEh0dLXPmzLEqRUdHR8vq1autytDR0dGyadMmq9p4VFSU7NixQ0JDQ6147969ViXv6OhoOXz4sFUp+urVq3Ly5Em5evWqiIjExMTI+fPnJS4uzupfaGioxMfHW3FYWJg1pUNsbKxcvnzZIb5y5Yo1kV9acXh4uPUzFhsbKxEREQ5xZGSkQ5zcVxGRuLg4iY6OdoiT+54cx8TEWHF8fLxVBTs5Tv5sIiIJCQmZOy0CsBpIBC4BnwD33WA9P6BRejrijJcmnuzNWYnn36go4fff5cPjx2XlypWyadMmp/QjK0VERDjMR7R9+3bZu3evFc+aNctKDCIin332mSxcuNCKBw4cKNOnT7fiDh06yBtvvGHFtWrVsqY8SExMlICAAPniiy9ExPaL39vbWz777DMREbl8+bIAMmbMGBERuXjxogAyYcIEERE5ffq0ANZcQUePHhVAfvrpJxERCQwMFMBKvHv27BFA5s+fb302QJYsWSIiIn/88YcAsmrVKhGxTXEAyIYNG0REZPny5QLI1q1bRURk0aJFAlhTKsydO1cAKzHOmDFDADl8+LCIiPzwww8CWHMTffPNNwJYiXLcuHECSEhIiPW9BSQyMlJERD7++GMBrOTw3nvviTHG+t4OHjxY8ubNa8VvvPGG+Pn5WfGAAQOkSJEiVvziiy/KfffdZ8W9evWSsmXLWnGPHj0yPfH8A7wIeKXnTVz1pYkne3Pm1Mlt9uyRgM2bJT7FdNNTpkyx5l9Rt3a3xy95UrjkX7QJCQly9uxZiYiIEBHbX/D79u2zflFHR0fLxo0brV/kERERsnLlSiu+cuWKLF68WM6ePSsiIiEhIfLLL79Y8YULF2TGjBlWfObMGfnhhx+s7U+ePClff/21nD9/XkRsiW7s2LFy4cIFEbHNBfTll19ac+8cOHBAxowZY51h7d27V7744gu5cuWKiNgS4RdffGElll27dsno0aOts5Zt27bJ559/bp2VbNmyRT777DPrDOqPP/6w5jFK/j5//vnnVrxmzRorqYvY5gL68ssvrXjp0qUOEwouWrTISuoiIvPnz0934tFpEXRahGzLmWX1FwUH8/i+fSysVo3ORYpw6dIlKlasSP78+Vm2bBnVqlVzSr+yE50WIfvSaRGUcoLHChemuKcnU8+dA6Bo0aKsX7+e+Ph46taty/vvv6+jxpS6AU08St0FDzc3nrvnHpaEhHDe/pBnrVq1+Ouvv+jevTvDhw/nwQcftCpXK6Wu0cSj1F16oXhxEoHpKYZTFy9enOnTp7NhwwYGDx6Mm5sbSUlJfPDBBxw8eNB5nVXKhWjiUeouPeDry0P58zPl/HlS3ytt2rQpzz33HGCr0/bpp59SuXJlWrZsybx584iPj3dGl5VyCZp4lEqHvvfey8HoaP4MD7/hOjVq1OD06dOMHDmSY8eO8eSTT1K6dGl0YIvKrXJU4jHGlDLGLDTGTDHGvOPs/qic78miRfF1c2PK+fM3Xa9YsWK88847HD16lN9++41WrVpRvnx5AIYNG0aXLl0YOXIka9eudZjlVKmcyGUSjz1ZXDTG7EvV3tYYc8gYc+Q2kkl1YJ6I9AFqZ1pnlbLz8/Cge7FizL54kajExFuu7+7uTocOHZg+fbpVdy0pKYn9+/fz7rvv0rp1awoWLEjbtm2tbZYuXcqqVas4dOiQjpRTOYKHszuQwo/ABGBacoO9EvZE4BEgCNhhjFkMuAMjU23fB/gTmGeM6QNMz4I+K0Wf4sWZev48v1y8SO97773j7T/66CM++ugjQkND2blzJzt27LCKfwIMHDjQmscHbGdPPXr0YPz48QB89913JCYm4uPjg6+vLz4+PpQrV47KlSsDcPz4cXx8fKxX8lQMSjmLyyQeEdlojCmTqrk+cEREjgEYY2YDnUVkJNAh9T6MMW8BH9j3NQ+YmsndVoqHCxSggrc30y5cuKvEk6xw4cK0adOGNm3aOLRv3LiREydOcPLkSeuVfJkuKSmJV1999brBCv379+ebb74hMTHRmnohWd68eXn77bcZPnw4kZGRNG3aFF9fXytp+fr60r17dzp27Eh4eDjffPON1Z78ql69OqVLlyYmJoZ///0XT09P8uTJY/1boEABvLy8rEEXKROpUi6TeG4gADidIg4CGtxk/RXAMGPM08CJtFYwxvQD+sG1h/5U9hQZGekyx+9hYNrVq/yyfj1FM+k9SpYsScmSJWncuDFge/JfRFi4cCExMTHExsYSExNDTEwM+fPnZ/369SQmJjJ48ODrlvv4+LB+/XoiIyPx8vIiKiqK4OBga3nhwoXx8/PjzJkzvPPO9Ve4Bw0aRJcuXTh69CgvvvjidcsHDx5M27Zt2bt3L4MGDcLNzY08efLg4eGBh4cHgwcPpnr16nz11Vd88cUXDsvy5MnDyy+/TKVKlVixYgUzZ850WObh4cGwYcMoVKgQ69evZ/369dYyd3d3a3svLy927NhBYGCgtczNzQ03Nze6dOmCMYa9e/cSFBRktSdv36RJE8A2IV9ISIjDck9PT6pUqQLAmTNniI6Ovm558eLFAdv0E4mJidby5HW8vb0BiImJcUjOxhjrewVY8y+lTNzJ62Rnrp547oiI7AO63WKd74HvwVYyR0t2ZF+uVHIlIDqan7Zv52S5cjxZqpSzu+OgVatWN13eocN1Fw8sIkKXLl2Ijo4mKiqKqKgooqOjKVGiBMWLF6dWrVoUKFCA+Ph44uPjiYuLIz4+npYtW1KpUiXKli1LSEiIw7L4+Hjat29PWFgYDz/8MBs2bLCWxcXFERcXx4MPPkjdunWJjo7mxIkTVnvyOk2bNsXf35/AwEAuXbpkLU9+TZ8+nXz58rF48WKmTr3+wse4ceMwxjBz5kwmTZrksMzX19eay+j7779n1qxZDsuLFy/OOXvFio4dOzpMFwFQoUIF/v33XwCaN2/Ohg0bHJbXqVOHXbt2AVC3bl12797tsLxZs2bWH1T333+/ta+Ux+u3334DoESJEpw9e9ZhbqennnqKGTNmALY/riMiIqzlxhj69OnDV199BUCBAgUAHJYPHDjQOhtOPltOufzNN9+87vt5p1yqVpv9UtsSEalmjxsBw0TkUXs8BMB+qS3dtFZb9uZKiQeg/q5dxIvw14N3XcIqV8mK4yciJCYmEhcXR0JCAomJiSQmJlKkSBEAQkNDCQ8Pt9oTEhIQEWtSvSNHjhAcHGwtT0xMxMPDwzoj2rZtG+fOnXNY7uvrS6dOnQBYvHgxZ86csZYlJSVRrFgxnnnmGQCmT5/OeftzYMmvkiVLWsu/+eYbQkJCrLMiEaFixYr07NkTgM8//5zw8HCH5TVq1KBHjx4AvP/++8TGxjrsv0GDBtbcTW+88cZ1BTxbtGhB165diYmJ4f/+7/+uW96hQwc6d+6crlptrp54PIDDQCvgDLADeFpE9mfE+2niyd5cLfGMCwri9SNH2F+vHlXsM4eqG3O146duX44pEmqMmQVsBSoZY4KMMX1FJAF4BVgJBAJzMyrpKJXRehQtihvwcybMSKpUTuIy93hEpOcN2pcBy9JappQrKZ43L60LFWLmxYt8XLasjuRS6gZc5oxHqZzgmXvu4URMDFtuUkJHqdzOZc54XEV8fDxBQUHExMQ4uyvqFgoUKEBgYGCm7d/Ly4sSJUpYQ1tvR5ciRejv5sbPFy7wsH3EkFLKkSaeVIKCgvDz86NMmTJ6qcTFRURE4Ofnlyn7FhFCQkIICgqibNmyt72dn4cHnfz9mXvxIuMqVCBPNn/eQqnMoP8rUomJicHf31+TTi5njMHf3/+uznx7FCtGSEICGy5fzviOKZUDaOJJgyYdBXf/c9C2cGHyubvzy6VLGdwjpXIGTTxKZTBvd3c6+PuzIDiYBJ36WqnraOJRKhN0K1qU4Ph4NurcOkpdRxOPCwoMDKRMmTIk2f9aTkxMpE2bNkybNu0WWypX0a5wYXzc3PRym1Jp0MTjgipXrkzlypWt4oNDhw6lUqVKPPfcc07umbpdPu7uPObvz4JLl0h0obJUSrkCHU7tot544w2+/PJL4uPj2bx5M+vWrXN2l9QderJoUX65dIk/Ll+meaFCzu6OUi5DE89NvP7vv+yJjMzQfdbKl4+xFSvecr02bdrw5ptvMmTIEDZu3HhHDzEq19De3x9v++U2TTxKXaOJx4U99NBD1K5d25pUSmUvvu7utC9cmPmXLjG+YkXcdZi+UoAmnpu6nTOTzHTgwAFeeOEFKz558iRffPEFIkL58uVp27Ytw4YN4/7772f79u2sWLHCib1VaXmyWDHmBwez+coVmhYs6OzuKOUSNPG4sP3791OtWjUr/vrrr/H29sbb25u9e/dy+vRpRo0aRUBAAI8++qgTe6pupH3hwngaw6LgYE08Stlp4nFRp0+fpmDBguTLl89qS0pKolevXtSoUQPAmoI25bS3yrX4eXjQslAhFgUHM7p8eT1OSqGJx2WVLFmSY8eOObS98sorvPvuu9x77734+fnx4osv8s4771C+fHk8PPRQuqpO/v4MCA0lMDpaZyZVCk082Urp0qX5+eefHdpmzZoFQLdu3ZzRJXUbOvr7M+Dff1kcHKyJRyn0AdIcY968ec7ugrqBEl5e1M2Xj0UhIc7uilIuQROPUlmgc5EibAsP53xsrLO7opTTaeJRKgt0LlIEAZboWY9SmniUygrVfX0pnTcvizXxKJV9E48xppwx5gdjzLwUbb7GmJ+MMZOMMc84s39KpWSMoXORIqwOCyMqMdHZ3VHKqZySeIwxU4wxF40x+1K1tzXGHDLGHDHGvHOzfYjIMRHpm6q5KzBPRF4COmVwt5VKl05FihCTlMSasDBnd0Upp3LWGc+PQNuUDcYYd2Ai0A6oAvQ0xlQxxlQ3xixJ9Sp2g/2WAE7bv9Y/K5VLaVqgAAU9PFgUHOzsrijlVE55jkdENhpjyqRqrg8cEZFjAMaY2UBnERkJdLjNXQdhSz57yMaXEVXOlMfNjTaFCrEiNBQR0SoGKtdypQdIA7h2tgK2JNLgRisbY/yBT4Daxpgh9gS1AJhgjHkM+O0G2/UD+gEULVqU9evXOywvUKAAERER6fgYKqskJiZm+rGKiYm57mckPcoCc4EfNmygQobtNXuKjIzM0O+tyj5cKfHcEREJAfqnaosCXkh7C2ud74HvASpVqiTNmzd3WB4YGIifn1+G9lVljoiIiEw/Vl5eXtSuXTvD9vdAbCyjtm7lUtmyvFi6dIbtNztav349qf//qdzBlS5HnQFKpohL2NtUNnbs2DH69u1rlfQ5deoUjz/+OH369OHTTz91at8WLlzISy+9RI8ePVi1alWWvGfxvHmpnS8fy0NDs+T9lHJFrpR4dgAVjTFljTGewFPAYif3yalatmxJQkLCHW8XFhZGly5d0vXeffr0oVixYg7TMgCsWLGCSpUqUaFChdtKHOXKleOHH36w4r1799KtWzemTJnCX3/9leY2Fy5c4Omnn6ZcuXLUrVuXRo0a8euvv970fVq0aMHKlSsd2saOHct//vOfG27z+OOPM2nSJL799lvmzJlzy8+SUdoWLsyWK1e4HB+fZe+plCtx1nDqWcBWoJIxJsgY01dEEoBXgJVAIDBXRPY7o3+uYP/+/fj7+99V1elChQoRGhpKSDoeVuzdu/d1E8slJiYycOBAli9fzoEDB5g1axYHDhwAbAmlQ4cODq+LFy9et9+GDRvyww8/0LJlS9q2bXvdchHh8ccfp2nTphw7doxdu3Yxe/ZsgoKCbtrfnj17Mnv2bIe22bNn07Nnz1v27eOPP2bgwIG3/b1Jr3aFC5MIOqxa5V4ikmtf999/v6R24MCB69qc4ZNPPpEZM2aIiMiePXukSZMmUrlyZTHGCCDvvfee7Ny5U5o1a2Zts3fvXmnUqJGIiIwaNUqmTp2arj4cP35cqlatasVbtmyRNm3aWPGIESNkxIgRt7WvJ554QkREPv/8c9mwYYNDW0pr1qyRpk2b3nRf06dPl3r16kn16tWlX79+kpCQICEhIVK0aFGJjY21+l6yZElJSkq64X6SkpLkv//9r6xevfqG62TGz0N8YqIU2LhR+gQGZvi+s5Pff//d2V1QdwnYKen43etKl9pUCsuWLeOxxx4jJiaGHj16MHr0aA4cOMDQoUN56623+PDDD6lcuTKHDx+2tnn//ff56KOPAOjcuTMLFy68br9NmjShVq1a173WrFlzyz6dOXOGkiWv3YYrUaIEZ87c/DZcSEgI/fv356+//mLkyJG0bduW8ePH079/f8qUKXPd+vv376dOnTo33F9gYCBz5sxh8+bNbN68GXd3d37++WcKFy5M/fr1Wb58OWA72+nevftNhyx/9dVXrFmzhnnz5vHtt9/e4tNnHA83Nx4pXNgaVq1UbpNtR7VlmbRG3XTvDgMGQHQ0tG9//fLevW2v4GBIPU/ObQwfjY6OJi4ujoIFC7JkyRLq1KlD/fr1AahRowYrVqzAGIOPjw/e3t5cvnyZY8eOERYWRuvWrQGoVKkShw4dum7ff/zxxy3fPyP5+/tf90v9TqZwGDhwIJs2bcLT05MdO3awdu1adu3aRb169UhKSiI2NpZixWzPEydfbuvcuTOzZ892uLeUlkGDBjFo0KA7/1AZoF3hwsy7dIl/oqKomWKWWaVyA008LsjHxwdjDJGRkezbt4/q1atby3bv3u1wRlClShUOHjzI8OHD+fjjj632kydPUrZs2ev23aRJkzSffRk9erSVtG4kICCA06evPWoVFBREQEDAHX22W6latSrz58+34okTJxIcHMyDDz4I2C4NP//884wcOfK64dSdO3fmjTfeYPfu3URHR1O3bt0M7VtGalu4MADLQ0I08ahcRxPPrdzsDMXH5+bLixS5rTOctDz66KOsWLECf39/1q1bB8Dhw4dZsGABW7ZssdarWrUqU6ZMQUR4+OGHrfZFixbRuXPn6/abnjOeevXq8e+//3L8+HECAgKYPXs2M2fOvOv9paVly5a8++67fPPNN9aItOjoaGt5q1atrATj7e1NaGgoERERlC5dmnz58tGiRQv69OlDz549M7RfGe2+vHmp6evL8tBQ3snlz/OoXCg9N4iy+8uVBxfs3LlTnnnmGYmIiJAOHTpI1apVpVmzZrJr1y6H9aZNmybGmOvaW7duLefOnbvr93/qqaekePHi4uHhIQEBATJ58mQREVm6dKlUrFhRypUrJx9//PFd7/9mzp49Kz169JAyZcpIvXr1pHnz5jJ79mxr+ezZs6VmzZpStWpVqVOnjmzdutVa9uuvvwoggRl04z4zfx7eOXpU3H//Xa7Ex2fae7gyHVyQfZHOwQVGcvHNzUqVKknq+yCBgYFUrlzZST1yVKNGDXbv3n3HQ6rDwsLo3LkzGzduzKSeuYasqFyQmT8P68LCaPX33/xWrRodihTJlPdwZVq5IPsyxuwSkQfvdnsd1ebC/vnnn7t+jienJ52c4KH8+fFyc9PneVSuo4lHKSfxcnenSYECmnhUrpO7L7WZSvId3zm0FVhegApFcnvdYJXsSPARrrS74uxuKOVSWtBCL7UppZTKPnL3cOr7ofmh5g5NgYGB+FXWaRGygyyZFiHQi9qScdMipJYoQrHNm+lUpAhTH3gg097HFenggmwsnXMY6hmPUk7kbgwtCxViTViYls9RuYYmHqWcrHWhQgTFxnL46lVnd0WpLKGJRykna12oEKDTJKjcQxOPUk5WzsuLMl5erNXEo3IJTTxKOZkxhtaFCrEuLIxEvc+jcgFNPCpTHTt2jL59+9LNPj3EwoULeemll+jRowerVq1yat9cqS+tCxXiSmIiO9OoHK5UTqOJx4W1bNmShISEO94uLCyMLl26pOu9+/TpQ7FixahWrZpD+4oVK6hUqRIVKlTg008/veV+ypUr5zAvzuOPP86kSZP49ttvmTNnTprbXLhwgaeffppy5cpRt25dGjVqxK+//nrT92nRogUrV650aBs7dqxV4Tott9OXrNKiYEEA1l++7NR+KJUVNPG4qP379+Pv73/XtdpCQ0MJCQm56/fv3bs3K1ascGhLTExk4MCBLF++nAMHDjBr1iwOHDgAwN69e+nQoYPD6+LFizfc/8cff8zAgQOvaxcRHn/8cZo2bcqxY8fYtWsXs2fPJigo6Kb9TZ4ELqXZs2fTs2fPW/btRn3JSsU8Pans48MGTTwqN0hPaevs/nLlaRE++eQTmTFjhoiI7NmzR5o0aSKVK1cWY4wA8t5778nOnTulWbNm1jZ79+6VRo0aiYjIqFGjZOrUqenqw/Hjx6Vq1apWvGXLFmnTpo0VjxgxQkaMGHFb+3riiSdERCQpKUn++9//yurVq9Ncb82aNdK0adOb7mv69OlSr149qV69uvTr108SEhIkJCREihYtKrGxsVbfS5YsKUlJSTfcz636IpK1Pw/9Dx0Sv40bJT4xMcve05l0WoTsi3ROi6BnPC5q2bJlPPbYY8TExNCjRw9Gjx7NgQMHGDp0KG+99RYffvghlStX5vDhw9Y277//Ph999BFgm41z4cKF1+23SZMm1KpV67rXmjVrbtmnM2fOULJkSSsuUaIEZ86cuek2ISEh9O/fn7/++ouRI0fy1VdfsWbNGubNm3fdlNhgO9NLOcNqaoGBgcyZM4fNmzezefNm3N3d+fnnnylcuDD169dn+fLlgO1sp3v37hhz40esb9WXrNasQAEiEhPZExnp7K4olamydckcY0w5YChQQES62dseBx4D8gM/iEi67hqnVdKje/fuDBgwgOjoaNq3b3/d8t69e9O7d2+Cg4Otm+rJ1t/GjKTR0dHExcVRsGBBlixZQp06dahfvz5gm6NnxYoVGGPw8fHB29uby5cvc+zYMcLCwqzpqytVqkTquYYgfTOQ3g1/f//rfqkPGjTotrcfOHAgmzZtwtPTkx07drB27Vp27dpFvXr1SEpKIjY2lmLFigHXLrd17tyZ2bNnO9xbSsugQYPuqC+ZrZn9Ps+GK1d4MH9+53ZGqUzktMRjjJkCdAAuiki1FO1tgXGAOzBZRG54B1tEjgF9jTHzUrQtBBYaYwoBowHnDle6Cz4+PhhjiIyMZN++fVSvXt1atnv3boczgipVqnDw4EGGDx/Oxx9/bLWfPHmSsmXLXrfvJk2aEJHGyKnRo0dbSetGAgICOH36tBUHBQUREBBwR5/tVqpWrcr8+fOteOLEiQQHB/Pgg7ZCuCLC888/z8iRI6+r1ZY8Jfbu3buJjo6mbt26Gdq3zHZv3rxU9PZmw+XLvJnizFKpnMaZZzw/AhOAackNxhh3YCLwCBAE7DDGLMaWhEam2r6PiNz47jX8z76vdLnZGYqPj89NlxcpUuS2znDS8uijj7JixQr8/f1Zt24dAIcPH2bBggVs2bLFWq9q1apMmTIFEeHhhx+22hctWkTnzp2v2296znjq1avHv//+y/HjxwkICGD27NnMnDnzrveXlpYtW/Luu+/yzTffWCPSoqOjreWtWrWyEoy3tzehoaFERERQunRp8uXLR4sWLejTpw89e/bM0H5llWYFC/LLxYskiuB+k8uESmVr6blBlN4XUAbYlyJuBKxMEQ8BhtzGfual+NoAo4DWt9rOlQcX7Ny5U5555hmJiIiQDh06SNWqVaVZs2aya9cuh/WmTZsmxpjr2lu3bi3nzp276/d/6qmnpHjx4uLh4SEBAQEyefJkERFZunSpVKxYUcqVKycff/zxXe//Zs6ePSs9evSQMmXKSL169aR58+Yye/Zsa/ns2bOlZs2aUrVqValTp45s3brVWvbrr78KIIGBgRnSl6z+eZh+7pzw+++yOzw8S9/XGXRwQfZFOgcXOHUiOGNMGWCJ2C+1GWO6AW1F5EV73AtoICKv3GB7f+ATbGdIk0VkpDFmEPA8sAPYIyLfptqmH9APoGjRonXnzp3rsM8CBQpQoYJrTATXqFEj/vjjjzseUh0WFkbPnj2vGw6d0yQmJuLu7p6p73HkyBGuXMm6ieAuAj2AgUC3W6yb3UVGRpIvXz5nd0PdhRYt0jcRXLZOPOlVqVIlSX0DPjAwkMqVK2fG26kMlhXz8Tjj56H8n39SI18+fk318G5Oo/PxZF/GmBw1A+kZIOVd1RL2NqVyjWYFC7Lx8mWSnPhHoVKZydUSzw6gojGmrDHGE3gKWOzkPimVpZoVLEhoQgL7o6Kc3RWlMoXTEo8xZhawFahkjAkyxvQVkQTgFWAlEAjMFZH9zuqjUs7QrEABAC2fo3Ispw2nFpE0x7uKyDJgWRZ3RymXUcbbm1J587LxyhVeKVHC2d1RKsO52qU2pRTQuEABNl+5gjMH/yiVWTTxKOWCHi5QgLNxcZyMiXF2V5TKcJp4lHJBD9lrtW0OD3dyT5TKeJp4lHJB1fPlw8/dnS1Z+PCqUllFE49SLsjdGBrmz89mTTwqB9LE44ICAwMpU6YMSUlJgK00TJs2bZg2bdottlQ5yUP587M3Korwu5j+XClXponHBVWuXJnKlSuzZMkSAIYOHUqlSpV47rnnnNwzlZUeLlCAJGCb3udROUy2ngguJ3vjjTf48ssviY+PZ/PmzdbUCCr3aJA/P27A5itXeKRwYWd3R6kMo4nnFtab9Rm6v+bS/LbWa9OmDW+++SZDhgxh48aN5MmTJ0P7oVxffg8Pqvv66sg2leNo4nFhDz30ELVr16Z48eIAREVFMWDAADw9PWnevDnPPPOMk3uoMtvDBQow7cIFnRhO5SiaeG7hds9QMsOBAwd44YUXrHjBggV069aNjh070qNHD008ucBDBQrw9dmz7I2MpFYmTwGhVFbRwQUubP/+/VRLMSdLUFAQJUvaZo3I7AnQlGt4WB8kVTmQJh4Xdfr0aQoWLOgwQ2OJEiUICgoCsIZaq5yttJcX93l66oOkKkfRS20uqmTJkhw7dsyhrWvXrrzyyissXbqUjh07OqlnKisZY3jIXjBUqZxCE0824uvry9SpU53dDZXFHsqfn3mXLnE2Npb78uZ1dneUSje91KaUi2tgv8+zXe/zqBxCE49SLq52vnx4GMP2iAhnd0WpDKGJRykX5+3uTk1fXy2do3IMTTxp0FkfFbjWz0H9/PnZERFBkgv1Sam7pYknFS8vL0JCQlzql47KeiJCSEgIXl5ezu4KYLvPE5GYyMHoaGd3Ral001FtqSQ/K3Pp0iVnd0XdQkxMTKYmBi8vL0qUKJFp+78T9e1VC7aHh1PF19fJvVEqfTTxpJInTx7Kli3r7G6o27B+/Xpq167t7G5kiUo+PhRwd2dbRAS9773X2d1RKl2y7aU2Y0w5Y8wPxph5qdp9jTE7jTEdnNU3pTKamzHUy59fh1SrHMEpiccYM8UYc9EYsy9Ve1tjzCFjzBFjzDs324eIHBORvmksGgzMzcj+KuUK6vv58U9UFFcTE53dFaXSxVlnPD8CbVM2GGPcgYlAO6AK0NMYU8UYU90YsyTVq1haOzXGPAIcAC5mbveVynoN8ucnQYTdkZHO7opS6eKUezwistEYUyZVc33giIgcAzDGzAY6i8hI4HYvmzUHfLElrqvGmGUi4lBN0xjTD+hnD2NTn3U5SQEgI4tx3e3+bne721nvVuvcaPmdtBcBgm/Rj6yQkcfvlvtqnP4+pPf43emyG63vCscvu/3fu511s+L/XqWbdfCWRMQpL6AMsC9F3A2YnCLuBUy4yfb+wLfAUWBIqmW9gQ630Yedzvr8qfrxvSvs73a3u531brXOjZbfSXtOPH6Zfewy4vjd6bKbHFOnH7/s9n/vdtbNDv/3su2oNhEJAfrfYNmPWdubdPvNRfZ3u9vdznq3WudGy++03RVkZN8y+9jd7ro3W+dOl+WWY5ee/WXk8XP5/3vGnr2ynP1S2xIRqWaPGwHDRORRezwEQGyX2jKrDztF5MHM2r/KXHr8sjc9ftlXeo+dKw2n3gFUNMaUNcZ4Ak8BizP5Pb/P5P2rzKXHL3vT45d9pevYOeWMxxgzC9tAgCLABeADEfnBGNMeGAu4A1NE5JMs75xSSqlM5bRLbUoppXInV7rUppRSKhfQxKOUUipLaeJJwV7n7SdjzCRjzDPO7o+6Mzeq36dcnzHmcfv/uznGmDbO7o+6M8aYysaYb40x84wx/7nV+jk+8dxhXbiuwDwReQnolOWdVde5k+MnN67fp5zgDo/dQvv/u/5AD2f0Vzm6w+MXKCL9ge7Aw7fad45PPNxBXTigBHDavppWYnQNP3L7x0+5lh+582P3P/ty5Xw/cgfHzxjTCVgKLLvVjnN84hGRjUBoqmarLpyIxAGzgc5AELbkA7nge5Md3OHxUy7kTo6dsRkFLBeR3VndV3W9O/2/JyKLRaQdcMvbFLn1l2sA185swJZwAoAFwBPGmG9w7TIfuV2ax88Y42+M+RaonVz5QrmcG/3fexVoDXQzxqRZCku5hBv932tujBlvjPmO2zjjyba12jKDiEQBLzi7H+ru3Kx+n3JtIjIeGO/sfqi7IyLrgfW3u35uPeM5A5RMEZewt6nsQY9f9qXHLnvLkOOXWxOPM+rCqYyjxy/70mOXvWXI8cvxicdeF24rUMkYE2SM6SsiCcArwEogEJgrIvud2U+VNj1+2Zceu+wtM4+f1mpTSimVpXL8GY9SSinXoolHKaVUltLEo5RSKktp4lFKKZWlNPEopZTKUpp4lFJKZSlNPEoppbKUJh6llFJZShOPUjmUMaa7Mea8MaaBMUZu8Uoyxvg5u88qd9Dq1ErlXI9hK1H/L9AoRfsgoAOQcorpeBGJyMK+qVxMS+YolQMZY9yA88B/RGR+qmVLgUIi8pBTOqdyPb3UptRtMsY0s1+Wap+irax9XvoMnUvGGONmjIk0xrxujBljf48wY8zb9uW9jDEH7OssMMZ4p9pFPaAgsDqN3VcH9mVkf5W6E5p4lLpNIrIB+B0YCmCMKQAsAbYDb6Re3z6ds8etXjd4u3KAL/A6EAs8jW0++8+MMROwlaN/G3gHeJzrJzB8DPhDRMJT9akAtvlU9t7xN0CpDKL3eJS6Mx8AG40xbYA3gXjgKRFJTGPd54Gpt7FPk0Zbdfu/Y+yzc2KM+RfbfPaVgdZiv05ujOkHVEq1/WPAz2nst5r9Xz3jUU6jiUepOyAifxhj1gC/ApeBBiISeYPVf8N2yetu1LDv/5sUbb72fz8Vx5uzvkBocmCMuReoje0sKbXkxKNnPMppNPEodeeOAK2B10Qk6CbrhQJX7vI9qgObRCQ+RVsNIAHYmNxgjPEByuB4BtMeOCYih26w3wsiEnyX/VIq3fQej1J3wH5Zqw/wN9D3Fqs/j+1S3K1eaakO7EnVVhM4KCKxqdZzA/5J0fYYtvtBaamGnu0oJ9MzHqVukzHmEWAitoRzGNhqjGknIstvsMldXWqzj1CrgC25pVTjBm1RwFH7tp7AI8ATN9h9NWD6nfZJqYykiUep22CMqQr8AowSkWn2tjXAMCDNxCMiIUDIXbxdVWxnMWklmdTDtmsA+0UkyR43xTZYYUMan+FewB8941FOppfalLoFY0wxbMOmVwHvpVg0HKhvjHksg9+yOinOYux9KASUwPGSGtgST+rLbGtSXY5LuV/QEW3KybRygVI5iDHmMPC5iExydl+UuhFNPEoppbKUXmpTSimVpTTxKKWUylKaeJRSSmUpTTxKKaWylCYepZRSWUoTj1JKqSyliUcppVSW0sSjlFIqS/0/mZT5XKYFx6UAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xx=np.logspace(np.log10(1),np.log10(1000),1000)\n",
    "plt.loglog(xx,Yeq(xx),'c',label=r'$Y_{\\rm eq}$')\n",
    "\n",
    "k=fk(σv=1E-18)\n",
    "x=np.logspace(np.log10(1),np.log10(1000),10000)\n",
    "Y=integrate.odeint(BE,Yeq(x[0]),x)\n",
    "plt.loglog(x,Y,'r--',label=r'$\\langle \\sigma v \\rangle = 10^{-18}\\ {\\rm GeV}^{-2}$')\n",
    "k=fk(σv=1E-17)\n",
    "x=np.logspace(np.log10(1),np.log10(1000),10000)\n",
    "Y=integrate.odeint(BE,Yeq(x[0]),x)\n",
    "plt.loglog(x,Y,'k:')\n",
    "k=fk(σv=1E-16)\n",
    "x=np.logspace(np.log10(1),np.log10(1000),10000)\n",
    "Y=integrate.odeint(BE,Yeq(x[0]),x)\n",
    "plt.loglog(x,Y,'k:')\n",
    "k=fk(σv=1E-15)\n",
    "x=np.logspace(np.log10(1),np.log10(1000),10000)\n",
    "Y=integrate.odeint(BE,Yeq(x[0]),x)\n",
    "plt.loglog(x,Y,'k:')\n",
    "k=fk(σv=1E-14)\n",
    "x=np.logspace(np.log10(1),np.log10(1000),10000)\n",
    "Y=integrate.odeint(BE,Yeq(x[0]),x)\n",
    "plt.loglog(x,Y,'k:')\n",
    "k=fk(σv=1E-13)\n",
    "x=np.logspace(np.log10(1),np.log10(1000),10000)\n",
    "Y=integrate.odeint(BE,Yeq(x[0]),x)\n",
    "plt.loglog(x,Y,'k:')\n",
    "k=fk(σv=1E-12)\n",
    "x=np.logspace(np.log10(1),np.log10(1000),10000)\n",
    "Y=integrate.odeint(BE,Yeq(x[0]),x)\n",
    "plt.loglog(x,Y,'k--',\n",
    "          label=r'$\\langle \\sigma v \\rangle = 10^{-12}\\ {\\rm GeV}^{-2}$')\n",
    "plt.hlines(1E-11,1,1000,linestyles='-',color='m',lw=2,label='$Y_0$')\n",
    "plt.ylim(1E-14,1)\n",
    "plt.xlim(1,1000)\n",
    "plt.xlabel('$x = m/T$', size= 15)\n",
    "plt.ylabel('$Y$', size= 15)\n",
    "plt.title('$m = 100$ GeV', size= 15)\n",
    "plt.legend()\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The required $Y_0$ seem to imply\n",
    "$$\n",
    "\\langle \\sigma v \\rangle = 10^{-8}(\\hbar c)^2\\text{GeV}^{-2}c\n",
    "$$\n",
    "instead of $2.7\\times 10^{-10}(\\hbar c)^2\\text{GeV}^{-2}c$"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "WIMP.ipynb",
   "provenance": [],
   "version": "0.3.2"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}