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



# Multiple system parameters

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'




r_com =  [ -2.09740608e+14   3.26343957e+14   0.00000000e+00] m
v_com =  10.2990918747 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)






# Is the entire system bound?

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




Ekin+Epot             Ekin              Epot
[[ -8.14852405e+33   7.89410351e+32  -8.93793440e+33]
[ -2.41169871e+33   1.85092855e+33  -4.26262726e+33]
[ -6.96334640e+33   3.70318303e+32  -7.33366470e+33]
[ -5.79933036e+33   7.38177344e+28  -5.79940417e+33]] J



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)




KE/PE ratio= 0.114



# Determination of binary orbits

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:
return 3. * (1. + m_out/m_in)**(2./3.) * (1. - e_out)**(-1./6.) * (7./4. - cos_inclination/2. - cos_inclination**2)**(1./3.)



## Case 1: Closest pair, B5-IRS1 and B5-Cond2

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)




Original
|v_ij| 20.0 m / s
|r_ij| 4.96552026021e+14 m
E_bin= -4.13037758634e+33 J
Angular momentum components [  9.80709904e+31   5.54575410e+29   0.00000000e+00] m4 / s2
Total angular momentum 8.6561738794e+44 kg m2 / s




In [11]:

print np.round(sma.to(u.au)), per.to(u.yr), np.round(e, decimals=3), Etot




1667.0 161248.43628 yr 0.992 -4.13037758634e+33 J



The binary system is bound

### Case 1b: Closest pair, B5-IRS1 and B5-Cond2 but with masses = 0.2M$_{\odot}$



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))




masses = 0.2Msun
|v_ij| 20.0 m / s
|r_ij| 4.96552026021e+14 m
E_bin= -2.12310392725e+34 J
Angular momentum components [  9.80709904e+31   5.54575410e+29   0.00000000e+00] m4 / s2
Total angular momentum 1.97538326992e+45 kg m2 / s




In [13]:

print np.round(sma_1b.to(u.au)), per_1b.to(u.yr), np.round(e_1b, decimals=3), Etot_1b




1663.0 107189.279734 yr 0.996 -2.12310392725e+34 J



## Case 2: Possible triple system, (B5-IRS1 + B5-Cond2) and B5-Cond1

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)




Calculation using total mass for stability criterion
|v_ij| 211.235955056 m / s
|r_ij| 1.58277859667e+15 m
E_bin= -2.31716597294e+32 J
Angular momentum components [  4.19893108e+34   6.97937567e+34   0.00000000e+00] m4 / s2
Total angular momentum 4.48557652603e+46 kg m2 / s
1.10141684661e+16 m 37312055.661 yr 0.856296133335 -2.31716597294e+32 J 3.58724168659

>>>stable if  44.18 > 3.59




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)




Calculation using the reduced mass for stability criterion
|v_ij| 211.235955056 m / s
|r_ij| 1.58277859667e+15 m
E_bin= 5.91642490696e+32 J
Angular momentum components [  4.19893108e+34   6.97937567e+34   0.00000000e+00] m4 / s2
Total angular momentum 2.07637818248e+46 kg m2 / s
e_out > 1, it is not a closed system
1.06195007186e+15 m 1531773.48091 yr 2.49044539722 5.91642490696e+32 J -1.0

>>>stable if  4.26 > -1.0



The system is bound only when using the total mass for the orbital analysis, but it is not stable following the Valtoren criterion.

## Case 3: Possible triple system, (B5-IRS1 + B5-Cond2) and B5-Cond3



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)




Calculation using total mass for stability criterion
|v_ij| 81.2359550562 m / s
|r_ij| 9.08540762482e+14 m
E_bin= -4.26360780597e+33 J
Angular momentum components [  3.35401791e+33   2.09333379e+33   0.00000000e+00] m4 / s2
Total angular momentum 8.77560451362e+45 kg m2 / s
Inclination of 0deg
4.96071518113e+14 m 368815.138941 yr 0.831471328849 -4.26360780597e+33 J 3.34040264856
Inclination of 45deg
4.96071518113e+14 m 368815.138941 yr 0.831471328849 -4.26360780597e+33 J 5.11281752061
i=0
>>>stable if  1.99 > 3.34
i=45
>>>stable if  1.99 > 5.11



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)




Calculation using reduced mass for stability criterion
|v_ij| 81.2359550562 m / s
|r_ij| 9.08540762482e+14 m
E_bin= -9.52775816759e+32 J
Angular momentum components [  3.35401791e+33   2.09333379e+33   0.00000000e+00] m4 / s2
Total angular momentum 4.32657941945e+45 kg m2 / s
5.46494020077e+14 m 603499.536713 yr 0.662489851864 -9.52775816759e+32 J 4.76743617974 7.29703385867

i=0
>>>stable if  2.19 > 4.77
i=45
>>>stable if  2.19 > 7.3



The system is bound, but the Valtonen criterion suggest the system is unstable.