The Kinematics of the Local Group in a Cosmological Context

This notebook containts all data analysis for the ApJ Letter The Kinematics of the Local Group in a Cosmological Context by Jaime E. Forero-Romero, Yehuda Hoffman, Sebastian Bustamante, Stefan Gottloeber and Gustavo Yepes.

Follow through the notebook to generate the figures and numbers used in the paper . You can also generate results that were mentioned but not explicitly reported. For instance, the results for pairs obtained for a Friend-of-Friends halo finder.

ar$\chi$iv: 1303.2690v1

See paper here: The kinematics of the Local Group in a cosmological context - ArXiv link

And github repo here: Kinematics github.

The following has been modified to highlight the possibilities of analysis for the Harley Wood Winter School talk on 2014-07-19 by Jonathan Whitmore.

In [1]:
%pylab inline
%load_ext autoreload
%autoreload 2
from make2DHistogram import *
from generateInfo import *
Populating the interactive namespace from numpy and matplotlib
In [2]:
# Global paths and constants
data_path = "../LG_Kinematics/data/"

G_GRAV = 4.54E-48 #units of Mpc^3 Msun^-1 s^-2
KM_TO_MPC = 3.2E-20
# HUBBLE = 0.70 # Original March 13, 2013
# Planck 2013 http://arxiv.org/abs/1303.5076  20 Mar 2013
# H0=67.3+/-1.2 km/s/Mpc
HUBBLE = 0.673
E_UNITS = 1.0E-36
L_UNITS = 1.0

Here one can select between two halo finders: "BDM" or "FOF".

The boolean narrow_data is True if we want to produced plots with a narrower selection for distances and total mass closer to the observations. narrow_data = False corresponds to the full sample and narrow_data = True corresponds to the reduced sample.

In [3]:
halo_finder="BDM" 
narrow_data = False

MonteCarlo for the observed values

Here we generate one million points from the observational constraints for the LG (separation, velocities, distances)

In [4]:
# The reference is van der Marel et al. ApJ,753:8
V_M31 = array([66.1,-76.3,45.1]) #km/s
V_M31_sigma = array([26.7, 19.0, 26.5])

R_M31 = array([-378.9,612.7,-283.1]) # kpc
R_M31_sigma = abs(R_M31/20.0) #This is an approximated value, not the actual observational uncertainty

print R_M31_sigma, sqrt(sum(R_M31*R_M31))

R_M31 = R_M31 / 1000.0 # in Mpc
R_M31_sigma = R_M31_sigma /1000.0

M31_Mass = 1.6E12
M31_Mass_sigma = 0.5E12
MW_Mass = 1.6E12
MW_Mass_sigma = 0.5E12

Mass_Total = 3.14E12 # in units of Msun
Mass_Total_sigma = 0.58E12
[ 18.945  30.635  14.155] 774.023326522
In [5]:
#generate the velocities, positions and masses
n_points = 1000000
V_M31_MonteCarlo = zeros((n_points,3))
R_M31_MonteCarlo = zeros((n_points,3))

M31_Mass_MonteCarlo = (random.normal(M31_Mass,M31_Mass_sigma, n_points))
MW_Mass_MonteCarlo = (random.normal(MW_Mass,MW_Mass_sigma, n_points))

negative_M31 = where(M31_Mass_MonteCarlo < 0.0)
negative_MW = where(MW_Mass_MonteCarlo < 0.0)

Total_Mass_MonteCarlo = M31_Mass_MonteCarlo + MW_Mass_MonteCarlo
Reduced_Mass_MonteCarlo = (M31_Mass_MonteCarlo * MW_Mass_MonteCarlo)/Total_Mass_MonteCarlo

for i in range(3):
    V_M31_MonteCarlo[:,i] = random.normal(V_M31[i],V_M31_sigma[i],n_points)
    R_M31_MonteCarlo[:,i] = random.normal(R_M31[i],R_M31_sigma[i],n_points)

#don't use the items with negative values of the mass
positive = where((M31_Mass_MonteCarlo > 0.0) & (MW_Mass_MonteCarlo > 0.0))[0]

V_M31_MonteCarlo = V_M31_MonteCarlo[positive,:]
R_M31_MonteCarlo = R_M31_MonteCarlo[positive,:]
M31_Mass_MonteCarlo = M31_Mass_MonteCarlo[positive]
MW_Mass_MonteCarlo = MW_Mass_MonteCarlo[positive]

Total_Mass_MonteCarlo = M31_Mass_MonteCarlo + MW_Mass_MonteCarlo
Reduced_Mass_MonteCarlo = (M31_Mass_MonteCarlo * MW_Mass_MonteCarlo)/Total_Mass_MonteCarlo
    
In [6]:
# Compute the angular momentum, energy and lambda

# array of angular momentum (vector) and the norm
J = cross(R_M31_MonteCarlo, V_M31_MonteCarlo) 
J_norm = sqrt(J[:,0]**2 + J[:,1]**2 + J[:,2]**2) 

# Mechanical energy
R_M31_MonteCarlo_norm  = sqrt(R_M31_MonteCarlo[:,0]**2 + R_M31_MonteCarlo[:,1]**2 + R_M31_MonteCarlo[:,2]**2)
E_kin = 0.5*(V_M31_MonteCarlo[:,0]**2 + V_M31_MonteCarlo[:,1]**2 + V_M31_MonteCarlo[:,2]**2)
E_kin = E_kin * KM_TO_MPC * KM_TO_MPC
E_pot = -G_GRAV * Total_Mass_MonteCarlo / R_M31_MonteCarlo_norm

E_kin = E_kin/E_UNITS
E_pot = E_pot/E_UNITS
E = E_kin + E_pot 

# Lambda
lambda_obs = (Reduced_Mass_MonteCarlo**1.5) * (J_norm * L_UNITS * KM_TO_MPC)* sqrt(abs(E * E_UNITS)) / (G_GRAV * (Total_Mass_MonteCarlo)**(5.0/2.0))
In [7]:
# Total mass histogram
n, bins, patches = hist(Total_Mass_MonteCarlo, 20, normed=0, color = 'silver')
In [8]:
# M31 mass histogram
n, bins, patches = hist(M31_Mass_MonteCarlo, 20, normed=0, color = 'silver')
In [9]:
#surface density plots in the E-J space. This defines the 1-sigma surface
H, xedges,yedges = histogram2d(J_norm,E, bins=(30,38))
X,Y=meshgrid(yedges[0:-1],xedges[0:-1])
th=4000 # 4000, corresponds to 1-sigma, 200 to 3-sigma
print 'This is the fraction of points included inside the contour:',sum(H[where(H>th)])/sum(H) 
levels=array([th])
C=contour(X,Y,H, levels, colors='k', linewidths=2)
ax = axes()
ax.set_xlabel("Energy")
ax.set_ylabel("Angular Momentum")
This is the fraction of points included inside the contour: 0.718192396819
Out[9]:
<matplotlib.text.Text at 0x10893a7d0>

N-body data

The next four cells put into a single file the data from the pairs (their IDs) and the halo catalogs. This is done for Bolshoi and the three constrained simulation

In [10]:
# Generate the pair data for Bolshoi
FOF_file=data_path+"Bolshoi_catalog_"+halo_finder+".dat"
ID_file=data_path+"IsoPairs_Bolshoi_"+halo_finder+"_ID.dat"
output_file=data_path+"pair_physical_vweb_Bolshoi"+halo_finder+".dat"
generate_info(FOF_file, ID_file, output_file=output_file, hubble=HUBBLE)
1924
selected in total 1823
In [11]:
# Generate the pair data for CLUES 1
FOF_file=data_path+"IC2710_halos_catalog.dat"
ID_file=data_path+"index_LG_IC2710"
output_file=data_path+"pair_physical_vweb_IC2710.dat"
generate_info(FOF_file, ID_file, output_file=output_file, hubble=0.73)
# generate_info(FOF_file, ID_file, output_file=output_file, hubble=HUBBLE)
1
selected in total 1
In [12]:
# Generate the pair data for CLUES 2
FOF_file=data_path+"IC16953_halos_catalog.dat"
ID_file=data_path+"index_LG_IC16953"
output_file=data_path+"pair_physical_vweb_IC16953.dat"
generate_info(FOF_file, ID_file, output_file=output_file, hubble=0.73)
# generate_info(FOF_file, ID_file, output_file=output_file, hubble=HUBBLE)
1
selected in total 1
In [13]:
# Generate the pair data for CLUES 3
FOF_file=data_path+"IC10909_halos_catalog.dat"
ID_file=data_path+"index_LG_IC10909"
output_file=data_path+"pair_physical_vweb_IC10909.dat"
generate_info(FOF_file, ID_file, output_file=output_file, hubble=0.73)
# generate_info(FOF_file, ID_file, output_file=output_file, hubble=HUBBLE)
1
selected in total 1

After saving all the data into a single file, that file is reloaded to keep in memory the data we care about (radial velocities, tangential velocities, separations, total_mass, angular momentum, energy, lambda). This is for Bolshoi.

In [14]:
#load the data from Bolshoi and make selections on distances and mass
bolshoi_file=data_path+"pair_physical_vweb_Bolshoi"+halo_finder+".dat"
physical_prop=np.loadtxt(bolshoi_file) 

#define the main quantities we will use later to prepare the plots
radial_vel = physical_prop[:,0]
tangential_vel = physical_prop[:,1]
distance = physical_prop[:,8]
total_mass = physical_prop[:,9]/ 1E12

energy = physical_prop[:,5] + physical_prop[:,6] 
angular_momentum = physical_prop[:,7] 

mu = physical_prop[:,10]

lambda_orbit = mu**1.5 * (angular_momentum * L_UNITS * KM_TO_MPC)* sqrt(abs(energy * E_UNITS)) / (G_GRAV * (total_mass *1E12 )**(5.0/2.0))

if(narrow_data):
    index_good_mass = where((total_mass>1.0) & (total_mass<4.0) & (distance>0.70) & (distance<0.80))
    index_good_mass = index_good_mass[0]
    physical_prop = physical_prop[index_good_mass,:]

    radial_vel = physical_prop[:,0]
    tangential_vel = physical_prop[:,1]
    distance = physical_prop[:,8]
    total_mass = physical_prop[:,9]/1E12

    energy = physical_prop[:,5] + physical_prop[:,6] 
    angular_momentum = physical_prop[:,7] 

    mu = physical_prop[:,10]

    lambda_orbit = mu**1.5 * (angular_momentum * L_UNITS * KM_TO_MPC)* sqrt(abs(energy * E_UNITS)) / (G_GRAV * (total_mass *1E12 )**(5.0/2.0))
In [15]:
#Figure with the separation distribution
fig = figure(1, figsize=(9.5,6.5))
ax = axes()
ax.set_xlabel("$\mathrm{Pair\ Separation\ [Mpc]}$",fontsize=25)
ax.set_ylabel("$\mathrm{N_{pairs}}$",fontsize=25)
ticklabels_x = ax.get_xticklabels()
ticklabels_y = ax.get_yticklabels()

for label_x in ticklabels_x:
    label_x.set_fontsize(18)
    label_x.set_family('serif')
for label_y in ticklabels_y:
    label_y.set_fontsize(18)
    label_y.set_family('serif')
n, bins, patches = ax.hist(distance, 20, normed=0, color = 'silver')
ax.set_xlim([0.4,1.0])
filename='separation_'+halo_finder
# if(narrow_data):
#     filename=filename+'_narrow'
# savefig(filename + '.eps',format = 'eps', transparent=True)
# savefig(filename + '.pdf',format = 'pdf', transparent=True)
# savefig(filename + '.png',format = 'png', transparent=True)
In [16]:
#Figure with the total mass distribution
print amin(total_mass), amax(total_mass)

fig = figure(1, figsize=(9.5,6.5))
ax = axes()
ax.set_xlabel("$\mathrm{Total\ Mass\ [10^{12}\ M_{\odot}]}$",fontsize=25)
ax.set_ylabel("$\mathrm{N_{pairs}}$",fontsize=25)
ticklabels_x = ax.get_xticklabels()
ticklabels_y = ax.get_yticklabels()

for label_x in ticklabels_x:
    label_x.set_fontsize(18)
    label_x.set_family('serif')
for label_y in ticklabels_y:
    label_y.set_fontsize(18)
    label_y.set_family('serif')

ax.set_xlim([1.0,12.0])
n, bins, patches = ax.hist(total_mass, 20, normed=0, color = 'silver')

filename='total_mass_'+halo_finder
if(narrow_data):
    filename=filename+'_narrow'
savefig(filename + '.eps',format = 'eps', transparent=True)
savefig(filename + '.pdf',format = 'pdf', transparent=True)
savefig(filename + '.png',format = 'png', transparent=True)
1.50099554235 13.844576523
In [17]:
#observational data
mean_point = [-110.0, 0.0]
sigma_point = [4.4, 34.3]
delta_radial = 4.5
delta_tangential = 17.5

#load the data from the CLUES pairs
n_clues = 3
clues1=data_path+"pair_physical_vweb_IC10909.dat"
clues2=data_path+"pair_physical_vweb_IC16953.dat"
clues3=data_path+"pair_physical_vweb_IC2710.dat"
clues=chararray([n_clues], itemsize=1024)
clues[0]=clues1;clues[1]=clues2; clues[2]=clues3
clues_data = np.zeros([n_clues,4])
if(narrow_data):
    n_clues = 1
for i in range(n_clues):
    tmp_dat = np.loadtxt(clues[i])
    clues_data[i,0] = tmp_dat[0]# radial 
    clues_data[i,1] = tmp_dat[1]# tangential
    clues_data[i,2] = tmp_dat[8]# separation
    clues_data[i,3] = tmp_dat[9]/1E12 #total mass
    
filename='test_rt_'+halo_finder

if(narrow_data):
    filename=filename+'_narrow'

n=MakeSimple2DPlot(radial_vel,tangential_vel,filename=filename, 
                      clues=clues_data,clues_label='$\mathrm{Constrained\ Sim.}$',
                      contours=0,xlims=[-240.0,0.0], ylims=[0.0,240.0],
                      nbins=13,nxbins=13,nybins=13,bw=1, 
                      xlabel='$\mathrm{Radial\ Velocity\ [km\ s^{-1}]}$', 
                      ylabel='$\mathrm{Tangential\ Velocity\ [km\ s^{-1}]}$', 
                      mean_obs=mean_point, sigma_obs=sigma_point)

In [18]:
#find the highest density peak
H2d, r_edges, t_edges = histogram2d(radial_vel, tangential_vel, bins=(12,12), range=[[-240.0,0.0],[0.0,240.0]])
max_radial_id = argmax(H2d, axis=0)
max_tangential_id = argmax(H2d, axis=1)
max_r_id = argmax(H2d[arange(12), max_tangential_id])
max_t_id = argmax(H2d[max_radial_id, arange(12)])

theory_point = array([r_edges[max_r_id]+10.0, t_edges[max_t_id]+10.0])
print 'total number of pairs', size(radial_vel)
print 'location of the highest density peak in Bolshoi', theory_point

#selection of the number of galaxies around the observational constraints
index_obs = where((radial_vel>(mean_point[0]-sigma_point[0])) & (radial_vel<(mean_point[0]+sigma_point[0])) 
    & (tangential_vel<sigma_point[1]))
print 'number of points in observations', size(index_obs)

#selection of the number of galaxies around the Bolshoi point
index_max = where((radial_vel>(theory_point[0]-sigma_point[0])) & (radial_vel<(theory_point[0]+sigma_point[0])) 
    & (tangential_vel<(theory_point[1]+sigma_point[1]/2.0)) & (tangential_vel>(theory_point[1]-sigma_point[1]/2.0)))
print 'number of points around the highest density peak', size(index_max)


#find the number of points with radial velocities consistent with observations
tang_to_radial = abs(tangential_vel / radial_vel)
index = where(tang_to_radial<0.32)
print 'number of points with the right radial to tang', size(index), '. In percentage:', 1.0*size(index)/radial_vel.size
total number of pairs 1823
location of the highest density peak in Bolshoi [-70.  50.]
number of points in observations 5
number of points around the highest density peak 21
number of points with the right radial to tang 236 . In percentage: 0.129456939111
In [19]:
#Figure in the E-J plane
n_clues=3
clues1=data_path+"pair_physical_vweb_IC10909.dat"
clues2=data_path+"pair_physical_vweb_IC16953.dat"
clues3=data_path+"pair_physical_vweb_IC2710.dat"
clues=chararray([n_clues], itemsize=1024)
clues[0]=clues1;clues[1]=clues2; clues[2]=clues3
clues_data = np.zeros([n_clues,2])
if(narrow_data):
    n_clues = 1
for i in range(n_clues):
    tmp_dat = np.loadtxt(clues[i])
    clues_data[i,0] = tmp_dat[5] + tmp_dat[6] # energy
    clues_data[i,1] = tmp_dat[7] #angular momentum

filename='test_EJ_'+halo_finder
if(narrow_data):
    filename=filename+'_narrow'
MakeSimple2DPlot(energy,angular_momentum,filename=filename, 
                      clues=clues_data, clues_label='$\mathrm{Constrained\ Sim.}$', xlims=[-25.0,0.0], ylims=[0.0,150.0],
                      nbins=13,nxbins=13,nybins=13,bw=1, 
                      xlabel='$\mathrm{e_{tot} [10^{-36}} \mathrm{Mpc^2\ s^{-2}]}$', 
                      ylabel='$\mathrm{l_{orb} [Mpc\ } \mathrm{km\ s^{-1}]}$', 
                      contours=0, X_field=X, Y_field=Y, Z_field=H, levels=levels)
HOLA

Out[19]:
0
In [20]:
#Figure in the E-J plane
# Previous plot with 12 bins rather than 13!
n_clues=3
jw_bin_number = 12
clues1=data_path+"pair_physical_vweb_IC10909.dat"
clues2=data_path+"pair_physical_vweb_IC16953.dat"
clues3=data_path+"pair_physical_vweb_IC2710.dat"
clues=chararray([n_clues], itemsize=1024)
clues[0]=clues1;clues[1]=clues2; clues[2]=clues3
clues_data = np.zeros([n_clues,2])
if(narrow_data):
    n_clues = 1
for i in range(n_clues):
    tmp_dat = np.loadtxt(clues[i])
    clues_data[i,0] = tmp_dat[5] + tmp_dat[6] # energy
    clues_data[i,1] = tmp_dat[7] #angular momentum

filename='test_EJ_'+halo_finder
if(narrow_data):
    filename=filename+'_narrow'
MakeSimple2DPlot(energy,angular_momentum,filename=filename, 
                      clues=clues_data, clues_label='$\mathrm{Constrained\ Sim.}$', xlims=[-25.0,0.0], ylims=[0.0,150.0],
                      nbins=jw_bin_number,nxbins=jw_bin_number,nybins=jw_bin_number,bw=1, 
                      xlabel='$\mathrm{e_{tot} [10^{-36}} \mathrm{Mpc^2\ s^{-2}]}$', 
                      ylabel='$\mathrm{l_{orb} [Mpc\ } \mathrm{km\ s^{-1}]}$', 
                      contours=0, X_field=X, Y_field=Y, Z_field=H, levels=levels)
HOLA

Out[20]:
0
In [21]:
index = where((angular_momentum<40.0) & (energy>-20.0) & (energy<-5.0))
print 'approx. number of points inside contour', size(index), 'in percentage', 1.0 * size(index)/energy.size
approx. number of points inside contour 276 in percentage 0.151398793198
In [22]:
#Figure with the spin parameter distribution
fig = figure(1, figsize=(9.5,6.5))
ax = axes()
ax.set_xlabel("$\log_{10}\mathrm{\lambda}$",fontsize=25)
ax.set_ylabel("$\mathrm{Normalized\ Histogram}$",fontsize=25)
ticklabels_x = ax.get_xticklabels()
ticklabels_y = ax.get_yticklabels()

for label_x in ticklabels_x:
    label_x.set_fontsize(18)
    label_x.set_family('serif')
for label_y in ticklabels_y:
    label_y.set_fontsize(18)
    label_y.set_family('serif')

n_clues=3
clues1=data_path+"pair_physical_vweb_IC10909.dat"
clues2=data_path+"pair_physical_vweb_IC16953.dat"
clues3=data_path+"pair_physical_vweb_IC2710.dat"
clues=chararray([n_clues], itemsize=1024)
clues[0]=clues1;clues[1]=clues2; clues[2]=clues3
clues_data = np.zeros([n_clues,4])

if(narrow_data):
    n_clues = 1
for i in range(n_clues):
    tmp_dat = np.loadtxt(clues[i])
    clues_data[i,0] = tmp_dat[5] + tmp_dat[6] # energy
    clues_data[i,1] = tmp_dat[7] #angular momentum
    clues_data[i,2] = tmp_dat[9] # total mass
    clues_data[i,3] = tmp_dat[10] # reduced_mass
    
    e_CLUES = clues_data[i,0]
    l_CLUES = clues_data[i,1]
    total_mass_CLUES = clues_data[i,2]
    mu_CLUES = clues_data[i,3]

    lambda_CLUES = mu_CLUES**1.5 * (l_CLUES * L_UNITS * KM_TO_MPC) * \
        sqrt(abs(e_CLUES * E_UNITS)) / (G_GRAV * (total_mass_CLUES )**(5.0/2.0))
    lambda_CLUES = log10(lambda_CLUES)
    y_value=1.9

    ax.axvline(x=lambda_CLUES,ymin=0.8,ymax=0.95,c="black",linewidth=2)
    if(i==0):
        ax.scatter(lambda_CLUES, y_value, color='white', edgecolors='black', s=200, label='$\mathrm{Constrained\ Sim.}$')  
    else:
        ax.scatter(lambda_CLUES, y_value, color='white', edgecolors='black', s=200)  
        
        

hist_A, bin_edges_A = histogram(log10(lambda_obs), bins=arange(-3.0,0.5,0.1), density=True)
hist_B, bin_edges_B = histogram(log10(lambda_orbit), bins=arange(-3.0,0.5,0.1), density=True)

plot(bin_edges_A[:-1],hist_A,color='black',ms=1,linewidth=3.0, label="$\mathrm{Obs.}$")
plot(bin_edges_B[:-1],hist_B,'r--',color='black',ms=2,linewidth=3.0, label="$\mathrm{Bolshoi}$")

ax.legend(loc=0, scatterpoints=1, prop={'size':20})

ax.set_xlim([-3.0,0.5])
ax.set_ylim([0.0,2.2])
filename='test_lambda_'+halo_finder
if(narrow_data):
    filename=filename+'_narrow'
savefig(filename + '.eps',format = 'eps', transparent=True)
savefig(filename + '.pdf',format = 'pdf', transparent=True)
savefig(filename + '.png',format = 'png', transparent=True)
In [23]:
mean_lambda_obs = mean(log10(lambda_obs))
variance_obs = var(log10(lambda_obs))
mean_lambda_sim = mean(log10(lambda_orbit))
variance_sim = var(log10(lambda_orbit))

print 'mean lambda, observations', mean_lambda_obs, variance_obs
print 'mean lambda, bolshoi', mean_lambda_sim, variance_sim
print 'lambda CLUES', lambda_CLUES

# calculate the number of points in the simulation that are consistent with mean value and sigma
index_sim = where((log10(lambda_orbit)<(mean_lambda_obs+variance_obs)) & (log10(lambda_orbit) > (mean_lambda_obs-variance_obs)))
print 'there are ', size(index_sim), 'pairs consistent with the observations'

# mean lambda, observations -1.72013855372 0.0797602386996
# mean lambda, bolshoi -1.470309254 0.137777720569
# lambda CLUES -1.72803661534
# there are  257 pairs consistent with the observations
mean lambda, observations -1.72061467416 0.0799156324454
mean lambda, bolshoi -1.47345497862 0.138217470615
lambda CLUES -1.72803661534
there are  247 pairs consistent with the observations