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]: