ex07_density_perturbations_evolution_tophat.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Astro 304 \"Galaxies\"\n",
    "\n",
    "### Density perturbations and their evolution in expanding universe: non-linear evolution in the top-hat collapse model\n",
    "\n",
    "<img size=700 src=\"http://astro.uchicago.edu/~andrey/classes/a304/fig/spherical_cow.png\"></img>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "You can turn this notebook into a slide show using nbconvert package:\n",
    "\n",
    "jupyter-nbconvert --to slides ex07_density_perturbations_evolution_tophat.ipynb --post serve\n",
    "\n",
    "For many of the image links to work, you need to be connected to internet. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# setup notebook for inline figures\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "\n",
    "# import pyplot and set some parameters to make plots prettier\n",
    "import matplotlib.pyplot as plt\n",
    "from code.plot_utils import plot_pretty\n",
    "plot_pretty()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### The top-hat spherical collapse model\n",
    "\n",
    "Consider a spherically-symmetric density fluctuation of initial physical radius $R_{\\rm\n",
    "  i}$, with a spherically-symmetric distribution of matter around it. Although this is an example of a \"spherical cow\" model, this particular model turns out to be remarkably useful. It does have limitations that we will consider in due time. \n",
    "\n",
    "If the density within $R=R_i$ is constant or increasing with decreasing $r\\leq R$, the mass within $r$ will be conserved during collapse because inner shells will collapse on the same or faster time scale. Recall that *free-fall* time scales as $t_{\\rm ff}\\propto \\rho^{-1/2}$. The fluctuation in which $\\rho(<r)=\\rm const$ for all $r\\leq R_i$ is called **top-hat fluctuation**.\n",
    "\n",
    "For any given $r\\leq R$, there should be no force from the matter at radii larger than a given radius $r$ because spherical gravitating shells do not exert a gravitational force inside of them according to the <a href=\"https://goo.gl/DTmepq\">the first Newton's theorem.</a>\n",
    "\n",
    "Thus, we can consider evolution of any $r\\leq R$, as a separate similar problem. \n",
    "Further, if we can neglect pressure of matter and associated force as well as any force associated with vacuum energy, a particle at radius $r$ is accelerated only by the mass inside of $r$, i.e.:\n",
    "\n",
    "$$\\ddot{r}=-\\frac{GM(<r)}{r^2}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Spherical collapse\n",
    "\n",
    "Let's consider evolution of gravitationally bound, $E={\\rm P.E.}+{\\rm K.E.}=W+K<0$, fluctuation with $\\delta_i>0.$ Total energy is actually the first integral of our equation of motion: \n",
    "\n",
    "$$E=\\frac{\\dot{r}^2}{2}-\\frac{GM(<r)}{r}.$$\n",
    "\n",
    "Note that  $M(<r)=\\rm const$ just like gravitating mass is constant in this problem due to our assumption about $\\delta(r)$ profile. Thus, equation $\\ddot{r}=-GM(<r)/r^2$ in this problem is the same as the equation of motion of a body moving in a central gravitational field (planet around a star, satellite around the Earth, etc.).\n",
    "\n",
    "Solution of this equation is well known, as you can imagine. I include the complete derivation of the solution for the general case of a body with angular momentum moving on elliptical orbit around a central gravitating mass in the <a href=\"\">ex07b_eqs_motion_central_field</a> notebook. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Assuming that matter has no angular momentum ($L=0$) and moves only radially (eccentricity $e=1$), the solution of $\\ddot{r}=-GM(<r)/r^2$ is \n",
    "\n",
    "$$r=\\frac{r_{\\rm ta}}{2}\\,(1-\\cos\\eta); \\ \\ \\ t=\\frac{t_{\\rm ta}}{\\pi}\\,(\\eta-\\sin \\eta),$$ \n",
    "\n",
    "where $0\\leq\\eta\\leq 2\\pi$ is parametric variable (called *eccentric anomaly* in orbital dynamics), and $r_{\\rm ta}$ is the *turnaround radius* -- the maximum radius during the evolution that is reached at the turnaround time, $t_{\\rm ta}$. \n",
    "\n",
    "The evolution $r(t)$ thus consists of an initially decelerating increase of the fluctuation radius, until it reaches the maximum value, $r_{\\rm ta}$, at $t_{\\rm ta}$, and subsequent decrease of $r(t)$ at $t>t_{\\rm ta}$ until the fluctuation collapses, virializes, and settles at the final\n",
    "radius $R_{\\rm f}$ at $t=t_{\\rm coll}$. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Turnaround radius and time\n",
    "\n",
    "At $r_{\\rm ta}$ velocity of the shell is zero, so its total energy is all in the form of potential energy: \n",
    "\n",
    "$$E=-\\frac{GM}{r_{\\rm ta}},\\ \\ \\ {\\rm or}\\ \\ \\  r_{\\rm ta}=\\frac{GM}{\\vert E\\vert}.$$\n",
    "\n",
    "At the same time, the full evolution until collapse corresponds to $t=2t_{\\rm ta}$ and so perturbation will collapse ($r=0$) from $t_{\\rm ta}$ in the time equal to $t_{\\rm ta}$. We know that homegeneous sphere collapses in *free-fall time*, $t_{\\rm ff}=\\sqrt{3\\pi/[32G\\rho_{\\rm ta}]}$ where $\\rho_{\\rm ta}=3M/(4\\pi r_{\\rm ta}^3)$, and so: \n",
    "\n",
    "$$t_{\\rm ta}=\\left(\\frac{3\\pi}{32G\\rho_{\\rm ta}}\\right)^{1/2}=\\pi\\,\\left(\\frac{r_{\\rm ta}^3}{32GM}\\right)^{1/2}=\\pi\\,\\frac{GM}{(2\\vert E\\vert)^{3/2}}.$$\n",
    "\n",
    "Thus, $r_{\\rm ta}$ and $t_{\\rm ta}$ are determined by the mass $M$ and total binding energy of the perturbation $\\vert E\\vert$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The mass, in turn, can be written as a function of dimensionless *overdensity*\n",
    "\n",
    "$$\\delta=\\frac{\\rho(<r)-\\bar{\\rho}}{\\bar{\\rho}}.$$\n",
    "\n",
    "$$M=\\frac{4\\pi}{3}\\,(1+\\delta)\\,\\bar{\\rho}r^3=\\frac{4\\pi}{3}\\,(1+\\delta_i)\\,\\bar{\\rho}r_i^3={\\rm const}.$$\n",
    "\n",
    "Thus, turnaround radius and time will depend on the initial overdensity and radius, $\\delta_i$ and $r_i$ and total binding energy $\\vert E\\vert$ at $t=t_i$. These quantities constitute *initial conditions* for a given density fluctuation. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "For collapse the initial mass and radius of the perturbation should be such that it  is bound: \n",
    "\n",
    "$$\\frac{M\\dot{R}^2_{i}}{2}-\\frac{GM^2}{R_i}<0\\ \\ \\mathrm{or}\\ \\ \\frac{2GM}{R^3_i}>H^2(t_i)$$ \n",
    "\n",
    "where $H^2(t_i)=(\\dot{a}/a)^2=8\\pi G\\rho_{\\rm crit}(t)/3\\approx (\\dot{R}/R)^2_i$. If we insert expression for the mass for a given overdensity $\\delta_i$ into the boundedness condition, we see that it becomes $(1+\\delta_i)>\\rho_{\\rm crit}/\\bar{\\rho}=1/\\Omega_{\\rm m}$.  \n",
    "\n",
    "Thus, for $\\Omega_{\\rm m}=\\bar{\\rho}/\\rho_{\\rm crit}=1$, any overdensity $\\delta_i>0$ is bound and will eventually collapse. For $\\Omega_{\\rm m}<1$, however, only perturbations with $\\delta_{i}>1/\\Omega_{\\rm m}-1$ will collapse. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### An example:  \n",
    "\n",
    "Evolution of the physical radius (left) and density contrast $\\delta(<r)$ (right) as a function of time for a perturbation with comoving initial radius $R_i=1.5h^{-1}$ Mpc, mass $M\\approx 1.75\\times 10^{12}\\ M_\\odot$ in the $\\Omega_{\\rm m0}=0.3$, $\\Omega_\\Lambda=0$ cosmology ($H_0=70\\ \\rm km\\,s^{-1}\\,Mpc$).  \n",
    "\n",
    "<img width=800 src=\"http://astro.uchicago.edu/~andrey/classes/a304/fig/th_collapse_r_t.PNG\"></img>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Returning to the spherical top-hat model: collapse and \"virialization\"\n",
    "\n",
    "Physically, the final radius of the collapsed fluctuation, $R_{\\rm f}$, is assumed to be \n",
    "set by the virial relation: \n",
    "$2K=-W$ and the total energy $E=W/2$. At the turnaround $K=0$ as velocity is zero, but total energy is the same $E=W_{\\rm ta}$. \n",
    "Thus, \n",
    "\n",
    "$$\\frac{GM^2}{R_{\\rm ta}}=\\frac{GM^2}{2R_{\\rm f}}\\ \\ \\mathrm{or}\\ \\ R_{\\rm f}=\\frac{R_{\\rm ta}}{2}.$$ \n",
    "\n",
    "This is not true in cosmologies with $\\Omega_{\\Lambda}\\ne 0$ because in this case vacuum energy is acting as an anti-gravity force and initial energy is not conserved.  \n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Linear overdensity at collapse\n",
    "\n",
    "If the initial density contrast $\\delta_{\\rm i}$ would grow\n",
    "only at the linear rate, $D_+(z)$, then the density contrast at the time of\n",
    "collapse would be more than a hundred times smaller.\n",
    "\n",
    "Its value can be derived starting from the density contrast linearly extrapolated to\n",
    "the turn around epoch, $\\delta_{\\rm ta}$. One could expect that at this point the perturbation is already in the *non-linear* regime and so linearly extrapolated overdensity should be of order unity $\\delta_{\\rm ta}\\sim 1$. Because $t_{\\rm coll}=2t_{\\rm ta}$, further linear\n",
    "evolution for $\\Omega_{\\rm m}(z)=1$ until the collapse time gives\n",
    "\n",
    "$$\\delta_c=\\delta_{\\rm ta} D_+(t_c)/D_+(t_{\\rm ta})=a(t_c)/a(t_{\\rm ta})=(t_c/t_{\\rm ta})^{2/3}=2^{2/3}\\approx 1.6$$\n",
    "\n",
    "In fact, the exact calculation of the evolution equations in the case of \n",
    "$\\Omega_{\\rm m}(z)=1$ at the redshift of turnaround gives $\\delta_{ta}=1.062$\n",
    "(<a href=\"http://adsabs.harvard.edu/abs/1972ApJ...176....1G\">Gunn & Gott 1972</a>), which gives\n",
    "$\\delta_c=\\delta_{\\rm ta} 2^{2/3}\\approx 1.686$.\n",
    "\n",
    "In the right panel of the figure showing evolution above, density contrast $\\delta(<r)$ is compared to the linear evolution $\\delta_i\\,D_+(t)/D_+(t_i)$ (dashed line) correct for this cosmology and to $\\delta_i(t/t_i)^{2/3}$ expected linear growth of the $\\Omega_{\\rm m}=1$ cosmology (dark dotted line). The horizontal gray dotted line shows $\\delta_c=1.69$ -- one can see that  $\\delta_i\\,D_+(t)/D_+(t_i)=1.69$ at the time of collapse ($r_{\\rm TH}=0$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### \"Virial\" overdensity of collapsed peak\n",
    "\n",
    "The final mean internal density of a collapsed object in a $\\Omega_{\\rm m}=1$ Universe (with and $\\Omega_{\\Lambda}=0$) can be estimated\n",
    "by noting that in a $\\Omega_{\\Lambda}=0$ Universe the time interval $t_{\\rm\n",
    "  coll}-t_{\\rm ta}=t_{\\rm ta}$ should be equal to the free-fall time\n",
    "of a uniform sphere $t_{\\rm ff}=\\sqrt{3\\pi/(32G\\rho_{\\rm ta})}$, which\n",
    "means that the mean density of perturbation at turnaround is\n",
    "$\\rho_{\\rm ta}=3\\pi/(32Gt^2_{\\rm ta})$ and $\\rho_{\\rm coll}=8\\rho_{\\rm\n",
    "  ta}=3\\pi/(Gt_{\\rm coll}^2)$. \n",
    "  \n",
    " These densities can be contrasted with the mean matter densities at the corresponding time to get\n",
    "mean internal density contrasts: $\\Delta=\\rho/\\bar{\\rho}$. Mean density\n",
    "evolves as $\\bar{\\rho}=\\bar{\\rho}_{\\rm 0}a^{-3}$ where $\\bar{\\rho}_{\\rm m0}=3H_0^2/(8\\pi G)$. \n",
    "\n",
    "For $\\Omega_{\\rm m}=1$ cosmology  $a=(t/t_{\\rm U})^{2/3}$, where $t_{\\rm U}=2/(3H_0)$ is the current ($a=1$) age of universe. These equations combine to give $\\bar{\\rho}=1/(6\\pi Gt^2)$, which means that density contrast at $t=t_{\\rm coll}$, if we *assume instant virialization*, is \n",
    "\n",
    "$$\n",
    "\\Delta_{\\rm vir}\\equiv \\frac{\\rho_{\\rm coll}}{\\bar{\\rho}_{\\rm m}}=18\\pi^2=177.653.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Virial density contrast  in $\\Omega_{\\rm m}\\ne 1$ cosmologies\n",
    "\n",
    "For general cosmologies, density contrast can be computed by\n",
    "estimating $\\rho_{\\rm coll}$ and $\\bar{\\rho}_m(t_{\\rm coll})$ in a\n",
    "similar fashion. For lower $\\Omega_{\\rm m}$ models, fluctuation of the same mass\n",
    "$M$ and $\\delta$ has a larger initial radius and smaller physical\n",
    "density and, thus, takes longer to collapse. The density contrasts of\n",
    "collapsed objects therefore are larger in lower density models because\n",
    "the mean density of matter at the time of collapse is\n",
    "smaller. Accurate (to $<1\\%$ for $\\Omega_{\\rm m}=0.1-1$) approximations\n",
    "for $\\Delta_{\\rm vir}$ in open ($\\Omega_\\Lambda=0$) and flat $\\Lambda$CDM\n",
    "($1-\\Omega_\\Lambda-\\Omega_{\\rm m}=0$) cosmologies are given by eq. 6 in Bryan & Norman 1998. For example, for the concordance $\\Lambda$CDM\n",
    "cosmology with $\\Omega_{\\rm m}=0.3$ and $\\Omega_\\Lambda=0.7$ (Komatsu et al 2011, Planck 2015},\n",
    "density contrast at $z=0$ is $\\Delta_{\\rm vir}\\approx 358$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Linearly extrapolated overdensity in $\\Omega_{\\rm m}\\ne 1$ cosmologies\n",
    "\n",
    "In the\n",
    "case of $\\Omega_{\\rm m}\\ne 1$ we expect that $\\delta_{\\rm ta}$ should have\n",
    "different values. For instance, for $\\Omega_{\\rm m}< 1$ density contrast at\n",
    "the turn-around should be higher to account for the higher rate of the\n",
    "Hubble expansion. However, linear growth from $t_{\\rm ta}$ to $t_{\\rm\n",
    "  coll}$ is smaller due to the slower redshift dependence of\n",
    "$D_+(z)$. As a matter of fact, these two factors nearly cancel, so that \n",
    "$\\delta_c$ has a weak dependence on\n",
    "$\\Omega_{\\rm m}$ and $\\Omega_\\Lambda$ (e.g., <a href=\"http://adsabs.harvard.edu/abs/2005A%26A...443..819P\">Percival 2005</a>). For the concordance\n",
    "$\\Lambda$CDM cosmology at $z=0$, for example, $\\delta_c\\approx 1.675.$ The exact dependence of $\\delta_c$ on cosmological parameters is thus often neglected, but one should be aware that in some cases this may matter (see, e.g., <a href=\"http://adsabs.harvard.edu/abs/2011MNRAS.410.1911C\">Alimi et al. 2011</a>). "
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}