In [1]:
import numpy
import tqdm
import wendy
%pylab inline
import matplotlib.animation as animation
from matplotlib import cm
from IPython.display import HTML
_SAVE_GIFS= False
rcParams.update({'axes.labelsize': 17.,
'font.size': 12.,
'legend.fontsize': 17.,
'xtick.labelsize':15.,
'ytick.labelsize':15.,
'text.usetex': _SAVE_GIFS,
'figure.figsize': [5,5],
'xtick.major.size' : 4,
'ytick.major.size' : 4,
'xtick.minor.size' : 2,
'ytick.minor.size' : 2,
'legend.numpoints':1})
numpy.random.seed(2)
Populating the interactive namespace from numpy and matplotlib

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:

In [2]:
# Initial disk
N= 10000
# compute zh based on sigma and totmass
totmass= 1. # Sigma above
sigma= 1.
zh= sigma**2./totmass # twopiG = 1. in our units
tdyn= zh/sigma
x= numpy.arctanh(2.*numpy.random.uniform(size=N)-1)*zh*2.
v= numpy.random.normal(size=N)*sigma
v-= numpy.mean(v) # stabilize
m= totmass*numpy.ones_like(x)/N
In [3]:
# Properties of colder disk grown within
N_cold= 10000
totmass_cold= 1. # Sigma above
sigma_cold= 0.1
zh_cold= sigma_cold**2./totmass_cold # twopiG = 1. in our units
x_cold= numpy.arctanh(2.*numpy.random.uniform(size=N_cold)-1)*zh_cold*2.
v_cold= numpy.random.normal(size=N_cold)*sigma_cold
v_cold-= numpy.mean(v_cold) # stabilize
m_cold= totmass_cold*numpy.ones_like(x_cold)/N_cold

Now we run the $N$-body integration, adding a small fraction of the cold disk in each time step:

In [4]:
nt= 2500
xt= numpy.empty((N,nt+1))
vt= numpy.empty((N,nt+1))
xt[:,0]= x
vt[:,0]= v
xt_cold= numpy.zeros((N_cold,nt+1))
vt_cold= numpy.zeros((N_cold,nt+1))
xt_cold[:,0]= x_cold
vt_cold[:,0]= v_cold
for ii in tqdm.trange(nt):
skip= N_cold//nt
tx= numpy.concatenate((x,x_cold[:(ii+1)*skip]))
tv= numpy.concatenate((v,v_cold[:(ii+1)*skip]))
tm= numpy.concatenate((m,m_cold[:(ii+1)*skip]))
g= wendy.nbody(tx,tv,tm,0.05*tdyn,approx=True,nleap=10)
tx,tv= next(g)
xt[:,ii+1]= tx[:N]
vt[:,ii+1]= tv[:N]
xt_cold[:(ii+1)*skip,ii+1]= tx[N:N+(ii+1)*skip]
vt_cold[:(ii+1)*skip,ii+1]= tv[N:N+(ii+1)*skip]
# Re-set x,v
x= tx[:N]
v= tv[:N]
x_cold[:(ii+1)*skip]= tx[N:N+(ii+1)*skip]
v_cold[:(ii+1)*skip]= tv[N:N+(ii+1)*skip]
100%|██████████| 2500/2500 [01:28<00:00, 23.43it/s]

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:

In [5]:
def init_anim_frame():
line1= plot([],[])
xlabel(r'$x$')
ylabel(r'$v$')
xlim(-7.99,7.99)
ylim(-5.99,5.99)
return (line1[0],)
figsize(6,4)
fig, ax= subplots()
c= wendy.energy(xt[:,0],vt[:,0],m,individual=True)
s= 5.*((c-numpy.amin(c))/(numpy.amax(c)-numpy.amin(c))*2.+1.)
line= ax.scatter(x[::N//5000],v[::N//5000],c=c[::N//5000],s=s[::N//1000],edgecolors='None',cmap=cm.jet)
txt= ax.annotate(r'$t=%.0f$' % (0.),
(0.95,0.95),xycoords='axes fraction',
horizontalalignment='right',verticalalignment='top',size=18.)
subsamp= 4
def animate(ii):
line.set_offsets(numpy.array([xt[::N//5000,ii*subsamp],vt[::N//5000,ii*subsamp]]).T)
txt.set_text(r'$t=%.0f$' % (ii*subsamp/20.))
return (line,)
anim = animation.FuncAnimation(fig,animate,init_func=init_anim_frame,
frames=nt//subsamp,interval=40,blit=True,repeat=True)
if _SAVE_GIFS: