In [1]:
%matplotlib inline
import numpy as np
import astropy.units as u
from astropy.units import cds
cds.enable
import matplotlib.pyplot as plt
Set relative positions and velocities for all 4 stellar objects. Mass of stellar outcome is estimated using an efficiency of $0.3$ and each condensation mass.
The mass unit is $M_\odot$, the position unit is $au$, and the velocity unit is km s$^{-1}$.
In [2]:
# YSO, Cond1, Cond2, Cond3
eff=0.3
mass=np.array([0.1, eff*0.362, eff*0.26, eff*0.3])
x_au=np.array([0.0,-8251.1, 248.9, 3873.9])
y_au=np.array([0.0, 7934.9,3309.9,-3315.1])
z_au=np.array([0.0, 0.0, 0.0, 0.0])
v_c=np.array([10.21, 10.43, 10.23, 10.30])
The center of mass ($\vec{r}_{com}$ and $\vec{v}_{com}$) is calculated as $$\vec{r}_{com} = \frac{ 1}{\Sigma_i m_i} \Sigma_i \vec{r}_i m_i$$ $$\vec{v}_{com} = \frac{ 1}{\Sigma_i m_i} \Sigma_i \vec{v}_i m_i.$$
In [3]:
x_com= np.sum(x_au*mass)/np.sum(mass)
y_com= np.sum(y_au*mass)/np.sum(mass)
z_com= np.sum(z_au*mass)/np.sum(mass)
v_com= np.sum(v_c*mass)/np.sum(mass)
print 'r_com = ', (np.array([x_com,y_com,z_com])*u.au).to(u.m)
print 'v_com = ', v_com, 'km/s'
Display system configuration
In [4]:
fig=plt.figure( figsize=( 4.0, 4.0) )
ax= plt.subplot()
color_sym=['red','blue','green','orange']
plt.scatter( x_au[0], y_au[0], color=color_sym[0], s=100, marker='*')
plt.scatter( x_au[1:4], y_au[1:4], color=color_sym[1:4], s=100, marker='o', edgecolor='k')
plt.scatter( x_com, y_com, color='black', marker='+', s=100)
plt.setp( ax, xlabel="x (au)")
plt.setp( ax, ylabel="y (au)")
plt.xlim(1e4, -1e4)
plt.ylim(-1e4, 1e4)
ax.locator_params(axis = 'both', nbins = 5)
Determine binding energy for entire system. It calculates the Kinetic and Potential energies to then determine if the system is bound. The total potential energy of the system is calculated as $$E_{pot, tot} = \Sigma_{i\ne j} \frac{-G m_i m_j}{\left|\vec{r}_{i}-\vec{r}_{j}\right|},$$ while the total kinetic energy of the system is calculated as $$E_{kin, tot}= \Sigma_i \frac{m_i \left| \vec{v}_i -\vec{v}_{com}\right|^2 }{2},$$ where $\vec{v}_{com}$ is the velocity in the center of mass frame.
The energies are all calculated and stored in MKS units (Joules).
In [5]:
nstars=len(x_au)
Epot_i=np.zeros(nstars)*u.J
Ekin_i=np.zeros(nstars)*u.J
#Kinetic energy for each object
Ekin_i = ( 0.5*mass*u.Msun* ((v_c-v_com)*u.km/u.s)**2 ).to(u.J)
for i in range(nstars):
for j in range(nstars):
if i != j:
rij=[x_au[i]-x_au[j], y_au[i]-y_au[j], z_au[i]-z_au[j]]
Epot_i[i] += -( mass[i]*mass[j]/( np.linalg.norm(rij))*cds.G*u.Msun*u.Msun/u.au).to(u.J)
print 'Ekin+Epot Ekin Epot'
print np.array([(Ekin_i+Epot_i).T, Ekin_i.T, Epot_i.T]).T * Epot_i.unit
Kinetic-Potential energy ratio for entire system (bound if <1)
In [6]:
print 'KE/PE ratio=',np.round(np.sum(Ekin_i)/np.abs(np.sum(Epot_i)), decimals=3)
Here we define some functions needed to calculate the binding energy: Calculate the binding energy of the system as $$E_\textrm{bind}=M_1 M_2\left(\frac{(V_1 - V_2)^2}{2(M_1+M_2)}-\frac{G}{|r_1-r_2|} \right).$$
The eccentricity is calculated as $$e=\sqrt{1 + 2 \frac{(M_1 + M_2) L_{tot}^2 E_\textrm{bind}}{G^2 M_1^3 M_2^3}},$$ the semi-major axis is calculated as $$sma = \frac{GM_1 M_2}{2 \left| E_\textrm{bind} \right|} $$ while the period of the system is determined by $$per= 2\pi \sqrt{ \frac{ sma^3}{G (M_1 + M_2)}}.$$
In [7]:
def binary_orbit(v1,v2,r1,r2,m1,m2):
""" Function to determine sma, per and e of a binary system
sma: semi-major axis. Units of m
per: period of the system. Units of yr.
e: eccentricity of the system.
Input parameters are
v1 and v2: velocities for objects 1 and 2. Units of m/s
r1 and r2: positions for objects 1 and 2. Units of m.
m1 and m2: masses for objects 1 and 2. Units of kg.
"""
sma=0.0
per=0.0
e=0.0
Etot=0.0
v_ij=np.linalg.norm(v1 -v2)*v1.unit
r_ij=np.linalg.norm(r1 -r2)*r1.unit
print '|v_ij|', v_ij
print '|r_ij|', r_ij
binary_Etot= ( m1*m2*( v_ij**2/(2*(m1+m2)) - cds.G/r_ij ) ).to(u.J)
print 'E_bin=',binary_Etot
L12 =np.cross( (r1-r2), (v1-v2)) * r1.unit * v1.unit
L12_tot=np.linalg.norm(L12)*L12.unit * m1*m2/(m1+m2)
print 'Angular momentum components', L12**2
print 'Total angular momentum', L12_tot
# Determine the eccentricity of the system
e=np.sqrt(1. + (2.*(m1 + m2)* L12_tot**2 *binary_Etot)/( cds.G**2 * m1**3 * m2**3 ))
sma=( (cds.G*m1*m2)/(2.*np.abs(binary_Etot)) ).to(u.m)
per=( 2*np.pi*np.sqrt((sma**3)/(cds.G*(m1 + m2))) ).to(u.yr)
return sma,per,e,binary_Etot
Calculate the Valtonen stability criterion: $$Q_v=3\left(1+\frac{m_{out}}{m_{in}}\right)^{2/3}\left(1-e_{out}\right)^{-1/6}\left( \frac{7}{4} -\frac{1}{2}\cos{i} - \cos^2{i} \right)^{1/3}.$$
Systems with $a_{out}/a_{in}>Q_v$ are stable for at least $10^4$ revolutions of the outer component.
In [8]:
def valt_crit( m_in, m_out, e_out, inclination):
"""
Calculate the Valtoren stability criterion.
m_in and m_out are the outer and inner component mass (in the same units).
e_out is the eccentricity of the outer component.
inclination is the relative inclination of the orbits, in degrees.
"""
if e_out > 1:
print 'e_out > 1, it is not a closed system'
return -1.0
else:
cos_inclination = np.cos( np.radians(inclination) )
return 3. * (1. + m_out/m_in)**(2./3.) * (1. - e_out)**(-1./6.) * (7./4. - cos_inclination/2. - cos_inclination**2)**(1./3.)
Define variables needed in the calculations based on the objects properties
In [9]:
# YSO
i=0
r1=(np.array([ x_au[i], y_au[i], z_au[i]])*u.au).to(u.m)
v1=(np.array([ 0, 0, v_c[i]])*u.km/u.s).to(u.m/u.s)
m1=(mass[i]*u.solMass).to(u.kg)
# B5-Cond2
i=2
r2=(np.array([ x_au[i], y_au[i], z_au[i]])*u.au).to(u.m)
v2=(np.array([ 0, 0, v_c[i]])*u.km/u.s).to(u.m/u.s)
m2=(mass[i]*u.solMass).to(u.kg)
Calculate binary orbit properties
In [10]:
print 'Original'
sma,per,e,Etot = binary_orbit(v1,v2,r1,r2,m1,m2)
In [11]:
print np.round(sma.to(u.au)), per.to(u.yr), np.round(e, decimals=3), Etot
The binary system is bound
In [12]:
print 'masses = 0.2Msun'
sma_1b,per_1b,e_1b,Etot_1b = binary_orbit(v1,v2,r1,r2,(0.2*u.solMass).to(u.kg),(0.2*u.solMass).to(u.kg))
In [13]:
print np.round(sma_1b.to(u.au)), per_1b.to(u.yr), np.round(e_1b, decimals=3), Etot_1b
First calculate the center of mass properties for pair B5-IRS1 and B5-Cond2
In [14]:
x_com_YSO_Cond2= np.sum((x_au*mass)[[0,2]])/np.sum(mass[[0,2]])
y_com_YSO_Cond2= np.sum((y_au*mass)[[0,2]])/np.sum(mass[[0,2]])
z_com_YSO_Cond2= np.sum((z_au*mass)[[0,2]])/np.sum(mass[[0,2]])
v_com_YSO_Cond2= np.sum((v_c *mass)[[0,2]])/np.sum(mass[[0,2]])
mass_YSO_Cond2 = np.sum(mass[[0,2]])
red_mass_YSO_Cond2 = 1./np.sum(1./mass[[0,2]])
In [15]:
# YSO + B5-Cond2
r1=(np.array([ x_com_YSO_Cond2, y_com_YSO_Cond2, z_com_YSO_Cond2])*u.au).to(u.m)
v1=(np.array([ 0, 0, v_com_YSO_Cond2])*u.km/u.s).to(u.m/u.s)
m1=(mass_YSO_Cond2*u.solMass).to(u.kg)
m1_red=(red_mass_YSO_Cond2*u.solMass).to(u.kg)
# B5-Cond1
i=1
r2=(np.array([ x_au[i], y_au[i], z_au[i]])*u.au).to(u.m)
v2=(np.array([ 0, 0, v_c[i]])*u.km/u.s).to(u.m/u.s)
m2=(mass[i]*u.solMass).to(u.kg)
In [16]:
print 'Calculation using total mass for stability criterion'
sma_triple,per_triple,e_triple,Etot_triple = binary_orbit(v1,v2,r1,r2,m1,m2)
valt_triple=valt_crit(m1, m2, e_triple, 0.0)
print sma_triple, per_triple, e_triple, Etot_triple, valt_triple
print ''
print '>>>stable if ', np.round(sma_triple/sma, decimals=2), '>', np.round(valt_triple, decimals=2)
In [17]:
print 'Calculation using the reduced mass for stability criterion'
sma_triple_red,per_triple_red,e_triple_red,Etot_triple_red = binary_orbit(v1,v2,r1,r2,m1_red,m2)
valt_triple_red=valt_crit(m1_red, m2, e_triple_red, 0.0)
print sma_triple_red,per_triple_red,e_triple_red,Etot_triple_red, valt_triple_red
print ''
print '>>>stable if ', np.round(sma_triple_red/sma, decimals=2), '>', np.round(valt_triple_red, decimals=2)
The system is bound only when using the total mass for the orbital analysis, but it is not stable following the Valtoren criterion.
In [18]:
# YSO + B5-Cond2
r1=(np.array([ x_com_YSO_Cond2, y_com_YSO_Cond2, z_com_YSO_Cond2])*u.au).to(u.m)
v1=(np.array([ 0, 0, v_com_YSO_Cond2])*u.km/u.s).to(u.m/u.s)
m1=(mass_YSO_Cond2*u.solMass).to(u.kg)
m1_red=(red_mass_YSO_Cond2*u.solMass).to(u.kg)
# B5-Cond3
i=3
r2=(np.array([ x_au[i], y_au[i], z_au[i]])*u.au).to(u.m)
v2=(np.array([ 0, 0, v_c[i]])*u.km/u.s).to(u.m/u.s)
m2=(mass[i]*u.solMass).to(u.kg)
In [19]:
print 'Calculation using total mass for stability criterion'
sma_triple,per_triple,e_triple,Etot_triple = binary_orbit(v1,v2,r1,r2,m1,m2)
print 'Inclination of 0deg'
valt_triple=valt_crit(m1, m2, e_triple, 0.0)
print sma_triple, per_triple, e_triple, Etot_triple, valt_triple
print 'Inclination of 45deg'
valt_triple_45=valt_crit(m1, m2, e_triple, 45.0)
print sma_triple, per_triple, e_triple, Etot_triple, valt_triple_45
print 'i=0'
print '>>>stable if ', np.round(sma_triple/sma, decimals=2), '>', np.round(valt_triple, decimals=2)
print 'i=45'
print '>>>stable if ', np.round(sma_triple/sma, decimals=2), '>', np.round(valt_triple_45, decimals=2)
The system is bound, but the Valtonen criterion suggest the system is unstable.
In [20]:
print 'Calculation using reduced mass for stability criterion'
sma_triple_red,per_triple_red,e_triple_red,Etot_triple_red = binary_orbit(v1,v2,r1,r2,m1_red,m2)
valt_triple_red=valt_crit(m1_red, m2, e_triple_red, 0.0)
valt_triple_red_45=valt_crit(m1_red, m2, e_triple_red, 45.0)
print sma_triple_red,per_triple_red,e_triple_red,Etot_triple_red, valt_triple_red, valt_triple_red_45
print ''
print 'i=0'
print '>>>stable if ', np.round(sma_triple_red/sma, decimals=2), '>', np.round(valt_triple_red, decimals=2)
print 'i=45'
print '>>>stable if ', np.round(sma_triple_red/sma, decimals=2), '>', np.round(valt_triple_red_45, decimals=2)
The system is bound, but the Valtonen criterion suggest the system is unstable.