{
"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
}