{
"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",
"import copy\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": [
"# The *Gaia* phase-space spiral\n",
"\n",
"One of the amazing early discoveries in the *Gaia* DR2 data set is the *Gaia* phase-space spiral. This is a spiral feature in the vertical phase-space distribution function $f(z,v_z)$ first found by [Antoja et al. (2018)](https://ui.adsabs.harvard.edu/abs/2018Natur.561..360A/abstract). In this example, we investigate how such a phase-space spiral can form from a simple perturbation to the Milky Way's disk. We also use it as an opportunity to showcase the support for arbitrary external forces and for different sorting algorithms in ``wendy``'s approximate *N*-body solution.\n",
"\n",
"A simple model for the phase-space spiral is that it results from the disk's equilibrium $f(z,v_z)$ being offset from $\\langle v_z \\rangle =0$ ([Darling & Widrow 2019](https://ui.adsabs.harvard.edu/abs/2019MNRAS.484.1050D/abstract)). As a highly simplified model for this, we initialize a self-gravitating $\\mathrm{sech}^2$ disk and only treat a fraction $\\alpha$ of that as self-gravitating masses, filling in the rest as a static, external potential. This simplification makes it straightforward to setup the equilibrium solution, because that is just that of a fully self-gravitating, $\\mathrm{sech}^2$ disk.\n",
"\n",
"First we sample $N$ points from the disk:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Initial disk\n",
"N= 100000\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": "markdown",
"metadata": {},
"source": [
"We then assign only $\\alpha$ of that disk to be self-gravitating and define the external force:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"alpha= 0.3 # \"live\" fraction\n",
"# Adjust masses to only represent alpha of the mass\n",
"m*= alpha\n",
"# 1-alpha in the mass is then given by the external force\n",
"sigma2= sigma**2.\n",
"def iso_force(x,t):\n",
" return -(1.-alpha)*sigma2*numpy.tanh(0.5*x/zh)/zh"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We then run the $N$-body simulation, using the approximate algorithm with an external force and using a fast ``quicksort`` implementation to calculate the $N$-body forces. For a simple external force implemented using ``numpy`` functions as we have done here, ``numba`` is used to compile C byte-code which is directly called in the underlying C code, allowing this external force to be very efficiently added to the system (don't worry, external forces that cannot be automatically compiled to C code are also supported, but they are slightly slower)."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"g= wendy.nbody(x,v,m,0.05*tdyn,nleap=10,approx=True,sort='quick',ext_force=iso_force)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|██████████| 1000/1000 [01:44<00:00, 9.54it/s]\n"
]
}
],
"source": [
"nt= 1000\n",
"xt= numpy.empty((N,nt+1))\n",
"vt= numpy.empty((N,nt+1))\n",
"xt[:,0]= x\n",
"vt[:,0]= v\n",
"x_init= copy.copy(x)\n",
"v_init= copy.copy(v)\n",
"for ii in tqdm.trange(nt):\n",
" tx,tv= next(g)\n",
" xt[:,ii+1]= tx\n",
" vt[:,ii+1]= tv"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We check that the original disk is indeed in equilibrium:"
]
},
{
"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",
"# 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": {},
"source": [
"Indeed, the disk is in equilibrium!\n",
"\n",
"Next, we offset all of the initial velocities by $1\\sigma$ and run the simulation again to study the effect of this perturbation. This time we use ``timsort``, a version of Python's own sorting algorithm (typically ``sort='quick'`` is in fact the fastest method; ``wendy`` also supports ``sort='merge'`` for a mergesort, ``sort='qsort'`` for the C standard library's own sorting algorithm, and ``sort='parallel'`` for a parallel implementation of mergesort)."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"x= copy.copy(x_init)\n",
"v= copy.copy(v_init)+sigma\n",
"g= wendy.nbody(x,v,m,0.05*tdyn,nleap=10,approx=True,sort='tim',ext_force=iso_force)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|██████████| 1000/1000 [01:41<00:00, 9.82it/s]\n"
]
}
],
"source": [
"nt= 1000\n",
"xt= numpy.empty((N,nt+1))\n",
"vt= numpy.empty((N,nt+1))\n",
"xt[:,0]= x\n",
"vt[:,0]= v\n",
"for ii in tqdm.trange(nt):\n",
" tx,tv= next(g)\n",
" xt[:,ii+1]= tx\n",
" vt[:,ii+1]= tv"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we plot the evolution of the phase-space distribution of time, color-coding the points by their initial energy in the unperturbed gravitational field, and we see that a strong spiral quickly develops and winds up over time. The spiral develops, because the starts of different energies have different frequencies and they therefore orbit on different times. The frequency goes down with increasing energy (or the period goes up), so the result is a winding spiral pattern:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 9,
"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(-4.99,4.99)\n",
" return (line1[0],)\n",
"figsize(6,4)\n",
"fig, ax= subplots()\n",
"# Directly compute the initial energy from the known sech^2 disk potential\n",
"c= v_init**2./2.+2.*sigma2*numpy.log(numpy.cosh(0.5*x_init/zh))\n",
"s= 5.*((c-numpy.amin(c))/(numpy.amax(c)-numpy.amin(c))*2.+1.)\n",
"line= ax.scatter(x,v,c=c,s=s,edgecolors='None',cmap=cm.jet_r)\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[:,ii*subsamp],vt[:,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('phasespiral_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": {},
"source": [
"When we look at the evolution of the density, the phase-space spiral shows up as a wave pattern that propagates up and down, as also seen in the data ([Bennett & Bovy 2019](https://ui.adsabs.harvard.edu/abs/2019MNRAS.482.1417B/abstract)):"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 10,
"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('phasespiral_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": "markdown",
"metadata": {},
"source": [
"Of course, this is a highly simplified model and it does not capture important aspects of the *Gaia* phase-space spiral, such as it's more obvious appearance when color-coding stars by their disk-planar rotational and radial velocities, which is due to the coupling between the vertical and planar motion of stars near the Sun. But this simple model illustrates that a spiral in the local vertical phase-space density naturally develops if something perturbs the phase-space distribution (for example, a passing Milky Way satellite galaxy)."
]
},
{
"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.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}