In [1]:

import numpy
from galpy.potential import LogarithmicHaloPotential
from galpy.orbit import Orbit
from galpy.util import bovy_plot, bovy_coords, bovy_conversion
%pylab inline




Populating the interactive namespace from numpy and matplotlib



# Initial conditions for $N$-body simulations to create the impact we want

Setup the potential and coordinate system



In [2]:

lp= LogarithmicHaloPotential(normalize=1.,q=0.9)
R0, V0= 8., 220.



Functions for converting coordinates between rectangular to cylindrical:



In [3]:

def rectangular_to_cylindrical(xv):
R,phi,Z= bovy_coords.rect_to_cyl(xv[:,0],xv[:,1],xv[:,2])
vR,vT,vZ= bovy_coords.rect_to_cyl_vec(xv[:,3],xv[:,4],xv[:,5],R,phi,Z,cyl=True)
out= numpy.empty_like(xv)
# Preferred galpy arrangement of cylindrical coordinates
out[:,0]= R
out[:,1]= vR
out[:,2]= vT
out[:,3]= Z
out[:,4]= vZ
out[:,5]= phi
return out
def cylindrical_to_rectangular(xv):
# Using preferred galpy arrangement of cylindrical coordinates
X,Y,Z= bovy_coords.cyl_to_rect(xv[:,0],xv[:,5],xv[:,3])
vX,vY,vZ= bovy_coords.cyl_to_rectvec(xv[:,1],xv[:,2],xv[:,4],xv[:,5])
out= numpy.empty_like(xv)
out[:,0]= X
out[:,1]= Y
out[:,2]= Z
out[:,3]= vX
out[:,4]= vY
out[:,5]= vZ
return out



At the time of impact, the phase-space coordinates of the GC can be computed using orbit integration:



In [4]:

xv_prog_init= numpy.array([30.,0.,0.,0.,105.74895,105.74895])
RvR_prog_init= rectangular_to_cylindrical(xv_prog_init[:,numpy.newaxis].T)[0,:]
prog_init= Orbit([RvR_prog_init[0]/R0,RvR_prog_init[1]/V0,RvR_prog_init[2]/V0,
RvR_prog_init[3]/R0,RvR_prog_init[4]/V0,RvR_prog_init[5]],ro=R0,vo=V0)
times= numpy.linspace(0.,10./bovy_conversion.time_in_Gyr(V0,R0),10001)
prog_init.integrate(times,lp)
xv_prog_impact= [prog_init.x(times[-1]),prog_init.y(times[-1]),prog_init.z(times[-1]),
prog_init.vx(times[-1]),prog_init.vy(times[-1]),prog_init.vz(times[-1])]




galpyWarning: Using C implementation to integrate orbits
galpyWarning: The behavior of Orbit member functions has changed in versions > 0.1 to return positions in kpc, velocities in km/s, energies and the Jacobi integral in (km/s)^2, the angular momentum o.L() and actions in km/s kpc, frequencies in 1/Gyr, and times and periods in Gyr if a distance and velocity scale was specified upon Orbit initialization with ro=...,vo=...; you can turn this off by specifying use_physical=False when calling the function (e.g., o=Orbit(...); o.R(use_physical=False)



The DM halo at the time of impact is at the following location:



In [5]:

xv_dm_impact= numpy.array([-13.500000,2.840000,-1.840000,6.82200571,132.7700529,149.4174464])
RvR_dm_impact= rectangular_to_cylindrical(xv_dm_impact[:,numpy.newaxis].T)[0,:]
dm_impact= Orbit([RvR_dm_impact[0]/R0,RvR_dm_impact[1]/V0,RvR_dm_impact[2]/V0,
RvR_dm_impact[3]/R0,RvR_dm_impact[4]/V0,RvR_dm_impact[5]],ro=R0,vo=V0)
dm_impact= dm_impact.flip()
times= numpy.linspace(0.,10./bovy_conversion.time_in_Gyr(V0,R0),1001)
dm_impact.integrate(times,lp)



The orbits over the past 10 Gyr for both objects are:



In [6]:

prog_init.plot()
dm_impact.plot(overplot=True)
plot(RvR_dm_impact[0],RvR_dm_impact[3],'ro')
xlim(0.,35.)
ylim(-20.,20.)




Out[6]:

(-20.0, 20.0)



## Initial condition for the King cluster

We start the King cluster at 10.25 WD time units, which corresponds to 10.25x0.9777922212082034 Gyr. The phase-space coordinates of the cluster are then:



In [7]:

prog_backward= prog_init.flip()
ts= numpy.linspace(0.,(10.25*0.9777922212082034-10.)/bovy_conversion.time_in_Gyr(V0,R0),1001)
prog_backward.integrate(ts,lp)
print [prog_backward.x(ts[-1]),prog_backward.y(ts[-1]),prog_backward.z(ts[-1]),
-prog_backward.vx(ts[-1]),-prog_backward.vy(ts[-1]),-prog_backward.vz(ts[-1])]




[29.577809361172694, -2.4079802843036058, -2.4053154109323742, 36.904368548683991, 104.25395168642832, 103.9042939586286]



## Starting 0.125 time units ago



In [8]:

dm_impact= Orbit([RvR_dm_impact[0]/R0,RvR_dm_impact[1]/V0,RvR_dm_impact[2]/V0,
RvR_dm_impact[3]/R0,RvR_dm_impact[4]/V0,RvR_dm_impact[5]],ro=R0,vo=V0)
dm_impact= dm_impact.flip()
ts= numpy.linspace(0.,0.125*0.9777922212082034/bovy_conversion.time_in_Gyr(V0,R0),10001)
dm_impact.integrate(ts,lp)
print [dm_impact.x(ts[-1]),dm_impact.y(ts[-1]),dm_impact.z(ts[-1]),
-dm_impact.vx(ts[-1]),-dm_impact.vy(ts[-1]),-dm_impact.vz(ts[-1])]




[5.2570183967301833, -7.8848620802774407, -4.5786210693408176, -195.34540991555195, -51.644977186288536, -155.93944808745755]



## Starting 0.25 time units ago



In [9]:

dm_impact= Orbit([RvR_dm_impact[0]/R0,RvR_dm_impact[1]/V0,RvR_dm_impact[2]/V0,
RvR_dm_impact[3]/R0,RvR_dm_impact[4]/V0,RvR_dm_impact[5]],ro=R0,vo=V0)
dm_impact= dm_impact.flip()
ts= numpy.linspace(0.,0.25*0.9777922212082034/bovy_conversion.time_in_Gyr(V0,R0),10001)
dm_impact.integrate(ts,lp)
print [dm_impact.x(ts[-1]),dm_impact.y(ts[-1]),dm_impact.z(ts[-1]),
-dm_impact.vx(ts[-1]),-dm_impact.vy(ts[-1]),-dm_impact.vz(ts[-1])]




[4.8965792014085778, 8.9181416272959186, 10.033551072654888, 149.50458227428786, -97.714578287010966, 21.235006395774576]



## Starting 0.375 time units ago



In [10]:

dm_impact= Orbit([RvR_dm_impact[0]/R0,RvR_dm_impact[1]/V0,RvR_dm_impact[2]/V0,
RvR_dm_impact[3]/R0,RvR_dm_impact[4]/V0,RvR_dm_impact[5]],ro=R0,vo=V0)
dm_impact= dm_impact.flip()
ts= numpy.linspace(0.,0.375*0.9777922212082034/bovy_conversion.time_in_Gyr(V0,R0),10001)
dm_impact.integrate(ts,lp)
print [dm_impact.x(ts[-1]),dm_impact.y(ts[-1]),dm_impact.z(ts[-1]),
-dm_impact.vx(ts[-1]),-dm_impact.vy(ts[-1]),-dm_impact.vz(ts[-1])]




[-8.5662119198502609, 0.50120970587644287, -6.4724441270571687, -50.169060756613071, 214.43730878238947, 112.96727941016003]



## Starting 0.50 time units ago



In [11]:

dm_impact= Orbit([RvR_dm_impact[0]/R0,RvR_dm_impact[1]/V0,RvR_dm_impact[2]/V0,
RvR_dm_impact[3]/R0,RvR_dm_impact[4]/V0,RvR_dm_impact[5]],ro=R0,vo=V0)
dm_impact= dm_impact.flip()
ts= numpy.linspace(0.,0.50*0.9777922212082034/bovy_conversion.time_in_Gyr(V0,R0),10001)
dm_impact.integrate(ts,lp)
print [dm_impact.x(ts[-1]),dm_impact.y(ts[-1]),dm_impact.z(ts[-1]),
-dm_impact.vx(ts[-1]),-dm_impact.vy(ts[-1]),-dm_impact.vz(ts[-1])]




[7.9948807918667892, -11.344426189362455, 0.9626499711978439, -100.4814119301404, -84.037056563181252, -151.74007155179646]



## Initial conditions for the Plummer DM subhalo with $\lambda$ scaled interaction velocities

To test the impulse approximation, we want to simulate interactions where the relative velocity ${\bf w}$ is changed by a factor of $\lambda$: ${\bf w} \rightarrow \lambda {\bf w}$. We start by computing the relative velocity for the impacts above and define a function that returns a dark-matter velocity after scaling the relative velocity by $\lambda$:



In [12]:

v_gc= numpy.array([xv_prog_impact[3],xv_prog_impact[4],xv_prog_impact[5]])
v_dm= numpy.array([6.82200571,132.7700529,149.4174464])
w_base= v_dm-v_gc
def v_dm_scaled(lam):
return w_base*lam+v_gc



## Starting 0.25 time units ago, scaled down by 0.5



In [13]:

lam= 0.5
xv_dm_impact= numpy.array([-13.500000,2.840000,-1.840000,v_dm_scaled(lam)[0],v_dm_scaled(lam)[1],v_dm_scaled(lam)[2]])
RvR_dm_impact= rectangular_to_cylindrical(xv_dm_impact[:,numpy.newaxis].T)[0,:]
dm_impact= Orbit([RvR_dm_impact[0]/R0,RvR_dm_impact[1]/V0,RvR_dm_impact[2]/V0,
RvR_dm_impact[3]/R0,RvR_dm_impact[4]/V0,RvR_dm_impact[5]],ro=R0,vo=V0)
dm_impact= dm_impact.flip()
ts= numpy.linspace(0.,0.25*0.9777922212082034/bovy_conversion.time_in_Gyr(V0,R0),10001)
dm_impact.integrate(ts,lp)
print [dm_impact.x(ts[-1]),dm_impact.y(ts[-1]),dm_impact.z(ts[-1]),
-dm_impact.vx(ts[-1]),-dm_impact.vy(ts[-1]),-dm_impact.vz(ts[-1])]




[0.90421868718568454, 3.7176316528507094, 11.079300548190444, 71.494711851673969, 1.0444192199789275, 136.43737446859959]




In [ ]: