Notebook

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¶

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):
    # Add mass
    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:
    anim.save('adiabcontract_phasespace.gif',writer='imagemagick',dpi=80)
# The following is necessary to just get the movie, and not an additional initial frame
plt.close()
out= HTML(anim.to_html5_video())
plt.close()
out

Out[5]:

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.

In [6]:

figsize(6,4)
fig, ax= subplots()
ii= 0
a= ax.hist(xt[:,ii],bins=31,histtype='step',lw=1.,color='k',range=[-8.,8.],weights=31./16./N*numpy.ones(N))
xs= numpy.linspace(-8.,8.,101)
ax.plot(xs,totmass/4./zh/numpy.cosh(xs/2./zh)**2.,'b--',lw=2.,zorder=0)
ax.set_xlim(-8.,8.)
ax.set_ylim(10.**-3.,1.)
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$\rho(x)$')
ax.set_yscale('log')
ax.annotate(r'$t=0$',(0.95,0.95),xycoords='axes fraction',
             horizontalalignment='right',verticalalignment='top',size=18.)
subsamp= 4
def animate(ii):
    ax.clear()
    a= ax.hist(xt[:,ii*subsamp],bins=31,histtype='step',lw=1.,color='k',range=[-8.,8.],weights=31./16./N*numpy.ones(N))
    xs= numpy.linspace(-8.,8.,101)
    ax.plot(xs,totmass/4./zh/numpy.cosh(xs/2./zh)**2.,'b--',lw=2.,zorder=0)
    ax.set_xlim(-8.,8.)
    ax.set_ylim(10.**-3.,1.)
    ax.set_xlabel(r'$x$')
    ax.set_ylabel(r'$\rho(x)$')
    ax.set_yscale('log')
    ax.annotate(r'$t=%.0f$' % (ii*subsamp/20.),
                (0.95,0.95),xycoords='axes fraction',
                horizontalalignment='right',verticalalignment='top',size=18.)
    return a[2]
anim = animation.FuncAnimation(fig,animate,#init_func=init_anim_frame,
                               frames=nt//subsamp,interval=40,blit=True,repeat=True)
if _SAVE_GIFS:
    anim.save('adiabcontract_density.gif',writer='imagemagick',dpi=80)
# The following is necessary to just get the movie, and not an additional initial frame
plt.close()
out= HTML(anim.to_html5_video())
plt.close()
out

Out[6]:

In [ ]: