In [1]:
import numpy
import tqdm
import warnings
warnings.simplefilter('ignore',DeprecationWarning)
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})
import copy
numpy.random.seed(2)
Populating the interactive namespace from numpy and matplotlib

Adiabatic versus non-adiabatic changes to an exponential disk

Let's investigate how resilient an exponential disk ($\rho(x) \propto e^{-|x|/H}$) is against non-adiabatic changes to the potential. We will first setup an exponential disk in approximate equilibrium, let it evolve for a while to let it fully equilibrate, and then inject energy into it in both an adiabatic and non-adiabatic manner. First we setup the disk and evolve it:

In [2]:
# Initial disk
N= 30000
totmass= 1. # Sigma above
sigma= 1.
zh= sigma**2./totmass # twopiG = 1. in our units
tdyn= zh/sigma
x= -numpy.log(1.-numpy.random.uniform(size=N))*zh
x[numpy.random.uniform(size=N) < 0.5]*= -1.
v= numpy.random.normal(size=N)*sigma*numpy.sqrt(1.-0.75*numpy.exp(-numpy.fabs(x)/zh))
v-= numpy.mean(v) # stabilize
m= totmass*numpy.ones_like(x)/N
In [3]:
g= wendy.nbody(x,v,m,0.05*tdyn,approx=True,nleap=1)
In [4]:
nt= 3000
xt= numpy.empty((N,nt+1))
vt= numpy.empty((N,nt+1))
Et= numpy.empty((nt+1))
xt[:,0]= x
vt[:,0]= v
Et[0]= wendy.energy(x,v,m)
for ii in tqdm.trange(nt):
    tx,tv= next(g)
    xt[:,ii+1]= tx
    vt[:,ii+1]= tv
    Et[ii+1]= wendy.energy(tx,tv,m)
x_start= xt[:,-1]
v_start= vt[:,-1]
100%|██████████| 3000/3000 [00:28<00:00, 104.44it/s]

The following movie shows that the exponential disk density does not significantly evolve:

In [5]:
figsize(6,4)
fig, ax= subplots()
ii= 0
a= ax.hist(xt[:,ii],bins=101,histtype='step',lw=1.,color='k',range=[-8.,8.],weights=101./16./N*numpy.ones(N))
xs= numpy.linspace(-8.,8.,101)
ax.plot(xs,totmass/2./zh*numpy.exp(-numpy.fabs(xs)/2./zh)**2.,'b--',lw=2.,zorder=0)
ax.set_xlim(-8.,8.)
ax.set_ylim(10.**-3.,.65)
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= 5
def animate(ii):
    ax.clear()
    norm= 1./N
    a= ax.hist(xt[:,ii*subsamp],bins=101,histtype='step',lw=1.,color='k',range=[-8.,8.],weights=101./16.*norm*numpy.ones(N))
    xs= numpy.linspace(-8.,8.,101)
    ax.plot(xs,totmass/2./zh*numpy.exp(-numpy.fabs(xs)/2./zh)**2.,'b--',lw=2.,zorder=0)
    ax.set_xlim(-8.,8.)
    ax.set_ylim(10.**-3.,.65)
    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)
# 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]:

Now we inject energy into about 40% of the orbits in an adiabatic manner, by changing the velocities over many dynamical times:

In [6]:
nt= 3000
x= copy.deepcopy(x_start)
v= copy.deepcopy(v_start)
xt= numpy.empty((N,nt+1))
vt= numpy.empty((N,nt+1))
xt[:,0]= x
vt[:,0]= v
launch_v= 2.*sigma
frac_lost= 0.4
time_lost= 100
nlost= int(round(frac_lost*N/4./time_lost))
g= wendy.nbody(x,v,m,0.05*tdyn,approx=True,nleap=1)
for ii in tqdm.trange(nt):
    if (ii >= 300 and ii < 300+time_lost) or (ii >= 400 and ii < 400+time_lost) \
        or (ii >= 500 and ii < 500+time_lost) or (ii >= 600 and ii < 600+time_lost):
        idx= numpy.random.permutation(N)[:nlost]#numpy.argpartition(numpy.fabs(x),nlost)[:nlost]
        if numpy.random.uniform() < 0.5:
            v[idx[:nlost//2]]= launch_v
            v[idx[nlost//2:]]= -launch_v
        else:
            v[idx[:nlost//2]]= -launch_v
            v[idx[nlost//2:]]= launch_v
        g= wendy.nbody(x,v,m,0.05*tdyn,approx=True,nleap=1)
    x,v= next(g)
    if (ii >= 300 and ii < 300+time_lost) or (ii >= 400 and ii < 400+time_lost) \
        or (ii >= 500 and ii < 500+time_lost) or (ii >= 600 and ii < 600+time_lost):
        xt[:,ii+1]= x
        vt[:,ii+1]= v
        vt[idx,ii+1]= numpy.nan
    else:
        xt[:,ii+1]= x
        vt[:,ii+1]= v
100%|██████████| 3000/3000 [00:32<00:00, 91.25it/s] 

The density remains approximately exponential throughout the evolution (energy is injected between $t = 15$ and $t = 35$; note that we remove the particles with inflated energies from the histogram to focus on the response of the system):

In [7]:
figsize(6,4)
fig, ax= subplots()
ii= 0
a= ax.hist(xt[:,ii],bins=101,histtype='step',lw=1.,color='k',range=[-8.,8.],weights=101./16./N*numpy.ones(N))
xs= numpy.linspace(-8.,8.,101)
ax.plot(xs,totmass/2./zh*numpy.exp(-numpy.fabs(xs)/2./zh)**2.,'b--',lw=2.,zorder=0)
ax.set_xlim(-8.,8.)
ax.set_ylim(10.**-3.,.65)
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.)
ax.annotate(r'$\mathrm{adiabatic}$',(0.05,0.95),xycoords='axes fraction',
             horizontalalignment='left',verticalalignment='top',size=18.)
subsamp= 5
def animate(ii):
    ax.clear()
    if ii*subsamp < 300:
        norm= 1./N
    elif ii*subsamp < 700.:
        norm= 1./(N-(ii*subsamp-300.)*nlost)
    else:
        norm= 1./(N-400.*nlost)
    idx= True-numpy.any(numpy.isnan(vt[:,:ii*subsamp]),axis=1)
    a= ax.hist(xt[idx,ii*subsamp]-numpy.mean(xt[idx,ii*subsamp]),
               bins=101,histtype='step',lw=1.,color='k',range=[-8.,8.],
               weights=101./16.*norm*numpy.ones(N)[idx])
    xs= numpy.linspace(-8.,8.,101)
    ax.plot(xs,totmass/2./zh*numpy.exp(-numpy.fabs(xs)/2./zh)**2.,'b--',lw=2.,zorder=0)
    ax.set_xlim(-8.,8.)
    ax.set_ylim(10.**-3.,.65)
    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.)
    ax.annotate(r'$\mathrm{adiabatic}$',(0.05,0.95),xycoords='axes fraction',
                 horizontalalignment='left',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('expdiskadiab_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[7]: