{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "source": [ "import numpy\n", "import tqdm\n", "import wendy\n", "%pylab inline\n", "import matplotlib.animation as animation\n", "from matplotlib import cm\n", "from IPython.display import HTML\n", "_SAVE_GIFS= False\n", "rcParams.update({'axes.labelsize': 17.,\n", " 'font.size': 12.,\n", " 'legend.fontsize': 17.,\n", " 'xtick.labelsize':15.,\n", " 'ytick.labelsize':15.,\n", " 'text.usetex': _SAVE_GIFS,\n", " 'figure.figsize': [5,5],\n", " 'xtick.major.size' : 4,\n", " 'ytick.major.size' : 4,\n", " 'xtick.minor.size' : 2,\n", " 'ytick.minor.size' : 2,\n", " 'legend.numpoints':1})\n", "numpy.random.seed(2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Adiabatic contraction\n", "\n", "Adiabatic contraction happens during galaxy formation when a cold gaseous and stellar disk grows slowly inside a dark-matter halo. This slow growth causes the dark-matter halo to adiabatically contract and thus become more dense near the center. We can illustrate this process in one dimension by starting from a self-gravitating disk and slowly growing an equal-mass disk inside this disk with a velocity dispersion that is 10x smaller than the original disk. First we initialize both of these disk:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Initial disk\n", "N= 10000\n", "# compute zh based on sigma and totmass\n", "totmass= 1. # Sigma above\n", "sigma= 1.\n", "zh= sigma**2./totmass # twopiG = 1. in our units\n", "tdyn= zh/sigma\n", "x= numpy.arctanh(2.*numpy.random.uniform(size=N)-1)*zh*2.\n", "v= numpy.random.normal(size=N)*sigma\n", "v-= numpy.mean(v) # stabilize\n", "m= totmass*numpy.ones_like(x)/N" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Properties of colder disk grown within\n", "N_cold= 10000\n", "totmass_cold= 1. # Sigma above\n", "sigma_cold= 0.1\n", "zh_cold= sigma_cold**2./totmass_cold # twopiG = 1. in our units\n", "x_cold= numpy.arctanh(2.*numpy.random.uniform(size=N_cold)-1)*zh_cold*2.\n", "v_cold= numpy.random.normal(size=N_cold)*sigma_cold\n", "v_cold-= numpy.mean(v_cold) # stabilize\n", "m_cold= totmass_cold*numpy.ones_like(x_cold)/N_cold" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we run the $N$-body integration, adding a small fraction of the cold disk in each time step:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████| 2500/2500 [01:28<00:00, 23.43it/s]\n" ] } ], "source": [ "nt= 2500\n", "xt= numpy.empty((N,nt+1))\n", "vt= numpy.empty((N,nt+1))\n", "xt[:,0]= x\n", "vt[:,0]= v\n", "xt_cold= numpy.zeros((N_cold,nt+1))\n", "vt_cold= numpy.zeros((N_cold,nt+1))\n", "xt_cold[:,0]= x_cold\n", "vt_cold[:,0]= v_cold\n", "for ii in tqdm.trange(nt):\n", " # Add mass\n", " skip= N_cold//nt\n", " tx= numpy.concatenate((x,x_cold[:(ii+1)*skip]))\n", " tv= numpy.concatenate((v,v_cold[:(ii+1)*skip]))\n", " tm= numpy.concatenate((m,m_cold[:(ii+1)*skip]))\n", " g= wendy.nbody(tx,tv,tm,0.05*tdyn,approx=True,nleap=10)\n", " tx,tv= next(g)\n", " xt[:,ii+1]= tx[:N]\n", " vt[:,ii+1]= tv[:N]\n", " xt_cold[:(ii+1)*skip,ii+1]= tx[N:N+(ii+1)*skip]\n", " vt_cold[:(ii+1)*skip,ii+1]= tv[N:N+(ii+1)*skip]\n", " # Re-set x,v\n", " x= tx[:N]\n", " v= tv[:N]\n", " x_cold[:(ii+1)*skip]= tx[N:N+(ii+1)*skip]\n", " v_cold[:(ii+1)*skip]= tv[N:N+(ii+1)*skip]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The evolution of the system in phase space clearly shows the contraction. Note how the total area enclosed by each orbit remains approximately the same, because the contraction is adiabatic:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def init_anim_frame():\n", " line1= plot([],[])\n", " xlabel(r'$x$')\n", " ylabel(r'$v$')\n", " xlim(-7.99,7.99)\n", " ylim(-5.99,5.99)\n", " return (line1[0],)\n", "figsize(6,4)\n", "fig, ax= subplots()\n", "c= wendy.energy(xt[:,0],vt[:,0],m,individual=True)\n", "s= 5.*((c-numpy.amin(c))/(numpy.amax(c)-numpy.amin(c))*2.+1.)\n", "line= ax.scatter(x[::N//5000],v[::N//5000],c=c[::N//5000],s=s[::N//1000],edgecolors='None',cmap=cm.jet)\n", "txt= ax.annotate(r'$t=%.0f$' % (0.),\n", " (0.95,0.95),xycoords='axes fraction',\n", " horizontalalignment='right',verticalalignment='top',size=18.)\n", "subsamp= 4\n", "def animate(ii):\n", " line.set_offsets(numpy.array([xt[::N//5000,ii*subsamp],vt[::N//5000,ii*subsamp]]).T)\n", " txt.set_text(r'$t=%.0f$' % (ii*subsamp/20.))\n", " return (line,)\n", "anim = animation.FuncAnimation(fig,animate,init_func=init_anim_frame,\n", " frames=nt//subsamp,interval=40,blit=True,repeat=True)\n", "if _SAVE_GIFS:\n", " anim.save('adiabcontract_phasespace.gif',writer='imagemagick',dpi=80)\n", "# The following is necessary to just get the movie, and not an additional initial frame\n", "plt.close()\n", "out= HTML(anim.to_html5_video())\n", "plt.close()\n", "out" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "The density itself clearly becomes more peaked near zero over time. The initial, self-gravitating $\\mathrm{sech}^2$ disk is shown as the dashed curve." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "figsize(6,4)\n", "fig, ax= subplots()\n", "ii= 0\n", "a= ax.hist(xt[:,ii],bins=31,histtype='step',lw=1.,color='k',range=[-8.,8.],weights=31./16./N*numpy.ones(N))\n", "xs= numpy.linspace(-8.,8.,101)\n", "ax.plot(xs,totmass/4./zh/numpy.cosh(xs/2./zh)**2.,'b--',lw=2.,zorder=0)\n", "ax.set_xlim(-8.,8.)\n", "ax.set_ylim(10.**-3.,1.)\n", "ax.set_xlabel(r'$x$')\n", "ax.set_ylabel(r'$\\rho(x)$')\n", "ax.set_yscale('log')\n", "ax.annotate(r'$t=0$',(0.95,0.95),xycoords='axes fraction',\n", " horizontalalignment='right',verticalalignment='top',size=18.)\n", "subsamp= 4\n", "def animate(ii):\n", " ax.clear()\n", " a= ax.hist(xt[:,ii*subsamp],bins=31,histtype='step',lw=1.,color='k',range=[-8.,8.],weights=31./16./N*numpy.ones(N))\n", " xs= numpy.linspace(-8.,8.,101)\n", " ax.plot(xs,totmass/4./zh/numpy.cosh(xs/2./zh)**2.,'b--',lw=2.,zorder=0)\n", " ax.set_xlim(-8.,8.)\n", " ax.set_ylim(10.**-3.,1.)\n", " ax.set_xlabel(r'$x$')\n", " ax.set_ylabel(r'$\\rho(x)$')\n", " ax.set_yscale('log')\n", " ax.annotate(r'$t=%.0f$' % (ii*subsamp/20.),\n", " (0.95,0.95),xycoords='axes fraction',\n", " horizontalalignment='right',verticalalignment='top',size=18.)\n", " return a[2]\n", "anim = animation.FuncAnimation(fig,animate,#init_func=init_anim_frame,\n", " frames=nt//subsamp,interval=40,blit=True,repeat=True)\n", "if _SAVE_GIFS:\n", " anim.save('adiabcontract_density.gif',writer='imagemagick',dpi=80)\n", "# The following is necessary to just get the movie, and not an additional initial frame\n", "plt.close()\n", "out= HTML(anim.to_html5_video())\n", "plt.close()\n", "out" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "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.5.3" } }, "nbformat": 4, "nbformat_minor": 2 }