This is a simple, pure-python notebook to reproduce the constraints on the Gaia BH1 orbit from joint fitting of RVs and the astrometric constraints.
In [1]:
import numpy as np
import matplotlib.pyplot as plt
import emcee
import corner
%matplotlib inline
In [2]:
# the correlation matrix from Gaia, and the vector of best-fit parameters.
# these can be retrieved from the Gaia archive; I just copy them here so
# this notebook can be run without an internet connection or TAP instance.
corr_vec_BH1 = np.array([-0.9385065 , -0.10964689, 0.09768217, 0.02599761, -0.09885585,
-0.04858593, 0.04631094, -0.18919176, 0.12997928, 0.08754877,
-0.7026825 , 0.7199984 , -0.1333831 , -0.16456777, 0.06737773,
0.976118 , -0.90516096, -0.0790209 , -0.0389514 , 0.14850111,
-0.64331377, 0.7340439 , -0.7446727 , -0.32354277, -0.02299099,
0.10600532, -0.12226893, 0.70452565, 0.31863543, -0.2744563 ,
0.2560132 , 0.1476755 , -0.23248082, -0.85903525, 0.25882402,
-0.35804698, 0.9656088 , -0.9219766 , -0.2000572 , 0.05323359,
0.00881674, -0.5612588 , 0.9152694 , 0.8657674 , 0.1358042 ,
-0.01977275, 0.06291951, -0.08401399, 0.08238114, -0.81962156,
-0.07304911, -0.14815669, 0.03363582, 0.13756153, 0.0689249 ,
-0.21668243, 0.16537377, -0.2642203 , -0.15781169, 0.33261746,
0.7964632 , -0.14971647, 0.4369154 , -0.98625296, -0.0418497 ,
-0.25472474])
# ra, dec, parallax, pmra, pmdec, a_thiele_innes, b_thiele_innes, f_thiele_innes, g_thiele_innes, ecc, period, t_peri
mu_vec_BH1 = np.array([262.171208162297, -0.5810920153840334, 2.0954515888832432, -7.702050443675317,
-25.85042074524512, -0.262289119199277, 2.9291159041806485, 1.5248071482113947,
0.5343685858549198, 0.48893589298452034, 185.7656578920495, -12.024680365644883])
err_vec_BH1 = np.array([0.49589708, 0.15092158, 0.01745574, 0.02040691, 0.02699432, 0.16984475,
0.17521776, 0.15387644, 0.54653114, 0.074341,0.30688563, 6.34896183])
# the RV data
jds_BH1 = np.array([2457881.2899, 2458574.3663, 2459767.6226, 2459791.9186, 2459795.6461, 2459796.4995,
2459798.8399, 2459805.5101, 2459808.7388224313, 2459813.6045, 2459814.58740,
2459815.5927, 2459817.5278, 2459818.5266, 2459818.78698893, 2459819.5543,
2459820.5465, 2459821.5669, 2459822.5745, 2459823.5422430662, 2459824.5305,
2459825.5361, 2459823.8525, 2459824.8516, 2459826.7920, 2459828.5677, 2459829.5373,
2459829.5768, 2459830.6452, 2459831.6223, 2459833.7523, 2459834.5509, 2459834.7691,
2459835.7678, 2459838.8082, 2459838.7208, 2459840.7729, 2459845.5069, 2459855.5012,
2459868.5128, 2459877.6978])
rvs_BH1 = np.array([20.0, 8.9, 63.8, 131.90, 141.4, 142.7, 140.6, 127.7, 118, 90.5, 86.1, 81.90, 74.5, 71.0,
67.8, 67.0, 64.0, 60.7, 57.8, 54.8, 52.1, 49.8, 53.76, 51.18, 46.59, 42.2, 42.1, 40.4,
38.5, 36.5, 33.23, 32.3, 31.74, 30.14, 26.35, 27.5, 24.20, 19.4, 14.2, 9.3, 10.5])
rv_errs_BH1 = np.array([4.1, 5.6, 3, 0.1, 3, 3, 4, 1.0, 4, 0.3, 0.3, 0.3, 0.3, 0.3, 4, 0.3, 0.3, 0.3, 0.3,
0.3, 0.3, 0.3, 0.1, 0.1, 0.1, 0.3, 3, 0.3, 0.3, 0.5, 0.1, 0.3, 0.1, 0.1, 0.1, 4,
0.1, 1.0, 1.0, 1.0, 1.5])
rv_names_BH1 = np.array(['LAMOST', 'LAMOST', 'MagE', 'HIRES', 'MagE', 'MagE', 'GMOS', 'XSHOOTER', 'GMOS', 'FEROS',
'FEROS', 'FEROS', 'FEROS', 'FEROS', 'GMOS', 'FEROS', 'FEROS', 'FEROS', 'FEROS', 'FEROS',
'FEROS', 'FEROS', 'HIRES', 'HIRES', 'HIRES', 'FEROS', 'MagE', 'FEROS', 'FEROS', 'FEROS',
'HIRES', 'FEROS', 'HIRES', 'HIRES', 'HIRES', 'GMOS', 'HIRES', 'XSHOOTER', 'XSHOOTER',
'XSHOOTER', 'ESI'])
# prior on the mass of the luminous star
m1_obs, m1_err = 0.95, 0.05
# the reference epoch of the gaia t_periastron; i.e., JD 2016.0
jd_ref = 2457389.0
In [3]:
# convert the correlation matrix to a covariance matrix
Nparam = 12
triangle_nums = (np.arange(1, Nparam)*(np.arange(1, Nparam) - 1)//2)
cov_mat_BH1 = []
count = 0
for i, sigx in enumerate(err_vec_BH1):
this_row = []
for j, sigy in enumerate(err_vec_BH1):
if i == j:
this_row.append(err_vec_BH1[i]**2)
elif j < i:
this_row.append(cov_mat_BH1[j][i])
else:
this_idx = triangle_nums[j-1] + i
this_row.append(corr_vec_BH1[this_idx]*sigx*sigy)
cov_mat_BH1.append(this_row)
cov_mat_BH1 = np.array(cov_mat_BH1)
In [4]:
# functions to predict RVs at a given time
def fsolve_newton(Mi, ecc, xtol = 1e-10):
'''
numerically solve with newton's method.
Mi: float, 2*pi/P*(t - T_p)
ecc: float, eccentricity
xtol: float, tolerance
'''
eps = 1
EE = Mi + ecc*np.sin(Mi) + ecc**2/2*np.sin(2*Mi) + ecc**3/8*(3*np.sin(3*Mi) - np.sin(Mi))
while eps > xtol:
EE0 = EE
EE = EE0 + (Mi + ecc*np.sin(EE0) - EE0)/(1 - ecc*np.cos(EE0))
eps = np.abs(EE0 - EE)
return EE
def solve_kep_eqn(t, T_p, P, ecc):
'''
Solve Keplers equation E - e*sin(E) = M for E
Here, M = 2*pi/P*(t - T_p)
t: array of times
T_p: float; periastron time
P: float, period
ecc: float, eccentricity
'''
M = 2*np.pi/P * (t - T_p)
E = np.zeros(t.shape)
for i,Mi in enumerate(M):
E[i] = fsolve_newton(Mi = Mi, ecc = ecc, xtol = 1e-10)
return E
def get_radial_velocities(t, P, T_p, ecc, K, omega, gamma):
'''
t: array of times
P: float, period
T_p: float, periastron time
ecc: float, eccentricity
K: float, RV semi-amplitude
omega: float, longitude of periastron (radians, not degrees)
gamma: float, center-of-mass RV
'''
E = solve_kep_eqn(t = t, T_p = T_p, P = P, ecc = ecc)
f = 2*np.arctan2(np.sqrt(1+ecc)*np.sin(E*.5),np.sqrt(1-ecc)*np.cos(E*.5))
vr = K*(np.cos(f + omega) + ecc*np.cos(omega)) + gamma
return vr
In [5]:
# likelihood for joint fit of astrometry and RVs
def get_a0_mas(period, m1, m2, parallax):
'''
predicts a0 (photocenter semi-major axis) in milliarcsec,
assuming a dark companion
period: float, days
m1: float, Msun, (mass of luminous star)
m2: float, Msun, (mass of companion)
parallax: float, milliarcsec
'''
G = 6.6743e-11 # SI units
Msun = 1.98840987069805e+30
AU = 1.4959787e+11
a_au = (((period*86400)**2 * G * (m1*Msun + m2*Msun)/(4*np.pi**2))**(1/3.))/AU
a_mas = a_au*parallax
q = m2/m1
a0_mas = a_mas*q/(1+q)
return a0_mas
def get_Kstar_kms(period, inc_deg, m1, m2, ecc):
'''
RV semi-amplitude of the luminous star
period: float, days
inc_deg: float, inclination in degrees
m1: float, mass of luminous star in Msun
m2: float, mass of companion in Msun
ecc: float, eccentricity
'''
G = 6.6743e-11 # SI units
Msun = 1.98840987069805e+30
Kstar = (2*np.pi*G*(m2*Msun) * (m2/(m1 + m2))**2 * np.sin(inc_deg*np.pi/180)**3 / \
(period*86400 * (1 - ecc**2)**(3/2)))**(1/3)
return Kstar/1000 # km/s
def lnL_with_covariance(y_vec, mu_vec, Sigma_mat):
'''
calculates the likelihood of a vector y_vec, given a multivariate Gaussian likelihood characterized by
a mean mu_vec and a covariance matrix Sigma_ma
'''
return -0.5*np.dot((y_vec - mu_vec).T, np.dot(np.linalg.inv(Sigma_mat), y_vec - mu_vec)) \
-0.5*np.log(np.linalg.det(Sigma_mat)) -len(mu_vec)/2*np.log(2*np.pi)
def ln_likelihood(theta, fit_rvs = False):
'''
theta = ra, dec, parallax, pmra, pmdec, period, ecc, inc_deg, omega_deg, w_deg, t_peri, v_com, m2, m1
fit_rvs: bool; whether to include the RVs or not
'''
ra, dec, parallax, pmra, pmdec, period, ecc, inc_deg, omega_deg, w_deg, t_peri, v_com, m2, m1 = theta
inc, omega, w = inc_deg*np.pi/180, omega_deg*np.pi/180, w_deg*np.pi/180
if fit_rvs:
Kstar_kms = get_Kstar_kms(period = period, inc_deg = inc_deg, m1 = m1, m2 = m2, ecc = ecc)
rv_pred = get_radial_velocities(t = jds_BH1-jd_ref, P = period, T_p = t_peri, ecc = ecc, K = Kstar_kms,
omega = w, gamma = v_com)
lnL = -0.5*np.sum((rv_pred - rvs_BH1)**2/rv_errs_BH1**2)
else:
lnL = 0
a0_mas = get_a0_mas(period = period, m1 = m1, m2 = m2, parallax = parallax)
A_pred = a0_mas*( np.cos(w)*np.cos(omega) - np.sin(w)*np.sin(omega)*np.cos(inc) )
B_pred = a0_mas*( np.cos(w)*np.sin(omega) + np.sin(w)*np.cos(omega)*np.cos(inc) )
F_pred = -a0_mas*( np.sin(w)*np.cos(omega) + np.cos(w)*np.sin(omega)*np.cos(inc) )
G_pred = -a0_mas*( np.sin(w)*np.sin(omega) - np.cos(w)*np.cos(omega)*np.cos(inc) )
y_vec = np.array([ra, dec, parallax, pmra, pmdec, A_pred, B_pred, F_pred, G_pred, ecc, period, t_peri])
lnL += lnL_with_covariance(y_vec = y_vec, mu_vec = mu_vec_BH1, Sigma_mat = cov_mat_BH1)
lnL += -0.5*(m1 - m1_obs)**2/m1_err**2 # prior on the mass of the luminous star
# put a flat prior on gamma to avoid walkers running to infinity
if np.abs(v_com - 50) > 50:
lnL += -np.inf
if np.isfinite(lnL):
return lnL
else:
return -np.inf
In [6]:
# first the "no-RVs" fit
ndim = 14
nwalkers = 64
p0 = np.array([262.17120816, -0.58109202, 2.09545159, -7.70205044,
-25.85042075, 185.76565789, 0.48893589, 126, 97.8, 12.8, -1.1, 46.6, 9.62, 0.93])
p0 = np.tile(p0, (nwalkers, 1))
p0 += p0 * 1.0e-10 * np.random.normal(size=(nwalkers, ndim))
sampler_norvs = emcee.EnsembleSampler(nwalkers, ndim, ln_likelihood, args=[False])
# 5000 step burn-in
state = sampler_norvs.run_mcmc(p0, 5000)
# 3000 more steps
sampler_norvs.reset()
state = sampler_norvs.run_mcmc(state, 3000)
/Users/kareem/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:16: RuntimeWarning: invalid value encountered in double_scalars app.launch_new_instance()
In [7]:
# next, the "with-rvs" fit. This will take a few minutes because predicting RVs is moderately expensive
# (in production runs we use a compiled C function to speed up the RV prediction)
sampler_withrvs = emcee.EnsembleSampler(nwalkers, ndim, ln_likelihood, args=[True])
# 5000 step burn-in
state = sampler_withrvs.run_mcmc(p0, 5000)
# 5000 more steps
sampler_withrvs.reset()
state = sampler_withrvs.run_mcmc(state, 3000)
In [8]:
labels = [r'$P_{\rm orb}\,\,[\rm days]$', r'$\rm ecc$', r'$\rm inc\,\,[\rm deg]$', r'$\Omega\,\,[\rm deg]$',
r'$\omega\,\,[\rm deg]$', r'$T_{p}\,\,[\rm days]$', r'$\gamma\,\,[\rm km\,s^{-1}]$',
r'$M_2\,\,[M_{\odot}]$', r'$M_{\star}\,\,[M_{\odot}]$']
# now compare the fits
half_sampler = sampler_withrvs.flatchain.T[5:].T
half_sampler_norvs = sampler_norvs.flatchain.T[5:].T
fig = corner.corner(half_sampler_norvs, labels=labels, show_titles=True, plot_datapoints = False,
plot_density = False, color = 'c')
fig = corner.corner(half_sampler, labels = labels, show_titles=True, plot_datapoints = False,
plot_density = False, color = 'k', fig = fig, label_kwargs = {'fontsize': 20})
fig.axes[4].plot([], [], 'c', label = r'$\rm astrometry\,\,only$')
fig.axes[4].plot([], [], 'k', label = r'$\rm astrometry+RVs$')
fig.axes[4].legend(loc = 'upper right', frameon = False, fontsize=30)
Out[8]:
<matplotlib.legend.Legend at 0x7f8eb91b2c10>
In [9]:
# posterior predictive check for the "with RVs" solution.
from astropy.time import Time
ra, dec, parallax, pmra, pmdec, period, ecc, inc_deg, omega_deg, w_deg, t_peri, \
v_com, m2, m1 = sampler_withrvs.flatchain.T
t_grid = np.linspace(np.min(jds_BH1) - 150, np.max(jds_BH1) + 150, 3000)
randints = np.random.randint(0, len(period), 50)
colors = ['#e41a1c', '#feb24c', '#e7298a', '#66a61e', '#a6761d', '#377eb8', '#756bb1']
f, ax = plt.subplots(2, 1, figsize = (14, 10))
plt.subplots_adjust(hspace=0.2)
for i in randints:
Kstar_kms = get_Kstar_kms(period = period[i], inc_deg = inc_deg[i], m1 = m1[i], m2 = m2[i], ecc = ecc[i])
rv_pred = get_radial_velocities(t = t_grid-jd_ref, P = period[i], T_p = t_peri[i], ecc = ecc[i], K = Kstar_kms,
omega = w_deg[i]*np.pi/180, gamma = v_com[i])
ax[0].plot(Time(t_grid, format= 'jd').decimalyear, rv_pred, 'k', alpha=0.1, rasterized=True)
ax[1].plot(Time(t_grid, format= 'jd').decimalyear, rv_pred, 'k', alpha=0.1, rasterized=True)
for i, name in enumerate(['LAMOST', 'MagE', 'GMOS', 'XSHOOTER', 'FEROS', 'HIRES', 'ESI']):
m = rv_names_BH1 == name
ax[0].errorbar(Time(jds_BH1[m], format='jd').decimalyear, rvs_BH1[m], yerr= rv_errs_BH1[m], fmt='o', color = colors[i], label = name)
ax[1].errorbar(Time(jds_BH1[m], format='jd').decimalyear, rvs_BH1[m], yerr= rv_errs_BH1[m], fmt='o', color = colors[i], label = name)
ax[0].set_xlabel(r'$\rm year$', fontsize=20)
ax[1].set_xlabel(r'$\rm year$', fontsize=20)
ax[0].set_ylabel(r'$\rm RV\,\,[\rm km\,s^{-1}]$', fontsize=20)
ax[1].set_ylabel(r'$\rm RV\,\,[\rm km\,s^{-1}]$', fontsize=20)
ax[1].legend(loc = 'upper left')
ax[0].legend(loc = 'upper left', ncol=4)
ax[1].set_xlim(2022.3, 2022.9)
ax[0].set_xlim(2017, 2023)
ax[1].set_xticks([2022.3, 2022.4, 2022.5, 2022.6, 2022.7, 2022.8, 2022.9])
ax[1].set_xticklabels([2022.3, 2022.4, 2022.5, 2022.6, 2022.7, 2022.8, 2022.9])
ax[0].set_ylim(-10, 165)
ax[1].set_ylim(-10, 165)
import matplotlib.patches as patches
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
def custom_mark_inset(axA, axB, fc='None', ec='k'):
xx = axB.get_xlim()
yy = axB.get_ylim()
xy = (xx[0], yy[0])
width = xx[1] - xx[0]
height = yy[1] - yy[0]
pp = axA.add_patch(patches.Rectangle(xy, width, height, fc=fc, ec=ec))
p1 = axA.add_patch(patches.ConnectionPatch(
xyA=(xx[0], yy[0]), xyB=(xx[0], yy[1]),
coordsA='data', coordsB='data',
axesA=axA, axesB=axB, linewidth=1, color='0.2'))
p2 = axA.add_patch(patches.ConnectionPatch(
xyA=(xx[1], yy[0]), xyB=(xx[1], yy[1]),
coordsA='data', coordsB='data',
axesA=axA, axesB=axB,linewidth=1, color='0.2'))
return pp, p1, p2
pp, p1, p2 = custom_mark_inset(ax[0], ax[1])
In [10]:
# posterior predictive check for the "no RVs" solution.
from astropy.time import Time
ra, dec, parallax, pmra, pmdec, period, ecc, inc_deg, omega_deg, w_deg, t_peri, \
v_com, m2, m1 = sampler_norvs.flatchain.T
t_grid = np.linspace(np.min(jds_BH1) - 150, np.max(jds_BH1) + 150, 3000)
randints = np.random.randint(0, len(period), 50)
colors = ['#e41a1c', '#feb24c', '#e7298a', '#66a61e', '#a6761d', '#377eb8', '#756bb1']
f, ax = plt.subplots(2, 1, figsize = (14, 10))
plt.subplots_adjust(hspace=0.2)
for i in randints:
Kstar_kms = get_Kstar_kms(period = period[i], inc_deg = inc_deg[i], m1 = m1[i], m2 = m2[i], ecc = ecc[i])
rv_pred = get_radial_velocities(t = t_grid-jd_ref, P = period[i], T_p = t_peri[i], ecc = ecc[i], K = Kstar_kms,
omega = w_deg[i]*np.pi/180, gamma = 46.6)
ax[0].plot(Time(t_grid, format= 'jd').decimalyear, rv_pred, 'c', alpha=0.1, rasterized=True)
ax[1].plot(Time(t_grid, format= 'jd').decimalyear, rv_pred, 'c', alpha=0.1, rasterized=True)
for i, name in enumerate(['LAMOST', 'MagE', 'GMOS', 'XSHOOTER', 'FEROS', 'HIRES', 'ESI']):
m = rv_names_BH1 == name
ax[0].errorbar(Time(jds_BH1[m], format='jd').decimalyear, rvs_BH1[m], yerr= rv_errs_BH1[m], fmt='o', color = colors[i], label = name)
ax[1].errorbar(Time(jds_BH1[m], format='jd').decimalyear, rvs_BH1[m], yerr= rv_errs_BH1[m], fmt='o', color = colors[i], label = name)
ax[0].set_xlabel(r'$\rm year$', fontsize=20)
ax[1].set_xlabel(r'$\rm year$', fontsize=20)
ax[0].set_ylabel(r'$\rm RV\,\,[\rm km\,s^{-1}]$', fontsize=20)
ax[1].set_ylabel(r'$\rm RV\,\,[\rm km\,s^{-1}]$', fontsize=20)
ax[1].legend(loc = 'upper left')
ax[0].legend(loc = 'upper left', ncol=4)
ax[1].set_xlim(2022.3, 2022.9)
ax[0].set_xlim(2017, 2023)
ax[1].set_xticks([2022.3, 2022.4, 2022.5, 2022.6, 2022.7, 2022.8, 2022.9])
ax[1].set_xticklabels([2022.3, 2022.4, 2022.5, 2022.6, 2022.7, 2022.8, 2022.9])
ax[0].set_ylim(-10, 195)
ax[1].set_ylim(-10, 195)
import matplotlib.patches as patches
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
def custom_mark_inset(axA, axB, fc='None', ec='k'):
xx = axB.get_xlim()
yy = axB.get_ylim()
xy = (xx[0], yy[0])
width = xx[1] - xx[0]
height = yy[1] - yy[0]
pp = axA.add_patch(patches.Rectangle(xy, width, height, fc=fc, ec=ec))
p1 = axA.add_patch(patches.ConnectionPatch(
xyA=(xx[0], yy[0]), xyB=(xx[0], yy[1]),
coordsA='data', coordsB='data',
axesA=axA, axesB=axB, linewidth=1, color='0.2'))
p2 = axA.add_patch(patches.ConnectionPatch(
xyA=(xx[1], yy[0]), xyB=(xx[1], yy[1]),
coordsA='data', coordsB='data',
axesA=axA, axesB=axB,linewidth=1, color='0.2'))
return pp, p1, p2
pp, p1, p2 = custom_mark_inset(ax[0], ax[1])
In [ ]: