lgarrison/scale_free_sigma8.ipynb

## scale_free_sigma8.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# $\\sigma_8$ in power-law cosmologies\n",
    "Author: Lehman Garrison (<https://lgarrison.github.io>)\n",
    "\n",
    "#### Update 10/25/2019:\n",
    "Thanks to Daniel Eisenstein for deriving an even simpler expression!  The notebook has been updated with this version."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Power-law cosmologies ($P(k) = A k^n$) have an analytic form for $\\sigma^2(R)$ with a spherical tophat window $W_T$:\n",
    "\n",
    "$$\\begin{align}\n",
    "\\sigma^2(R) &= \\int \\frac{d^3k}{(2\\pi)^3}\\tilde W_T^2(kR)P(k) \\\\\n",
    "            &= A \\int_0^\\infty \\frac{k^2dk}{2\\pi^2} \\left[\\frac{3(\\sin(kR) - kR\\cos(kR))}{(kR)^3}\\right]^2 k^n\\\\\n",
    "            &= \\frac{9 A (n+1) 2^{-n-1} R^{-n-3} \\sin \\left(\\frac{\\pi n}{2}\\right) \\Gamma (n-1)}{\\pi^2 (n-3)}\n",
    "\\end{align}\n",
    "$$\n",
    "where we used the Fourier transform $\\tilde W_T(kR)$ of the spherical tophat window.  The result in the last line (obtained via Mathematica) is only valid for $-3 < n < 1$; otherwise, the variance is undefined.\n",
    "\n",
    "However, the result is numerically unstable if $n$ is an integer in this range, since the function $\\Gamma(n-1)$ in the numerator goes to infinity for 0 and the negative integers.\n",
    "\n",
    "The $\\sigma^2(R)$ function itself is perfectly smooth and finite in this range, as seen below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEICAYAAABxiqLiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hc1Z3/8fdXvVqyqnuVXLCNccGAMb0ZNqaGJEBIKAkhgZCQZLOw7IawyW568gtsIHESAiy9BAihd0wzLrhgMO6y5SLJkiWrWP38/pgREY5kz4xndKd8Xs8zj6Zc3fuZa3m+c+499xxzziEiIhKsJK8DiIhIbFIBERGRkKiAiIhISFRAREQkJCogIiISkhSvAwyUoqIiN2bMGK9jiIjElGXLlu12zhX39VrcFxAzWwAsKCsrY+nSpV7HERGJKWZW0d9rcX8Iyzn3lHPuqry8PK+jiIjElbgvICIiEhkqICIiEhIVEBERCYkKiIiIhEQFREREQqICIiIiIYn7AmJmC8xsYUNDg9dRREQG3JMrtnPPO1sisu64LyC6DkREEtmTK3bw8NJtEVl33BcQEZFEVtPYRnFOekTWrQIiIhLHqhtbKc5VARERkSB0dzt2N7WrgIiISHD2tLTT1e10CEtERIJT09QGQHFuRkTWrwIiIhKnahp7CohaICHRdSAikqhUQA6RrgMRkUTVU0BKVEBERCQYNY1tZKUlk50emclnVUBEROJUTVNbxA5fgQqIiEjciuRV6KACIiISt2oa1QIREZEQ6BCWiIgEra2zi/qWDh3CEhGR4Oxuagcidw0IqICIiMSlSF9ECCogIiJxSQUkDDSUiYgkIhWQMNBQJiKSiHoKSGG2CoiIiAShpqmVwVmppKVE7mNeBUREJA5F+iJCUAEREYlLKiAiIhKSmqY2SiI0E2EPFRARkTjjnFMLREREgtfU1klrR3dEhzEBFRARkbgzENeAgAqIiEjcUQEREZGQVKuAiIhIKD5pgegciIiIBKOmqY3UZCMvMzWi21EBERGJMzWNbRTlpJOUZBHdTtwXEI3GKyKJZiCuAYEEKCAajVdEEs3WuhaG52dGfDtxX0BERBJJa0cXFbXNlJfkRHxbKiAiInFkS20z3Q7KSnMjvq2QCoiZZZtZcrjDiIjIoVlf1QRAWXGUtEDMLMnMLjazp82sGlgL7DSzNWb2CzMrj2xMEREJxIbqJsxgXHF2xLcVaAvkVWA8cCMwxDk30jlXAhwHvAv81My+GKGMIiISoA01TYwqyCIjNfIHiVICXO5U51zH/k865+qAx4DHzCyyV6yIiMhBbahqGpDDVxBgC6Sv4tHDzC4/2DIiIhJ5nV3dbN7dTFlpFBWQg7glDOsQEZFDtLWuhfau7gFrgQR0CMvMVvX3ElAavjgiIhKqDdX+HlgDcA0IBH4OpBQ4A9iz3/MGvB3WRCIiEpINNdFZQP4O5DjnVuz/gpm9FtZEIiISkg1VTQwZlEFuxsD0aQqogDjnrjzAaxeHL46IiIRqQ00T5QN0Ah00lImISFzo7nZsqG5i/ACdQAcVEBGRuLBzbyst7V0Ddv4DVEBEROJCTw+sgRiFt0fQBcTMTu79U0REvLe+qhEYuB5YEFoL5Jf7/RQREY9trGmiIDuNwpzIz0TY41AOYUV2st0w0ZS2IpII1g/gGFg94v4ciKa0FZF419HVzZodezls2KAB3W7cFxARkXi3dmcj+zq6mDV68IBuVwVERCTGLa2oA4iJAtLk/9kYziAiIhKaZRV7GJaXwbD8zAHdbtAFxDl3fO+fIiLirWUVe5g5wK0P0CEsEZGYtr1+HzsbWpmtAiIiIsFYVuGbZWPW6IIB37YKiIhIDFu2pY6stGQmD80d8G0HVEDM7FIzqzGzSjP7sv+5o83sx2a2LLIRRUSkP8u27uGIkfmkJA98eyDQLf4AOAs4AhhrZi8CjwBpwLcjlE1ERA6gua2Tj3Y2Dnj33R6BzkjY5JxbAmBmtwBVwATnXH3EkomIyAGt3FZPV7eL+gIyxMyuAj723ypVPEREvLW0Yg9mMGNUdBeQHwKHA5cA04BcM3sJeB943zl3f2TiiYhIf5ZW7GFCSS55mQMzB/r+Ai0gC51zrueBmY3AV1CmAWcC95uZ9V5GREQip72zm+UVezj7iGGeZQi0gLxqZo8BTzrntjrnKoFKfyvkODO7G3gVuCtCOUVEpJclW+poauvkpIklnmUItIDMB64AHjCzsUA9kAEkAy8Av3HOrYhMRBER2d/LH1WTlpLEsWWFnmUIqIA451qB24HbzSwVKAL26US6iIg3XllbxdzxhWSlBdoOCL9QBlPscM7tVPEQEfHGppomttS2cMok7w5fQQAFxMxOM7M/mtkR/sdXRT5WdGnt6GJfe5fXMUREAHhlbTUAJ0V7AQG+Afwr8EUzOxnf1egJo6Glg0n/+Rz3La7wOoqICOA7/zGxNJcRg7M8zRFIAalxztU7574HnA4cGeFMUWVQZgppKUnUNLZ5HUVEhIZ9HSzZUsfJk71tfUBgBeTpnjvOuRuAeyIXJ/qYGcU56SogIhIVFq2vobPbeX7+AwIoIM65J3vum9nNzrnbIhsp+hTnplOtAiIiUeCVj6rJz0r1bPiS3oLt/3WzmWUBBcBy4EHn3J7wx4ouxbnpbK1t8TqGiCS4zq5uXltXw0kTS0hOMq/jBN2N1wGtwPPASOBtM5se9lRRpiQ3nerGVq9jiEiCW7R+N3XN7cyfOsTrKEDwLZC1zrmb/fcfNbO7gN8DJ4c1VZQpzk1nT0sH7Z3dpKVoEkcR8cZf399Oflaqp8OX9Bbsp+FuM5vV88A5tw4oDm+k6FOSmwHA7iadBxERbzS2dvDCml185vChUfNFNtgWyHXAg/5pbFfjG5F3c9hTRZni3HQAahrbGJaf6XEaEUlEz32wi7bObs6bMcLrKJ8Iqow551biu5DwAf9TrwIXhTtUtCnxFxD1xBIRrzz+/nbGFGYxc1S+11E+EfQoXM65NnzXhjx9sGUjycyy8Q3w2A685py7L1Lb6t0CEREZaDsb9vHOplq+dUo5Zt73vuoRHQfS/MzsTjOrNrMP9nt+vpl9bGYbzOwG/9PnA486574KnB3JXEU5KiAi4p0n3t+Bc3DejOFeR/mUqCog+Cakmt/7CTNLBn6Hb+bDw4CLzOwwYASwzb9YREc6TEtJYnBWqrryisiAc87x+PuVzBo9mNGF2V7H+ZSgCoj5fNHMfuB/PMrM5oQrjHPuDaBuv6fnABucc5ucc+3Ag8A5QCW+IgL9vA8zu8rMlprZ0pqamkPKVpKboRaIiAy45Vv3sK6qifNnRlfrA4JvgdwOHMM/Tpw34msdRNJw/tHSAF/hGA78FbjAzO4AnurrF51zC51zs51zs4uLD623sYYzEREv3PV2BbkZKZx7RPQVkGBPoh/lnJtpZu8DOOf2mFlaBHL11tcZI+ecawYuj/C2P1Gcm87m3c0DtTkREXY1tPLs6p1cNncM2enezTzYn2BbIB3+cxIOwMyKge6wp/q0SnzDpvQYAeyI8Db/SUmub0Re59xAb1pEEtR9iyvoco4vHTPG6yh9CraA3Ao8DpSY2X8DbwL/E/ZUn7YEKDezsf7WzheAv0V4m/+kODed9q5u9u7rHOhNi0gCau3o4v7FWzllUimjCr2dOKo/QbWJnHP3+a9CPwXfoaVznXMfhSuMmT0AnAgUmVklcLNz7s9mdi2+ARyTgTudc2uCWOcCYEFZWdkhZfvkWpCmVvKyUg9pXSIiB/P3VTupbW7n8mPHeB2lX6FcSLgWWBuBLDjn+ryq3Tn3DPBMiOt8Cnhq9uzZXz2UbD0FpHpvG2UluYeyKhGRA3LO8Ze3NlNeksPc8YVex+lXsN147zaz/F6PB5vZneGPFX16BlSs0YCKIhJh72ysZc2OvVx27JiouvJ8f8GeAzncOVff88A/mdSM8EaKTr1bICIikfTbl9dTOiidC2ZGz8CJfQm2gCSZ2SfzKJpZASEcBotFgzJSSE9JUgtERCLq3U21LN5cx9UnjCcjNdnrOAcU7If/r4B3zOwR/+MLiXwvrKhgZr6LCfdqOBMRiZzfvrSe4tx0LpozyusoBxXscO734BvEsMp/O9//XNQyswVmtrChoeGQ11WSm64WiIhEzHub63hnU21MtD4g+JPo6fjmAxkEFACf7RkXK1o5555yzl2Vl5d3yOsq9l9MKCISCbe+vJ6inHQujoHWBwR/DuRJfAMZdgLNvW4JQeNhiUikLN5Uy5sbdvO148eRmRb9rQ8I/hzICOfc/IMvFp9KcjOob+mgrbOL9JTY+AcWkejX3e348dMfMSwvg0uPGe11nIAF2wJ528ymRSRJDOjpyru7qd3jJCIST55cuZ3V2xv41/kTY+LcR49gC8g8YJl/dsBVZrbazFZFIlg0KtHUtiISZvvau/j5cx9z+Ig8zpkefUO2H0iwh7DOjEiKCArXWFjQ+2JCdeUVkfD485ub2NnQyv/7/BEkJUXvVed9CbYbbwWwFygFRve6Ra1w9sLScCYiEk7Ve1u5/bWNnH5YKUeNi94xr/oTVAvEzL4CfAvfnBwrgKOBd4CTwx8t+hTmpGGm4UxEJDx+9PRHdHY5bjxrstdRQhLsOZBvAUcCFc65k/CNg3Vok43HkNTkJIblZbKlNmF6LotIhLy+roanVu7gGyeNZ2xRttdxQhJsAWl1zrWC76JC/9DuE8MfK3qVl+awrqrJ6xgiEsP2tXfxH0+sZlxRNl8/cbzXcUIW7En0Sv9w7k8AL5rZHjyYXtZL5SU5vL2xlq5uR3KMnfASkehw2yvr2Va3jwe+enRMX1MW7IyE5/nv/tDMXgXygOfCniqKlZfm0t7Zzda6lphtdoqId9bu2svCNzbx2VkjOCaKJ4sKRMhDsTvnXg9nkFgxodQ3G+G6qkYVEBEJSntnN995aCV5man8e4yeOO8toHMgZvam/2ejme3tdWs0s72RjXhowjkaL/gOYQGsr2oMy/pEJHHc+vJ6Pty5l5+cP42C7DSv4xyygAqIc26e+eZVnOKcG9TrluucGxThjIcknNeBAGSnpzA8P1Mn0kUkKMu37uH21zZw4awRnD5liNdxwiLgXljOOQc8HsEsMaO8NIf11SogIhKYlvZOvvvwSobmZfKDBYd5HSdsgu3G+66ZHRmRJDFkQmkuG2ua6Op2XkcRkRjwX099yObdzfzywunkZqR6HSdsgi0gJ+Gb0nZjIg6m2KO8JIf2zm4qdEGhiBzE4+9X8uCSbVxz0viY73W1v7gfTDES/tETq4lxxTkepxGRaLWhuombHv+AOWMKuP7UCV7HCbtgrwOpMLPBQDmQ0eulirCminJl/p5YG6obgfg4GSYi4bWvvYtr7ltORmoyt140g5TkYA/4RD8NphgC9cQSkQNxznHTE6v5uKqRu6+Yw5C8jIP/UgzSYIohmlCawzpdCyIifbjzrS38dfl2vn1qOSdMKPY6TsTE/WCK4b6QsMeE0lw21TTT2dUd1vWKSGx7c/1u/vvpD5k/ZQjXnVzudZyICraA7D+Y4pNE+WCK4b6QsEd5aS7tXb4xsUREACpqm7nm/uWUl+Tyq89Nj7kZBoOlwRRDNKHUdyJdPbFEBKC+pZ3L71qCGfzxS7PJTg95qMGYEVQLxMyuN7MR4BtM0Tn3N+dce2SiRbfxxRoTS0R82jq7uOr/llFZt4+Fl85mVGGW15EGRLCHsAYBz5vZIjO7xsxKIxEqFmSnpzC6MIsPdoT33IqIxJbubse/PrKK9zbX8cvPTWfO2AKvIw2YoAqIc+4W59wU4BpgGPC6mb0UkWQxYPboApZs2YNvmDARSUQ/e34tf1u5g+/Pn8jZ04d5HWdAhXplSzWwC6gFSsIXJ7YcNbaAuuZ2NtboehCRRPT71zfyh9c38cWjR/H1E2J3atpQBXsO5Otm9hrwMlAEfNU5d3gkgsWCI/1N1fc27/E4iYgMtPsXb+Wnz65lwfRh3HL2VHwzXiSWYFsgo4FvO+emOOduds59GIlQsWJMYRbFuem8t7nW6ygiMoCeWrmDm55YzUkTi/n156aTHOfddfsTbDfeGyIVJBaZGXPGFrB4cx3OuYT8BiKSaJ5dvZNvP7SCI0cXcPsls0iNwzGuAhX37zxSV6L3mDOmgJ0NrVTu2ReR9YtI9Hjug11884H3OWJkPndefiSZacleR/JU3BeQSF2J3qOny96SLXURWb+IRIfn1+zi2vuXc/iIPO66/EhyEuBCwYMJ9iT6zwJ5LpFMLM1lUEYK721WARGJV0+t3ME19y1n6vA87rpiTlzNKngogm2BnNbHcwk9yVRSku88iAqISHx6eMk2rnvwfWaOHsz/XTmHQSoenwiogPi7764GJvqnsu25bQZWRzZi9DtyTAGbdjdT09jmdRQRCaO/vLWZ7z+2inllRdx9uVoe+wv0IN79wLPAT4DePbEanXMJ/9W793mQs6YN9TiNiBwq5xy/fnEdt72ygdMPK+W2i2eQnpLYJ8z7ElALxDnX4Jzb4py7CNgLlOK7JmSqmR0fyYCxYOrwPDJTk1m8SdeDiMS6zq5ubvzram57ZQNfOHIkt18yU8WjH5rSNgxSk5M4alwBr35cww91PYhIzGpp7+S6B1bw0kdVfPPkMr5z2gT9fz4ATWkbJmdMGcLWuhbW7tLw7iKxqLqxlS8sfJdX1lZxy9lT+O7pE1U8DiLup7QdKKdOLsXMd6GRiMSWj3c1ct7v3mZ9VRMLL53Nl+eO8TpSTIj7KW0HSnFuOkeOLuD5NSogIrHk1bXVfPaOt+no6uaRq4/h1MMSdpqjoAU7H8h5zrl659wPgf8E/gycG4lgseiMqUNYu6uRLbubvY4iIgfhnOMPr2/kiruXMKowiyeuOZapwyMzYkW8Cnkok0Sf0rYvp/u/uagVIhLdWju6+O7DK/nJs2s5a9pQHr16LsPyM72OFXOC7YWVDlwAjOn9u865/wpvrNg0siCLqcMH8dyaXXwtASeXEYkF2+pauPreZXy4cy/fPW0C155cppPlIQq2BfIkcA7QCTT3ukWtSI/Gu78zDhvC+1vrqdrbOiDbE5HAvb6uhs/c9ibb6lq488tH8s1TylU8DkGwBWSEc+7zzrmfO+d+1XOLSLIwifRovPubP3UIAC/oMJZI1Ojq9l1Zftlf3mNoXgZPfXMeJ01K2Nm4wybYAvK2mU2LSJI4UVaSw7jibJ5evdPrKCKC7/qOS/+8mFtfXs/5M0bw+DeOZXRhttex4kJA50D8Ayk6//KXm9kmoA0wwCXyvOj7MzPOnzGcX76wjs27mxlbpD9UEa+8uX431z+8gsbWDn7+2cP53OyRXkeKK4GeRP9MRFPEmQtnj+Q3L63nwfe2cuNZk72OI5JwOrq6+dUL6/jDGxsZX5zDvVcexcQhuV7HijsBFRDnXAWAmWUA3wDm4WuRvAncEbF0Map0UAanTi7hkWWVfOf0CRqITWQAbd7dzLcffJ+VlQ1cNGcUP/jMYQk/9WykBHsO5B5gCnAb8L/AZOD/wh0qHlx81Gjqmtt5fk2V11FEEoJzjvsWV3DWbxexpbaFOy6ZyU/On6biEUHBTuo70Tk3vdfjV81sZTgDxYvjyooYMTiTBxZv5ezpw7yOIxLXqhtbueGx1byytpp5ZUX88sLpDMnL8DpW3Au2BfK+mR3d88DMjgLeCm+k+JCUZFw0ZxTvbKplU02T13FE4pJzjidXbOf037zBWxt2c/OCw7jnijkqHgMk2AJyFL6uvFvMbAu+uUBOMLPVZrYq7Oli3IWzR5CSZDzw3lavo4jEnd1NbXzjvuV868EVjCnM5unrjuPyY8eSlKQLAwdKsIew5kckRZwqyc3g9CmlPLRkG9edUq75lEXCwDnHEyu2c8tTH9LS1sUNZ07iK/PGkpIc8tB+EqJgR+OtAPKBBf5bvnOuoucWiYCx7usnlLG3tZN73tHuETlU2+v3cfldS7j+oZWMK8rmmW/N4+oTxqt4eCSovW5m3wLuA0r8t3vN7JuRCBYvpo3I48SJxfxp0Saa2zq9jiMSkzq7uvnTok2c9uvXeW9zHT9ccBiPXD2XshJd2+GlYMv2lcBRzrkfOOd+gG9O9K+GP1Z8+ebJ5exp6eD+xToXIhKsVZX1nHv7W/z46Y84elwhL1x/PJcdO5ZknevwXLDnQAzo6vW4y/+cHMCs0YOZV1bEH97YxKXHjCYjVf3SRQ6moaWDX7ywlvsWb6UoJ53fXTyTs6YN0ei5USTYAvIXYLGZPe5/fC6+WQnlIL55chmfX/guD763lcuOHet1HJGo1d3teGx5JT99di17Wtq5bO4Yrj9tAoPUCSXqBDqYYopzrtM592szew3fUCYGXO6cez+SAePFUeMKmTO2gDte38iFs0eSnR5s7RaJfyu31XPz39awYls9M0blc8+Vc5gyTNPMRqtAP8XeA2YCOOeWA8sjliiO/dv8SVxwx9vc8dpGvnfGRK/jiESN6sZWfvHcxzyyrJKinHR+eeF0zp8xXNd0RLlAC4j+FcNg1ujBnDdjOAsXbeJzs0cyqjDL60ginmrt6OLPb27m9lc30N7VzVePG6trpmJIoAWk2My+09+LzrlfhylP2JnZAmBBWVmZ11EAuOHMSTy/Zhc/fvpDFn5pttdxRDzR3e14atUOfv7cx2yv38dph5Xy72dN1vw5MSbQbrzJQA6Q288tag30lLYHUzoog2tOKuOFD6t4c/1ur+OIDLh3N9Vy7u1v8a0HV5CXmcp9XzmKP35ptopHDAq0BbLTOfdfEU2SQK6cN5aHlmzjh0+t4enr5mm+EEkIH+3cy8+fW8urH9cwNC+DX104nfN0niOmBdoC0b9wGGWkJnPL2VPYUN3Eb19a73UckYjaWtvC9Q+t4KxbF7GsYg83nDmJV793IhfMGqHiEeMCbYGcEtEUCeikSSV8bvYIfv/6Rk6ZXMqs0YO9jiQSVrsaWrntlfU8tGQbyUnGVceN4xsnlpGXpRPk8SLQKW3rIh0kEf3nZw7jrQ21fO+RlTx93Tyy0nRtiMS+6sZW7nhtI/ct3opzjovmjOLak8soHaQ5OuKNPrE8lJuRyi8uPJyL/7iYnz27llvOmep1JJGQVTe2svD1Tdy7uIKOLsf5M4Zz3SnljCxQd/V4pQLisbnji7hs7hjuensLc8uKOGPKEK8jiQSlam8rv399I/cv3kpHVzfnzhjOdSeXM0a9quKeCkgUuOHMSby/dQ/ffXgl46/Joawkx+tIIge1ra6FP7yxkYeXVtLV7WtxfOOkMnXHTSDmnPM6w4CYPXu2W7p0qdcx+rWjfh8LbnuT/KxUnrx2HjkaK0ui1LqqRn7/+kaeXLGDZDMumDWcr59QppEV4pSZLXPO9XnVsz6losSw/Exuu3gGl/75Pb738Epuv2SmujhK1HDOsbRiD79/bSMvr60mMzWZy+aO4avHjWNInk6OJyoVkCgyd3wRN545iR8//RE/fW4t/37WZK8jSYLr6nY8v2YXC9/YxIpt9QzOSuX6UyfwpWNGMzg7zet44jEVkChz5byxbK1rYeEbmyjKSeOq48d7HUkSUFNbJ48s3cZf3trC1roWRhVkccvZU7hw9gh1N5dP6C8hypgZNy+YQm1TO//zzFqKctI5f+YIr2NJgthW18Jdb2/h4SXbaGzrZOaofP79rEmcdtgQTSEr/0QFJAolJxm//vx09rS08/1HV5GZmsyZ04Z6HUvilHOONzfs5u63t/Dy2mqSzThr2lCumDeWI0bmex1PopgKSJRKT0nmD5fO4rK/LOHaB97n113dnHPEcK9jSRxp2NfBY8squXdxBZtqminMTuOaE8u45OhRDM3L9DqexAAVkCiWm5HKPVfM4cq7l/Dth1bQ1tnN52aP9DqWxDDnHCsrG3hg8VaeXLmd1o5ujhiZz68unM6/HD6UjFSNDC2BUwGJctnpKfzlsjlc9X9L+f6jq6hrbudrx4/DTMejJXB7Wzv424odPPDeVtbs2EtmajLnTB/OpceMZurw6JgrR2KPCkgMyExL5k9fns13H17JT59dy+aaZn507lTSUgIdjV8SkXOOxZvreHjpNp5ZvZPWjm4mDx3Ej86dyrlHDNO0sXLIVEBiRHpKMrd+YQbjirK59ZUNbK1r4fZLZqovvvyT7fX7eGxZJY8uq2RrXQs56SmcP3MEXzhyJNOG56n1KmGjAhJDkpKM75w+kdGF2dz419X8y62L+N9LZjJzlOYSSXSNrR08+8Eu/rq8knc3+WZfOGZcId8+tZz5U4fo2g2JCP1VxaALZo2gvDSHa+5fzud+/w7/Nn8SXzlurL5ZJpj2zm5eX1fDEyu289KHVbR1djOmMIvvnDaB82YM1zDqEnEaTDGGNezr4PuPruT5NVUcV17ETy84nOH56n4Zzzq7unl3Ux1/X7WDZz/YRcO+Dgqy0/iXaUM5d8YwZo4arC8SElYHGkxRBSTGOee4d/FWfvLMRySZcdO/TOYLR47Uh0gc6ezqZvHmOp5ZvZPn1+xid1M72WnJnHZYKeccMZx55UWkJqtDhUSGCgjxW0B6bKtr4fuPruKdTbXMGVvALWdPYfLQQV7HkhC1dXbx9oZanl+zixc+rKKuuZ3M1GROnlzCgsOHcuLEEl2zIQNCBYT4LyAA3d2Oh5Zu4+fPraVhXweXHj2a60+bQH6WemrFgoaWDl5bV80LH1bx+sc1NLV1kpuewkmTSjhr2hBOmFBCZpqKhgwsFRASo4D0qG9p59cvruPedyvITk/ha8eP4/Jjx5KtSaqiinOOjTXNvLq2mpfXVrF0yx46ux1FOemcOrmEM6YOYe74QtJTVDTEOyogJFYB6fHxrkZ++cLHvPhhFYXZaXzthHFcNGeULiDzUFNbJ+9urOW1ddW8vq6GbXX7AJg0JJeTJpVw2mGlHDEiX5OJSdRQASExC0iP5Vv38KsXPuatDbXkZqTwxaNHc/ncMZQM0kxykdbZ1c2q7Q28tX43izbsZnmFr5WRlZbM3PFFnDCxmJMmFjNisLrcSnRSASGxC0iPldvq+cMbG3n2g10km3HGlCFccqnmv2MAAA0lSURBVPQojhlXqF5bYdLV7fhwx14Wb67lnY21LN5cR1NbJwBThw9iXlkxx5UXMXvMYB2akpgQlwXEzMYBNwF5zrnPHmx5FZB/2LK7mXvfreCRZZU07OtgbFE25x4xnHNnDGN0YbbX8WJKa0cXK7fVs7RiD0u21LFsyx4a/QVjTGEWc8uKmDu+kGPGFVKYk+5xWpHgRV0BMbM7gc8A1c65qb2enw/8FkgG/uSc+2kA63pUBSQ0rR1d/H3VTh5bVsm7m2txDo4Ymc8ZU4ZwxpRSxhXneB0xqnR3OzbXNrOqsp4VW+tZvrWej3bupbPb939oQmkOs8cUcNTYAo4aW8iQPB0ilNgXjQXkeKAJuKengJhZMrAOOA2oBJYAF+ErJj/ZbxVXOOeq/b+nAhIGO+r38eSKHTyzeiertzcAML44m+MnFHN8eTFHjStIqPGUOru62VLbzJode/lgewNrduxl9fYGGlt9rYustGSmj8hn5uh8Zo4azKzRg9VdWuJS1BUQADMbA/y9VwE5Bvihc+4M/+MbAZxz+xeP/dfTbwExs6uAqwBGjRo1q6KiImz549n2+n28sGYXr6yt5r3NdbR1dpOSZEwdnsecsQXMHj2Yw0fkUzooPebPnXR1O3bU72NDdRPrqhpZV9XEx1V7WVfVRHtnNwBpKUlMGpLLtOF5TB+Rz+Ej8ygrziFFV39LAoiVAvJZYL5z7iv+x5cCRznnru3n9wuB/8bXYvnTwQqNWiChae3oYsmWOt7eWMuSzXWsqmygvcv3wVqcm8604XlMKM1lQmkOE0pzGV2YFXXdhFs7uthRv4/t9fuoqG1ha10LFbXNbN7dzJbalk8KBfje06QhuUwaksvEIYOYMmwQZSU5GipEEtaBCkg0HZPo66tsv9XNOVcLXB25OAKQkZrMceXFHFdeDPg+jNfsaGB1ZQOrtjewZvteFq2voaPrH/9URTlpjCrIYlh+JkPzMhiSl0lRThqF2ekUZKeRl5VKbkYKOWkpIV3v0NbZRUtbF01tndS3dNCwr4M9Le3UNbdT29RGTVMbVXvbqNrbyq6GVmqb2z/1+2kpSYwYnMm4omxOnFjC2KJsyktyKCvJ0WEokSBEUwGpBHpP+D0C2OFRFulHRmoys0YXMGt0wSfPdXR1s2V3M+urm6io9X27r6htYc2OvbzoH2a8L2aQkZJMRmoSGanJpCQbyWYk+4tKt4Nu5+jscrR1dtHW0U1rZ9enitX+kgwKstMoHZRB6aAMDh+Rx/D8TIblZzI8P5NRhVmU5mboQj2RMIimArIEKDezscB24AvAxd5GkkCkJidRXppLeWnuP73mnKO+pYPa5jZqm3ythIZ9HTS2dtLY2sG+ji5aO7pp7eiis9vR5b9hkGRGskFyUhLpqUmkp/gKTU56Cllpvp95mankZaaSn5VGUU4a+VlpnxQgEYksTwqImT0AnAgUmVklcLNz7s9mdi3wPL6eV3c659Z4kU/Cx8wYnJ3G4Ow0ykq8TiMi4eRJAXHOXdTP888Az4RzW2a2AFhQVlYWztWKiCS8uO9a4px7yjl3VV5entdRRETiStwXEBERiQwVEBERCYkKiIiIhCTuC4iZLTCzhQ0NDV5HERGJK3FfQHQSXUQkMuK+gIiISGTE7IRSwTKzGiDY4XiLgN0RiHOoojUXKFuolC000ZotWnNB8NlGO+eK+3ohYQpIKMxsaX+jUHopWnOBsoVK2UITrdmiNReEN5sOYYmISEhUQEREJCQqIAe20OsA/YjWXKBsoVK20ERrtmjNBWHMpnMgIiISErVAREQkJCogIiISEhWQXszsR2a2ysxWmNkLZjasn+W+bGbr/bcvD0CuX5jZWn+2x80sv5/ltpjZan/+pZHOFWS2+Wb2sZltMLMbBijbhWa2xsy6zazfbose7bdAs3mx3wrM7EX/3/eLZja4n+W6/PtshZn9LYJ5DrgPzCzdzB7yv77YzMZEKksI2S4zs5pe++krA5TrTjOrNrMP+nndzOxWf+5VZjYzpA0553Tz34BBve5fB/y+j2UKgE3+n4P99wdHONfpQIr//s+An/Wz3BagaID32UGz4ZthciMwDkgDVgKHDUC2ycBE4DVg9gGW82K/HTSbh/vt58AN/vs3HODvrWkAshx0HwDf6Pm/im8q7IcG6N8wkGyXAf87kH9b/u0eD8wEPujn9bOAZwEDjgYWh7IdtUB6cc7t7fUwG+irh8EZwIvOuTrn3B7gRWB+hHO94Jzr9D98FxgRye0FI8Bsc4ANzrlNzrl24EHgnAHI9pFz7uNIbycUAWbzZL/5t3G3//7dwLkDsM3+BLIPeud9FDjFzCxKsnnCOfcGUHeARc4B7nE+7wL5ZjY02O2ogOzHzP7bzLYBlwA/6GOR4cC2Xo8r/c8NlCvwfXPoiwNeMLNlZnbVAGbq0V82r/fZwXi93/rj1X4rdc7tBPD/7G82+wwzW2pm75pZpIpMIPvgk2X8X2YagMII5Qk2G8AF/sNEj5rZyAHIFYiw/G15Mie6l8zsJWBIHy/d5Jx70jl3E3CTmd0IXAvcvP8q+vjdQ+4LfbBc/mVuAjqB+/pZzbHOuR1mVgK8aGZr/d9EvM4WkX0WaLYAeLbfDraKPp6L+H4LYjWj/PttHPCKma12zm0MR75eAtkHEdtPBxHIdp8CHnDOtZnZ1fhaSidHPNnBhWWfJVwBcc6dGuCi9wNP888FpBI4sdfjEfiOY0c0l/9k/WeAU5z/IGYf69jh/1ltZo/ja2If8gdhGLJVAr2/eY0AdhxqrkCyBbgOT/ZbADzZb2ZWZWZDnXM7/Yc1qvtZR89+22RmrwEz8J0TCKdA9kHPMpVmlgLkceDDNwOWzTlX2+vhH/GdJ4wGYfnb0iGsXsysvNfDs4G1fSz2PHC6mQ3290453f9cJHPNB/4NONs519LPMtlmlttz35+rzx4YA50NWAKUm9lYM0vDd6IzYr12guHVfguQV/vtb0BP78IvA//UWvL//af77xcBxwIfRiBLIPugd97PAq/09yVroLPtd17hbOCjAcgViL8BX/L3xjoaaOg5bBmUge4dEM034DF8Hx6r8DU9h/ufnw38qddyVwAb/LfLByDXBnzHK1f4bz09ToYBz/jvj8PXC2QlsAbfYZKB2GcHzeZ/fBawDt831IHKdh6+b1ptQBXwfBTtt4Nm83C/FQIvA+v9Pwv8z3/y/wCYC6z277fVwJURzPNP+wD4L3xfWgAygEf8f4vvAeMGYj8FmO0n/r+rlcCrwKQByvUAsBPo8P+dXQlcDVztf92A3/lzr+YAvRQPdNNQJiIiEhIdwhIRkZCogIiISEhUQEREJCQqICIiEhIVEBERCYkKiIiIhEQFREREQqICInGp11wVH5jZI2aW5UGGpgis8+0gl/+hmX0v3DlEQAVE4tc+59wRzrmpQDu+q3Ajxj8kRMT/Pznn5kZ6GyKBUgGRRLAIKAMws+/4WyUfmNm3/c9938yu89//jZm94r9/ipnd67//RTN7z9+q+YOZJZvZGDP7yMxuB5bz6cHpPsXMnvAPF7+mZ8j4QLbbx3qaem33j/71vWBmmb2Wucl8s+S9hG/Sqp7n+3oPR/qHGs/wjwu2xsym9rHdx83sx2a2yMx2mdkhD2IpsU8FROKaf3TWM4HVZjYLuBw4Ct8sbF81sxn4Rt49zv8rs4EcM0sF5gGLzGwy8Hl8w74fAXThmy8GfB/Q9zjnZjjnKg4Q5Qrn3Cz/+q8zs8KDbfcgb60c+J1zbgpQD1zgf7+z8A3qNwM4HzjS/3yf78E5twTfwHo/xjcT4b3Oub4Gk5wK1DvnjsM3A+AlfSwjCSbhhnOXhJFpZiv89xcBfwa+DjzunGsGMLO/4vsAvwOY5R+Vtw1fa2K2/7XrgFOAWcAS8010l4lviPM3gArnm9HtYK4zs/P890fiKwDLDrLdA9nsnOt5f8uAMf77x/nfY4v/PfaMDtvfewDf4H9LgNa+tus/f5QH/Mb/VAq+oiUJTgVE4tU+/zftT5j1Pc2pc67DzLbga528jW805pOA8fiG3z4FuNs5d+N+6xsDNB8siJmdCJwKHOOca/HPnZERwHYPpK3X/S58BeGTt9RXjL7eg18BkAOk4hvZdv/3NAVY5pzr8j8+nOgZ8l48pENYkkjeAM41syz/3B/n8Y9DRW8A3/P/XITvpPsK5xuu+mXgs+absRAzKzCz0UFsNw/Y4y8ek/AdPuudqb/thvoezzOzTH/LZoH/+QO9h4XAf+KbTbKvCY+m4huqv8fh+IqdJDgVEEkYzrnlwF345oxYjG9ui/f9Ly8ChgLvOOeq8B3OWeT/vQ+B/8A3b/oq4EX/soF6Dkjx/+6PgN6HvPrdbij87/EhfB/4jx3sPZjZl4BO59z9wE+BI81s/ylXp/HpAjIVtUAENB+IiIiERi0QEREJiQqIiIiERAVERERCogIiIiIhUQEREZGQqICIiEhIVEBERCQk/x/IaaPWT9oyewAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.special import gamma\n",
    "def sigma2_mathematica(n,R):\n",
    "    return 9*(n + 1)*2**(-n - 1)*R**(-n - 3)*np.sin(np.pi*n/2)*gamma(n - 1)/(np.pi**2*(n-3))\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "nrange = np.linspace(-3,1,101, endpoint=False)\n",
    "plt.plot(nrange, sigma2_mathematica(nrange,1))\n",
    "plt.xlabel('Power law index $n$')\n",
    "plt.ylabel('Tophat variance $\\sigma^2(R=1)$');\n",
    "plt.yscale('log')\n",
    "\n",
    "plt.savefig('scale_free_sigma8.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For plotting purposes, we avoided exact integers in the $n$ range and used $A=1$.\n",
    "\n",
    "So, to find a numerically stable variant of the above expression, we need to re-write it to avoid evaluating $\\Gamma(n-1)$ at negative integers.\n",
    "\n",
    "Using some properties of the $\\Gamma$ function (<https://dlmf.nist.gov/5.5>), we can arrive at the following expression:\n",
    "\n",
    "$$\n",
    "\\sigma^2(R) = \\frac{9 A R^{-(3+n)}}{2 \\pi^{3/2}} \\frac{\\Gamma[(3+n)/2]}{\\Gamma[(2-n)/2] (n-3)(n-1)}\n",
    "$$\n",
    "\n",
    "It's worth noting that this expression is exact; we have not used any approximations.\n",
    "\n",
    "Note: a previous version of this notebook had two expressions, one for even poles and one for odd, which are now combined into one expression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigma2_robust(n,R):\n",
    "    return 9*R**-(3 + n)/(2*np.pi**(3/2)) * gamma((3 + n)/2)/(gamma((2 - n)/2)*(n - 3)*(n - 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "def test_all():\n",
    "    n = np.array([-2.9, -2.5, -2, -1.3, -1, -0.01, 0, 0.001, 0.75])\n",
    "    R = 8.\n",
    "    \n",
    "    res = pd.DataFrame(index=n)\n",
    "    res.index.name = 'n'\n",
    "    res['Naive'] = sigma2_mathematica(n, R)\n",
    "    res['Stable'] = sigma2_robust(n, R)\n",
    "    \n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/lgarrison/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:4: RuntimeWarning: invalid value encountered in multiply\n",
      "  after removing the cwd from sys.path.\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Naive</th>\n",
       "      <th>Stable</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>n</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>-2.900</th>\n",
       "      <td>0.432508</td>\n",
       "      <td>0.432508</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>-2.500</th>\n",
       "      <td>0.047497</td>\n",
       "      <td>0.047497</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>-2.000</th>\n",
       "      <td>-inf</td>\n",
       "      <td>0.011937</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>-1.300</th>\n",
       "      <td>0.002945</td>\n",
       "      <td>0.002945</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>-1.000</th>\n",
       "      <td>NaN</td>\n",
       "      <td>0.001781</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>-0.010</th>\n",
       "      <td>0.000471</td>\n",
       "      <td>0.000471</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.000</th>\n",
       "      <td>NaN</td>\n",
       "      <td>0.000466</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.001</th>\n",
       "      <td>0.000466</td>\n",
       "      <td>0.000466</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.750</th>\n",
       "      <td>0.000392</td>\n",
       "      <td>0.000392</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           Naive    Stable\n",
       "n                         \n",
       "-2.900  0.432508  0.432508\n",
       "-2.500  0.047497  0.047497\n",
       "-2.000      -inf  0.011937\n",
       "-1.300  0.002945  0.002945\n",
       "-1.000       NaN  0.001781\n",
       "-0.010  0.000471  0.000471\n",
       " 0.000       NaN  0.000466\n",
       " 0.001  0.000466  0.000466\n",
       " 0.750  0.000392  0.000392"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res = test_all()\n",
    "res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\\begin{tabular}{lrr}\n",
      "\\toprule\n",
      "{} &       Naive &      Stable \\\\\n",
      "n      &             &             \\\\\n",
      "\\midrule\n",
      "-2.900 &    0.432508 &    0.432508 \\\\\n",
      "-2.500 &   0.0474965 &   0.0474965 \\\\\n",
      "-2.000 &        -inf &   0.0119366 \\\\\n",
      "-1.300 &  0.00294465 &  0.00294465 \\\\\n",
      "-1.000 &         nan &  0.00178104 \\\\\n",
      "-0.010 &  0.00047106 &  0.00047106 \\\\\n",
      " 0.000 &         nan & 0.000466274 \\\\\n",
      " 0.001 & 0.000465801 & 0.000465801 \\\\\n",
      " 0.750 & 0.000392073 & 0.000392073 \\\\\n",
      "\\bottomrule\n",
      "\\end{tabular}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(res.to_latex(float_format=lambda s:'{:.6g}'.format(s)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}