mcfacts.physics.dynamics

Module for handling dynamical interactions.

Contains multiple functions which are each mocked up versions of a dynamical mechanism. Of varying fidelity to reality. Also contains GW orbital evolution for BH in the inner disk, which should probably move elsewhere.

mcfacts.physics.dynamics.bh_near_smbh(smbh_mass, disk_bh_pro_orbs_a, disk_bh_pro_masses, disk_bh_pro_orbs_ecc, timestep_duration_yr, inner_disk_outer_radius, disk_inner_stable_circ_orb, r_g_in_meters)

Evolve semi-major axis of single BH near SMBH according to Peters64

Test whether there are any BH near SMBH. Flag if anything within min_safe_distance (default=50r_g) of SMBH. Time to decay into SMBH can be parameterized from Peters(1964) as: .. math:: t_{gw} =38Myr (1-e^2)(7/2) (a/50r_{g})^4 (M_{smbh}/10^8M_{sun})^3 (m_{bh}/10M_{sun})^{-1}

Parameters:

  • smbh_mass (float) – Mass [M_sun] of supermassive black hole

  • disk_bh_pro_orbs_a (numpy.ndarray) – Orbital semi-major axes [r_{g,SMBH}] of prograde singleton BH at start of a timestep (math:r_g=GM_{SMBH}/c^2) with float type

  • disk_bh_pro_masses (numpy.ndarray) – Masses [M_sun] of prograde singleton BH at start of timestep with float type

  • disk_bh_pro_orbs_ecc (numpy.ndarray) – Orbital eccentricity [unitless] of singleton prograde BH with float type

  • timestep_duration_yr (float) – Length of timestep [yr]

  • inner_disk_outer_radius (float) – Outer radius of the inner disk [r_{g,SMBH}]

  • disk_inner_stable_circ_orb (float) – Innermost stable circular orbit around the SMBH [r_{g,SMBH}]

  • r_g_in_meters (float) – Gravitational radius of the SMBH in meters

Returns:

disk_bh_pro_orbs_a – Semi-major axis [r_{g,SMBH}] of prograde singleton BH at end of timestep assuming only GW evolution

Return type:

numpy.ndarray

mcfacts.physics.dynamics.bin_recapture(bin_mass_1_all, bin_mass_2_all, bin_orb_a_all, bin_orb_inc_all, timestep_duration_yr)

Recapture BBH that has orbital inclination >0 post spheroid encounter

Parameters:

  • blackholes_binary (AGNBinaryBlackHole) – binary black holes

  • timestep_duration_yr (float) – Length of timestep [yr]

Returns:

blackholes_binary – Binary black holes with binary orbital inclination [radian] updated

Return type:

AGNBinaryBlackHole

Notes

Purely bogus scaling does not account for real disk surface density. From Fabj+20, if i<5deg (=(5deg/180deg)*pi=0.09rad), time to recapture a BH in SG disk is 1Myr (M_b/10Msun)^-1(R/10^4r_g) if i=[5,15]deg =(0.09-0.27rad), time to recapture a BH in SG disk is 50Myrs(M_b/10Msun)^-1 (R/10^4r_g) For now, ignore if i>15deg (>0.27rad)

mcfacts.physics.dynamics.bin_spheroid_encounter(smbh_mass, timestep_duration_yr, bin_mass_1_all, bin_mass_2_all, bin_orb_a_all, bin_sep_all, bin_ecc_all, bin_orb_ecc_all, bin_orb_inc_all, time_passed, nsc_bh_imf_powerlaw_index, delta_energy_strong, nsc_spheroid_normalization, harden_energy_delta_mu, harden_energy_delta_sigma, r_g_in_meters)

Perturb orbits due to encounters with spheroid (NSC) objects

Parameters:

  • smbh_mass (float) – Mass [M_sun] of supermassive black hole

  • disk_bins_bhbh (numpy.ndarray) – [21, bindex] mixed array containing properties of binary BBH, see add_to_binary_array function for complete description

  • time_passed (float) –

    Current time set [yr]

    nsc_bh_imf_powerlaw_index : float Powerlaw index of nuclear star cluster BH IMF (e.g. M^-2) [unitless]. User set (default = 2).

  • timestep_duration_yr (float) – Length of timestep [yr]

  • nsc_bh_imf_powerlaw_index (float) – Powerlaw index of nuclear star cluster BH IMF (e.g. M^-2) [unitless]. User set (default = 2).

  • delta_energy_strong (float) – Average energy change [units??] per strong encounter

  • nsc_spheroid_normalization (float) – Normalization factor [unitless] determines the departures from sphericity of the initial distribution of perturbers (1.0=spherical)

  • harden_energy_delta_sigma (float) – Average energy exchanged in a strong 2 + 1 interaction that hardens the binary

  • harden_energy_delta_mu (float) – Variance of the energy exchanged in a strong 2 + 1 interaction that hardens the binary

  • r_g_in_meters (float) – Gravitational radius of the SMBH in meters

Returns:

disk_bins_bhbh – updated version of input after dynamical perturbations

Return type:

[21, bindex] mixed array

Notes

Warning: the powerlaw index for the mass of perturbers is for BH but should be for stars, and the mode mass is hardcoded inside the fn

Use Leigh+18 to figure out the rate at which spheroid encounters happen to binaries embedded in the disk Binaries at small disk radii encounter spheroid objects at high rate, particularly early on in the disk lifetime However, orbits at those small radii get captured quickly by the disk.

From Fig.1 in Leigh+18, Rate of sph. encounter = 20/Myr at t=0, normalized to a_bin=1AU, R_disk=10^3r_g or 0.2/10kyr timestep. Introduce a spheroid normalization factor nsc_spheroid_normalization=1 (default) allowing for non-ideal NSC (previous episodes; disky populations etc). Within 1Myr, for a dense model disk (e.g. Sirko & Goodman), most of those inner stellar orbits have been captured by the disk. So rate of sph. encounter ->0/Myr at t=1Myr since those orbits are gone (R<10^3r_g; assuming approx circular orbits!) for SG disk model For TQM disk model, rate of encounter slightly lower but non-zero.

So, inside R_com<10^3r_g: (would be rt of enc =0.2 if nsc_spheroid_normalization=1) Assume: :math:` ext{Rate of encounter} = 0.2 ( ext{nsc_spheroid_normalization}/1)( ext{timestep}/10kyr)^{-1} (R_{com}/10^3r_g)^{-1} (a_{bin}/1r_gM8)^{-2}` Generate random number from uniform [0,1] distribution and if <0.2 (normalized to above condition) then encounter

Encounter rt starts at = :math:` 0.2 ( ext{nsc_spheroid_normalization}/1)( ext{timestep}/10kyr)^{-1} (R_{com}/10^3r_g)^{-1} (a_{bin}/1r_gM8)^{-2}` at t=0 decreases to = :math:` 0. ( ext{nsc_spheroid_normalization}/1)( ext{timestep}/10kyr)^{-1} (R_{com}/10^3r_g)^{-1} (a_{bin}/1r_gM8)^{-2}` (time_passed/1Myr) at R<10^3r_g. Outside: R_com>10^3r_g Normalize to rate at (R_com/10^4r_g) so that rate is non-zero at R_com=[1e3,1e4]r_g after 1Myr. Decrease rate with time, but ensure it goes to zero at R_com<1.e3r_g.

So, rate of sph. encounter = 2/Myr at t=0, normalized to a_bin=1AU, R_disk=10^4r_g which is equivalently Encounter rate = 0.02 (nsc_spheroid_normalization/1)(timestep/10kyr)^-1 (R_com/10^4r_g)^-1 (a_bin/1r_gM8)^2 Drop this by an order of magnitude over 1Myr. Encounter rate = 0.02 (timestep/10kyr)^-1 (R_com/10^4r_g)^-1 (a_bin/1r_gM8)^2 (time_passed/10kyr)^-1/2 so ->0.002 after a Myr For R_com < 10^3r_g:

if time_passed <=1Myr

Encounter rt = 0.02*(nsc_spheroid_normalization/0.1)*(1-(1Myr/time_passed))(timestep/10kyr)^{-1}(R_com/10^3r_g)^-1 (a_bin/1r_gM8)^2 ….(1)

if time_passed >1Myr

Encounter rt = 0

For R_com > 10^3r_g:

Encounter rt = 0.002 (nsc_spheroid_normalization/0.1) (timestep/10kyr)^-1 (R_com/10^4r_g)^-1 (a_bin/1r_gM8)^2 (time_passed/10kyr)^-1/2 ….(2)

Return corrected binary with spin angles projected onto new L_bin. So can calculate chi_p (in plane components of spin) Return new binary inclination angle w.r.t disk Harden/soften/ionize binary as appropriate

Orbital angular momentum: Binary orbital angular momentum is

L_bin =M_bin*v_orb_bin X R_com

Spheroid orbital angular momentum is

L3=m3*v3 X R3

where m3,v3,R3 are the mass, velocity and semi-major axis of tertiary encounter.

Draw m3 from IMF random distrib. BUT mostly stars early on! TO DO: Switch from spheroid stars to spheroid BH at late time Draw a3 from uniform distribution a3=[10^-0.5,0.5]a_bbh say. v_3= c/sqrt(R_3) Ratio of L3/Lbin =(m3/M_bin)*sqrt(R3/R_com)….(3) so L3 = ratio*Lbin

Resultant L_bin must be the resultant in a parallelogram of L3 (one side) and L_bin(other side)

Angle of encounter: If angle of encounter between BBH and M3 (angle_enc)<|90deg|, ie angle_enc is in [0-90deg,270-360deg] then: L_bin_new = sqrt(L3^2 + L_bin^2 + 2L3L_bin cos(angle_enc)) ….(4)

= sqrt( (1+ratio^2)L_bin^2 + 2ratioL_bin^2 cos(angle_enc)) = sqrt((1+ratio^2) + 2ratio*cos(angle_enc)) L_bin_old

else if angle_enc is in [90deg,270 deg] L_bin_new = sqrt(L3^2 + L_bin^2 - 2L3L_bin cos(angle_enc)) ….(5)

= sqrt((1+ratio^2) - 2ratio*cos(angle_end))L_bin_old

and L_bin_new/L_bin = v_b_new x R_com_new/ v_b_old x R_com_old and for Keplerian vels v_b_com = sqrt(GM_smbh/a_com) so

L_bin_new/L_bin_old = sqrt(a_com_new/a_com_old) ….(6)

So new BBH semi-major axis: a_com_new = (L_bin_new/L_bin_old)^2 *(a_com_old) ….(7)

Angle of encounter: M3 has some random angle (i3) in the spheroid wrt disk (i=0deg) & BBH (also presumed i=0deg). But, over time, spheroid population (of STARS) with small inclination angles wrt disk (i=0 deg) are captured by disk (takes ~1Myr in SG disk; Fabj+20) So, at t=0, start with drawing from uniform distribution of i3=[0,360] After 1Myr in a SG disk, we want all the spheroid (star!) encounters inside R=1000r_g to go to zero. Over time remove e.g. i3=[0,+/-15], so draw from [15,345] next timestep Then remove i3 =+/-[15,30] so draw from [30,330] etc. So, if crit_time =1.e6 #1Myr then excluded_angles =(time_passed/crit_time)*180 select from i3 = [excluded angles,360-excluded angles] So:

crit_time=1.e6 if time_passed < crit_time

excluded_angles = (time_passed/crit_time)*180 if R<10^3r_g

#Draw random integer in range [excluded_angles,360-(excluded_angles)] i3 = rng.randint(excluded_angles, 360-(excluded_angles))….(8)

Calculate velocity of encounter compared to a_bin. Ignore what happens to m3, since it’s a random draw from the NSC and we are not tracking individual NSC components. If binary is hard ie GM1M2/a_bin > m3v_rel^2 then:

harden binary to a_bin = a_bin -da_bin and new binary eccentricity e_bin = e_bin + de around com and new binary orb eccentricity e_orb_com = e_orb_com + de

If binary is soft ie GM_bin/a_bin <m3v_rel^2 then:

soften binary to a_bin = a_bin + da_bin and new binary eccentricity e_bin = e_bin + de

mcfacts.physics.dynamics.circular_binaries_encounters_circ_prograde(smbh_mass, disk_bh_pro_orbs_a, disk_bh_pro_masses, disk_bh_pro_orbs_ecc, bin_mass_1, bin_mass_2, bin_orb_a, bin_sep, bin_ecc, bin_orb_ecc, timestep_duration_yr, disk_bh_pro_orb_ecc_crit, delta_energy_strong, disk_radius_outer, harden_energy_delta_mu, harden_energy_delta_sigma, r_g_in_meters)

“Adjust orb ecc due to encounters btw BBH and circularized singles

Parameters:

  • smbh_mass (float) – Mass [M_sun] of supermassive black hole

  • disk_bh_pro_orbs_a (numpy.ndarray) – Orbital semi-major axes [r_{g,SMBH}] of prograde singleton BH at start of a timestep (math:r_g=GM_{SMBH}/c^2) with float type

  • disk_bh_pro_masses (numpy.ndarray) – Masses [M_sun] of prograde singleton BH at start of timestep with float type

  • disk_bh_pro_orbs_ecc (numpy.ndarray) – Orbital eccentricity [unitless] of singleton prograde BH with float type

  • timestep_duration_yr (float) – Length of timestep [yr]

  • disk_bh_pro_orb_ecc_crit (float) – Critical orbital eccentricity [unitless] below which orbit is close enough to circularize

  • delta_energy_strong (float) – Average energy change [units??] per strong encounter

  • disk_bins_bhbh (numpy.ndarray) – [21, bindex] mixed array containing properties of binary BBH, see add_to_binary_array function for complete description

  • disk_radius_outer (float) – Outer radius of the inner disk (Rg)

  • harden_energy_delta_sigma (float) – Average energy exchanged in a strong 2 + 1 interaction that hardens the binary

  • harden_energy_delta_mu (float) – Variance of the energy exchanged in a strong 2 + 1 interaction that hardens the binary

  • r_g_in_meters (float) – Gravitational radius of the SMBH in meters

Returns:

disk_bins_bhbh – [21, bindex] mixed array, updated version of input after dynamical perturbations

Return type:

numpy.ndarray

Notes

Return array of modified binary BH separations and eccentricities perturbed by encounters within f*R_Hill, for circularized singleton population, where f is some fraction/multiple of Hill sphere radius R_H Right now assume f=1. Logic:

  1. Find number of binaries in this timestep given by bindex

  2. Find the binary center of mass (c.o.m.) and corresponding orbital velocities & binary total masses. disk_bins_bhbh[9,:] = bin c.o.m. = [R_bin1_orb_a,R_bin2_orb_a,…]. These are the orbital radii of the bins. disk_bins_bhbh[8,;] = bin_separation =[a_bin1,a_bin2,…] disk_bins_bhbh[2,:]+disk_bins_bhbh[3,:] = mass of binaries disk_bins_bhbh[13,:] = ecc of binary around com disk_bins_bhbh[18,:] = orb. ecc of binary com around SMBH Keplerian orbital velocity of the bin c.o.m. around SMBH: v_bin,i= sqrt(GM_SMBH/R_bin,i_com)= c/sqrt(R_bin,i_com)

  3. Calculate the binary orbital time and N_orbits/timestep For example, since T_orb =2pi sqrt(R_bin_com^3/GM_smbh) and R_bin_com^3/GM_smbh = (10^3r_g)^3/GM_smbh = 10^9 (R_bin_com/10^3r_g)^3 (GM_smbh/c^2)^3/GM_smbh

    = 10^9 (R_bin_com/10^3r_g)^3 (G M_smbh/c^3)^2

    So, .. math:

    T_{orb}
= 2\pi 10^{4.5} (R_{bin,orb a}/10^3r_g)^{3/2} GM_{smbh}/c^3
= 2\pi 10^{4.5} (R_{bin,orb a}/10^3r_g)^{3/2} (6.7e-11*2e38/(3e8)^3)
= 2\pi 10^{4.5} (R_{bin,orb a}/10^3r_g)^{3/2} (13.6e27/27e24)
= \pi 10^{7.5}  (R_{bin,orb a}/10^3r_g)^{3/2}
~ 3.15 yr (R_{bin,orb a}/10^3r_g)^3/2 (M_smbh/10^8Msun)

    i.e. Orbit~3.15yr at 10^3r_g around a 10^8M_{sun} SMBH. Therefore in a timestep=1.e4yr, a binary at 10^3r_g orbits the SMBH N_orbit/timestep =3,000 times.

  4. Calculate binding energy of bins = [GM1M2/sep_bin1, GMiMi+1,sep_bin2, ….] where sep_bin1 is in meters and M1,M2 are binary mass components in kg.

  5. Find those single BH with e>e_crit and their associated semi-major axes a_ecc =[a_ecc1, a_ecc2, ..] and masses m_ecc =[m_ecc1,m_ecc2, ..] and calculate their average velocities v_ecc = [GM_smbh/a_ecc1, GM_smbh/a_ecc2,…]

  6. Where (1-ecc_i)*a_ecc_i < R_bin_j_com < (1+ecc_i)*a_ecc_i, interaction possible

  7. Among candidate encounters, calculate relative velocity of encounter.

    \(v_{peri,i}=\sqrt(Gm_{ecc,i}/a_{ecc,i}[1+ecc,i/1-ecc,i])\) \(v_{apo,i} =\sqrt(Gm_{ecc,i}/a_{ecc,i}[1-ecc,i/1+ecc,i])\) \(v_{ecc,i} =\sqrt(GM/a_{ecc_i})\) ..average Keplerian vel.

    \(v_{rel} = abs(v_{bin,i} - v_{ecc,i})\)

  8. Calculate relative K.E. of tertiary, (1/2)m_ecc_i*v_rel_^2

  9. Compare binding en of binary to K.E. of tertiary.
    Critical velocity for ionization of binary is v_crit, given by:

    :math:`v_{crit} = sqrt(GM_1M_2(M_1+M_2+M_3)/M_3(M_1+M_2)a_{bin})

    If binary is hard ie GM_1M_2/a_bin > m3v_rel^2 then:
    harden binary

    a_bin -> a_bin -da_bin and

    new binary eccentricity

    e_bin -> e_bin + de

    and give +da_bin worth of binding energy (GM_bin/(a_bin -da_bin) - GM_bin/a_bin) to extra eccentricity ecc_i and a_ecc,i of m_ecc,i. Say average en of encounter is de=0.1 (10%) then binary a_bin shrinks by 10%, ecc_bin is pumped by 10% And a_ecc_i shrinks by 10% and ecc_i also shrinks by 10%

    If binary is soft ie GM_bin/a_bin <m3v_rel^2 then:
    if v_rel (effectively v_infty) > v_crit
    ionize binary

    update singleton array with 2 new BH with orbital eccentricity e_crit+de remove binary from binary array

    else if v_rel < v_crit
    soften binary

    a_bin -> a_bin + da_bin and

    new binary eccentricity

    e_bin -> e_bin + de

    and remove -da_bin worth of binary energy from eccentricity of m3.

Note1: Will need to test binary eccentricity each timestep.

If bin_ecc> some value (0.9), check for da_bin due to GW bremsstrahlung at pericenter.

  1. As 4, except now include interactions between binaries and circularized BH. This should give us primarily

    hardening encounters as in Leigh+2018, since the v_rel is likely to be small for more binaries.

Given array of binaries at locations [a_bbh1,a_bbh2] with binary semi-major axes [a_bin1,a_bin2,…] and binary eccentricities [e_bin1,e_bin2,…], find all the single BH at locations a_i that within timestep

either pass between a_i(1-e_i)< a_bbh1 <a_i(1+e_i)

Calculate velocity of encounter compared to a_bin. If binary is hard ie GM1M2/a_bin > m3v_rel^2 then:

harden binary to a_bin = a_bin -da_bin and new binary eccentricity e_bin = e_bin + de around com and new binary orb eccentricity e_orb_com = e_orb_com + de and now give da_bin worth of binding energy to extra eccentricity of m3.

If binary is soft ie GM_bin/a_bin <m3v_rel^2 then:

soften binary to a_bin = a_bin + da_bin and new binary eccentricity e_bin = e_bin + de and take da_bin worth of binary energy from eccentricity of m3.

If binary is unbound ie GM_bin/a_bin << m3v_rel^2 then:

remove binary from binary array add binary components m1,m2 back to singleton arrays with new orbital eccentricities e_1,e_2 from energy of encounter. Equipartition energy so m1v1^2 =m2 v_2^2 and generate new individual orbital eccentricities e1=v1/v_kep_circ and e_2=v_2/v_kep_circ Take energy put into destroying binary from orb. eccentricity of m3.

mcfacts.physics.dynamics.circular_binaries_encounters_circ_prograde_star(smbh_mass, disk_star_pro_orbs_a, disk_star_pro_masses, disk_star_pro_orbs_ecc, disk_star_pro_id_nums, bin_mass_1, bin_mass_2, bin_orb_a, bin_sep, bin_ecc, bin_orb_ecc, bin_id_nums, rstar_rhill_exponent, timestep_duration_yr, disk_bh_pro_orb_ecc_crit, delta_energy_strong, disk_radius_outer, harden_energy_delta_mu, harden_energy_delta_sigma, r_g_in_meters)

“Adjust orb ecc due to encounters btw BBH and circularized singles

Parameters:

  • smbh_mass (float) – Mass [M_sun] of supermassive black hole

  • disk_bh_pro_orbs_a (numpy.ndarray) – Orbital semi-major axes [r_{g,SMBH}] of prograde singleton BH at start of a timestep (math:r_g=GM_{SMBH}/c^2) with float type

  • disk_bh_pro_masses (numpy.ndarray) – Masses [M_sun] of prograde singleton BH at start of timestep with float type

  • disk_bh_pro_orbs_ecc (numpy.ndarray) – Orbital eccentricity [unitless] of singleton prograde BH with float type

  • timestep_duration_yr (float) – Length of timestep [yr]

  • disk_bh_pro_orb_ecc_crit (float) – Critical orbital eccentricity [unitless] below which orbit is close enough to circularize

  • delta_energy_strong (float) – Average energy change [units??] per strong encounter

  • disk_bins_bhbh (numpy.ndarray) – [21, bindex] mixed array containing properties of binary BBH, see add_to_binary_array function for complete description

  • disk_radius_outer (float) – Outer radius of the inner disk (Rg)

  • harden_energy_delta_sigma (float) – Average energy exchanged in a strong 2 + 1 interaction that hardens the binary

  • harden_energy_delta_mu (float) – Variance of the energy exchanged in a strong 2 + 1 interaction that hardens the binary

  • r_g_in_meters (float) – Gravitational radius of the SMBH in meters

Returns:

disk_bins_bhbh – [21, bindex] mixed array, updated version of input after dynamical perturbations

Return type:

numpy.ndarray

Notes

Return array of modified binary BH separations and eccentricities perturbed by encounters within f*R_Hill, for circularized singleton population, where f is some fraction/multiple of Hill sphere radius R_H Right now assume f=1. Logic:

  1. Find number of binaries in this timestep given by bindex

  2. Find the binary center of mass (c.o.m.) and corresponding orbital velocities & binary total masses. disk_bins_bhbh[9,:] = bin c.o.m. = [R_bin1_orb_a,R_bin2_orb_a,…]. These are the orbital radii of the bins. disk_bins_bhbh[8,;] = bin_separation =[a_bin1,a_bin2,…] disk_bins_bhbh[2,:]+disk_bins_bhbh[3,:] = mass of binaries disk_bins_bhbh[13,:] = ecc of binary around com disk_bins_bhbh[18,:] = orb. ecc of binary com around SMBH Keplerian orbital velocity of the bin c.o.m. around SMBH: v_bin,i= sqrt(GM_SMBH/R_bin,i_com)= c/sqrt(R_bin,i_com)

  3. Calculate the binary orbital time and N_orbits/timestep For example, since T_orb =2pi sqrt(R_bin_com^3/GM_smbh) and R_bin_com^3/GM_smbh = (10^3r_g)^3/GM_smbh = 10^9 (R_bin_com/10^3r_g)^3 (GM_smbh/c^2)^3/GM_smbh

    = 10^9 (R_bin_com/10^3r_g)^3 (G M_smbh/c^3)^2

    So, .. math:

    T_{orb}
= 2\pi 10^{4.5} (R_{bin,orb a}/10^3r_g)^{3/2} GM_{smbh}/c^3
= 2\pi 10^{4.5} (R_{bin,orb a}/10^3r_g)^{3/2} (6.7e-11*2e38/(3e8)^3)
= 2\pi 10^{4.5} (R_{bin,orb a}/10^3r_g)^{3/2} (13.6e27/27e24)
= \pi 10^{7.5}  (R_{bin,orb a}/10^3r_g)^{3/2}
~ 3.15 yr (R_{bin,orb a}/10^3r_g)^3/2 (M_smbh/10^8Msun)

    i.e. Orbit~3.15yr at 10^3r_g around a 10^8M_{sun} SMBH. Therefore in a timestep=1.e4yr, a binary at 10^3r_g orbits the SMBH N_orbit/timestep =3,000 times.

  4. Calculate binding energy of bins = [GM1M2/sep_bin1, GMiMi+1,sep_bin2, ….] where sep_bin1 is in meters and M1,M2 are binary mass components in kg.

  5. Find those single BH with e>e_crit and their associated semi-major axes a_ecc =[a_ecc1, a_ecc2, ..] and masses m_ecc =[m_ecc1,m_ecc2, ..] and calculate their average velocities v_ecc = [GM_smbh/a_ecc1, GM_smbh/a_ecc2,…]

  6. Where (1-ecc_i)*a_ecc_i < R_bin_j_com < (1+ecc_i)*a_ecc_i, interaction possible

  7. Among candidate encounters, calculate relative velocity of encounter.

    \(v_{peri,i}=\sqrt(Gm_{ecc,i}/a_{ecc,i}[1+ecc,i/1-ecc,i])\) \(v_{apo,i} =\sqrt(Gm_{ecc,i}/a_{ecc,i}[1-ecc,i/1+ecc,i])\) \(v_{ecc,i} =\sqrt(GM/a_{ecc_i})\) ..average Keplerian vel.

    \(v_{rel} = abs(v_{bin,i} - v_{ecc,i})\)

  8. Calculate relative K.E. of tertiary, (1/2)m_ecc_i*v_rel_^2

  9. Compare binding en of binary to K.E. of tertiary.
    Critical velocity for ionization of binary is v_crit, given by:

    :math:`v_{crit} = sqrt(GM_1M_2(M_1+M_2+M_3)/M_3(M_1+M_2)a_{bin})

    If binary is hard ie GM_1M_2/a_bin > m3v_rel^2 then:
    harden binary

    a_bin -> a_bin -da_bin and

    new binary eccentricity

    e_bin -> e_bin + de

    and give +da_bin worth of binding energy (GM_bin/(a_bin -da_bin) - GM_bin/a_bin) to extra eccentricity ecc_i and a_ecc,i of m_ecc,i. Say average en of encounter is de=0.1 (10%) then binary a_bin shrinks by 10%, ecc_bin is pumped by 10% And a_ecc_i shrinks by 10% and ecc_i also shrinks by 10%

    If binary is soft ie GM_bin/a_bin <m3v_rel^2 then:
    if v_rel (effectively v_infty) > v_crit
    ionize binary

    update singleton array with 2 new BH with orbital eccentricity e_crit+de remove binary from binary array

    else if v_rel < v_crit
    soften binary

    a_bin -> a_bin + da_bin and

    new binary eccentricity

    e_bin -> e_bin + de

    and remove -da_bin worth of binary energy from eccentricity of m3.

Note1: Will need to test binary eccentricity each timestep.

If bin_ecc> some value (0.9), check for da_bin due to GW bremsstrahlung at pericenter.

  1. As 4, except now include interactions between binaries and circularized BH. This should give us primarily

    hardening encounters as in Leigh+2018, since the v_rel is likely to be small for more binaries.

Given array of binaries at locations [a_bbh1,a_bbh2] with binary semi-major axes [a_bin1,a_bin2,…] and binary eccentricities [e_bin1,e_bin2,…], find all the single BH at locations a_i that within timestep

either pass between a_i(1-e_i)< a_bbh1 <a_i(1+e_i)

Calculate velocity of encounter compared to a_bin. If binary is hard ie GM1M2/a_bin > m3v_rel^2 then:

harden binary to a_bin = a_bin -da_bin and new binary eccentricity e_bin = e_bin + de around com and new binary orb eccentricity e_orb_com = e_orb_com + de and now give da_bin worth of binding energy to extra eccentricity of m3.

If binary is soft ie GM_bin/a_bin <m3v_rel^2 then:

soften binary to a_bin = a_bin + da_bin and new binary eccentricity e_bin = e_bin + de and take da_bin worth of binary energy from eccentricity of m3.

If binary is unbound ie GM_bin/a_bin << m3v_rel^2 then:

remove binary from binary array add binary components m1,m2 back to singleton arrays with new orbital eccentricities e_1,e_2 from energy of encounter. Equipartition energy so m1v1^2 =m2 v_2^2 and generate new individual orbital eccentricities e1=v1/v_kep_circ and e_2=v_2/v_kep_circ Take energy put into destroying binary from orb. eccentricity of m3.

mcfacts.physics.dynamics.circular_binaries_encounters_ecc_prograde(smbh_mass, disk_bh_pro_orbs_a, disk_bh_pro_masses, disk_bh_pro_orbs_ecc, bin_mass_1, bin_mass_2, bin_orb_a, bin_sep, bin_ecc, bin_orb_ecc, timestep_duration_yr, disk_bh_pro_orb_ecc_crit, delta_energy_strong, disk_radius_outer, r_g_in_meters)

“Adjust orb eccentricities due to encounters between BBH and eccentric single BHs

Return array of modified binary BH separations and eccentricities perturbed by encounters within f*R_Hill, for eccentric singleton population, where f is some fraction/multiple of Hill sphere radius R_H Right now assume f=1.

Parameters:

  • smbh_mass (float) – Mass [M_sun] of supermassive black hole

  • disk_bh_pro_orbs_a (numpy.ndarray) – Orbital semi-major axes [r_{g,SMBH}] of prograde singleton BH at start of a timestep (math:r_g=GM_{SMBH}/c^2) with float type

  • disk_bh_pro_masses (numpy.ndarray) – Masses [M_sun] of prograde singleton BH at start of timestep with float type

  • disk_bh_pro_orbs_ecc (numpy.ndarray) – Orbital eccentricity [unitless] of singleton prograde BH with float type

  • timestep_duration_yr (float) – Length of timestep [yr]

  • disk_bh_pro_orb_ecc_crit (float) – Critical orbital eccentricity [unitless] below which orbit is close enough to circularize

  • delta_energy_strong (float) – Average energy change [units??] per strong encounter

  • disk_bins_bhbh (numpy.ndarray) – [21, bindex] mixed array containing properties of binary BBH, see add_to_binary_array function for complete description

  • disk_radius_outer (float) – Outer radius of the inner disk (Rg)

  • r_g_in_meters (float) – Gravitational radius of the SMBH in meters

Returns:

disk_bins_bhbh – [21, bindex] mixed array, updated version of input after dynamical perturbations

Return type:

numpy.ndarray

Notes

Logic:

  1. Find number of binaries in this timestep given by bindex

  2. Find the binary center of mass (c.o.m.) and corresponding orbital velocities & binary total masses. disk_bins_bhbh[9,:] = bin c.o.m. = [R_bin1_orb_a,R_bin2_orb_a,…]. These are the orbital radii of the bins. disk_bins_bhbh[8,;] = bin_separation =[a_bin1,a_bin2,…] disk_bins_bhbh[2,:]+disk_bins_bhbh[3,:] = mass of binaries disk_bins_bhbh[13,:] = ecc of binary around com disk_bins_bhbh[18,:] = orb. ecc of binary com around SMBH Keplerian orbital velocity of the bin c.o.m. around SMBH: v_bin,i= sqrt(GM_SMBH/R_bin,i_com)= c/sqrt(R_bin,i_com)

  3. Calculate the binary orbital time and N_orbits/timestep For example, since T_orb =2pi sqrt{bin,orb a}^3/GM_smbh) and {bin,orb a}^3/GM_smbh = (10^3r_g)^3/GM_smbh = 10^9 ({bin,orb a}/10^3r_g)^3 (GM_smbh/c^2)^3/GM_smbh

    = 10^9 ({bin,orb a}/10^3r_g)^3 (G M_smbh/c^3)^2

    So, .. math:

    T_{orb}
= 2\pi 10^{4.5} (R_{bin,orb a}/10^3r_g)^{3/2} GM_{smbh}/c^3
= 2\pi 10^{4.5} (R_{bin,orb a}/10^3r_g)^{3/2} (6.7e-11*2e38/(3e8)^3)
= 2\pi 10^{4.5} (R_{bin,orb a}/10^3r_g)^{3/2} (13.6e27/27e24)
= \pi 10^{7.5}  (R_{bin,orb a}/10^3r_g)^{3/2}
~ 3.15 yr (R_{bin,orb a}/10^3r_g)^3/2 (M_smbh/10^8Msun)

    i.e. Orbit~3.15yr at 10^3r_g around a 10^8M_{sun} SMBH. Therefore in a timestep=1.e4yr, a binary at 10^3r_g orbits the SMBH N_orbit/timestep =3,000 times.

  4. Calculate binding energy of bins = [GM1M2/sep_bin1, GMiMi+1,sep_bin2, ….] where sep_bin1 is in meters and M1,M2 are binary mass components in kg.

  5. Find those single BH with e>e_crit and their associated semi-major axes a_ecc =[a_ecc1, a_ecc2, ..] and masses m_ecc =[m_ecc1,m_ecc2, ..] and calculate their average velocities v_ecc = [GM_smbh/a_ecc1, GM_smbh/a_ecc2,…]

  6. Where (1-ecc_i)*a_ecc_i < R_bin_j_com < (1+ecc_i)*a_ecc_i, interaction possible

  7. Among candidate encounters, calculate relative velocity of encounter.

    \(v_{peri,i}=\sqrt(Gm_{ecc,i}/a_{ecc,i}[1+ecc,i/1-ecc,i])\) \(v_{apo,i} =\sqrt(Gm_{ecc,i}/a_{ecc,i}[1-ecc,i/1+ecc,i])\) \(v_{ecc,i} =\sqrt(GM/a_{ecc_i})\) ..average Keplerian vel.

    \(v_{rel} = abs(v_{bin,i} - v_{ecc,i})\)

  8. Calculate relative K.E. of tertiary, (1/2)m_ecc_i*v_rel_^2

  9. Compare binding en of binary to K.E. of tertiary.
    Critical velocity for ionization of binary is v_crit, given by:

    :math:`v_{crit} = sqrt(GM_1M_2(M_1+M_2+M_3)/M_3(M_1+M_2)a_{bin})

    If binary is hard ie GM_1M_2/a_bin > m3v_rel^2 then:
    harden binary

    a_bin -> a_bin -da_bin and

    new binary eccentricity

    e_bin -> e_bin + de

    and give +da_bin worth of binding energy (GM_bin/(a_bin -da_bin) - GM_bin/a_bin) to extra eccentricity ecc_i and a_ecc,i of m_ecc,i. Say average en of encounter is de=0.1 (10%) then binary a_bin shrinks by 10%, ecc_bin is pumped by 10% And a_ecc_i shrinks by 10% and ecc_i also shrinks by 10%

    If binary is soft ie GM_bin/a_bin <m3v_rel^2 then:
    if v_rel (effectively v_infty) > v_crit
    ionize binary

    update singleton array with 2 new BH with orbital eccentricity e_crit+de remove binary from binary array

    else if v_rel < v_crit
    soften binary

    a_bin -> a_bin + da_bin and

    new binary eccentricity

    e_bin -> e_bin + de

    and remove -da_bin worth of binary energy from eccentricity of m3.

Note1: Will need to test binary eccentricity each timestep.

If bin_ecc> some value (0.9), check for da_bin due to GW bremsstrahlung at pericenter.

  1. As 4, except now include interactions between binaries and circularized BH. This should give us primarily

    hardening encounters as in Leigh+2018, since the v_rel is likely to be small for more binaries.

Given array of binaries at locations [a_bbh1,a_bbh2] with binary semi-major axes [a_bin1,a_bin2,…] and binary eccentricities [e_bin1,e_bin2,…], find all the single BH at locations a_i that within timestep

either pass between a_i(1-e_i)< a_bbh1 <a_i(1+e_i)

Calculate velocity of encounter compared to a_bin. If binary is hard ie GM1M2/a_bin > m3v_rel^2 then:

harden binary to a_bin = a_bin -da_bin and new binary eccentricity e_bin = e_bin + de around com and new binary orb eccentricity e_orb_com = e_orb_com + de and now give da_bin worth of binding energy to extra eccentricity of m3.

If binary is soft ie GM_bin/a_bin <m3v_rel^2 then:

soften binary to a_bin = a_bin + da_bin and new binary eccentricity e_bin = e_bin + de and take da_bin worth of binary energy from eccentricity of m3.

If binary is unbound ie GM_bin/a_bin << m3v_rel^2 then:

remove binary from binary array add binary components m1,m2 back to singleton arrays with new orbital eccentricities e_1,e_2 from energy of encounter. Equipartition energy so m1v1^2 =m2 v_2^2 and generate new individual orbital eccentricities e1=v1/v_kep_circ and e_2=v_2/v_kep_circ Take energy put into destroying binary from orb. eccentricity of m3.

mcfacts.physics.dynamics.circular_binaries_encounters_ecc_prograde_star(smbh_mass, disk_star_pro_orbs_a, disk_star_pro_masses, disk_star_pro_orbs_ecc, disk_star_pro_id_nums, bin_mass_1, bin_mass_2, bin_orb_a, bin_sep, bin_ecc, bin_orb_ecc, bin_id_nums, rstar_rhill_exponent, timestep_duration_yr, disk_bh_pro_orb_ecc_crit, delta_energy_strong, disk_radius_outer, r_g_in_meters)

“Adjust orb eccentricities due to encounters between BBH and eccentric single BHs

Return array of modified binary BH separations and eccentricities perturbed by encounters within f*R_Hill, for eccentric singleton population, where f is some fraction/multiple of Hill sphere radius R_H Right now assume f=1.

Parameters:

  • smbh_mass (float) – Mass [M_sun] of supermassive black hole

  • disk_bh_pro_orbs_a (numpy.ndarray) – Orbital semi-major axes [r_{g,SMBH}] of prograde singleton BH at start of a timestep (math:r_g=GM_{SMBH}/c^2) with float type

  • disk_bh_pro_masses (numpy.ndarray) – Masses [M_sun] of prograde singleton BH at start of timestep with float type

  • disk_bh_pro_orbs_ecc (numpy.ndarray) – Orbital eccentricity [unitless] of singleton prograde BH with float type

  • timestep_duration_yr (float) – Length of timestep [yr]

  • disk_bh_pro_orb_ecc_crit (float) – Critical orbital eccentricity [unitless] below which orbit is close enough to circularize

  • delta_energy_strong (float) – Average energy change [units??] per strong encounter

  • disk_bins_bhbh (numpy.ndarray) – [21, bindex] mixed array containing properties of binary BBH, see add_to_binary_array function for complete description

  • disk_radius_outer (float) – Outer radius of the inner disk (Rg)

  • r_g_in_meters (float) – Gravitational radius of the SMBH in meters

Returns:

disk_bins_bhbh – [21, bindex] mixed array, updated version of input after dynamical perturbations

Return type:

numpy.ndarray

Notes

Logic:

  1. Find number of binaries in this timestep given by bindex

  2. Find the binary center of mass (c.o.m.) and corresponding orbital velocities & binary total masses. disk_bins_bhbh[9,:] = bin c.o.m. = [R_bin1_orb_a,R_bin2_orb_a,…]. These are the orbital radii of the bins. disk_bins_bhbh[8,;] = bin_separation =[a_bin1,a_bin2,…] disk_bins_bhbh[2,:]+disk_bins_bhbh[3,:] = mass of binaries disk_bins_bhbh[13,:] = ecc of binary around com disk_bins_bhbh[18,:] = orb. ecc of binary com around SMBH Keplerian orbital velocity of the bin c.o.m. around SMBH: v_bin,i= sqrt(GM_SMBH/R_bin,i_com)= c/sqrt(R_bin,i_com)

  3. Calculate the binary orbital time and N_orbits/timestep For example, since T_orb =2pi sqrt{bin,orb a}^3/GM_smbh) and {bin,orb a}^3/GM_smbh = (10^3r_g)^3/GM_smbh = 10^9 ({bin,orb a}/10^3r_g)^3 (GM_smbh/c^2)^3/GM_smbh

    = 10^9 ({bin,orb a}/10^3r_g)^3 (G M_smbh/c^3)^2

    So, .. math:

    T_{orb}
= 2\pi 10^{4.5} (R_{bin,orb a}/10^3r_g)^{3/2} GM_{smbh}/c^3
= 2\pi 10^{4.5} (R_{bin,orb a}/10^3r_g)^{3/2} (6.7e-11*2e38/(3e8)^3)
= 2\pi 10^{4.5} (R_{bin,orb a}/10^3r_g)^{3/2} (13.6e27/27e24)
= \pi 10^{7.5}  (R_{bin,orb a}/10^3r_g)^{3/2}
~ 3.15 yr (R_{bin,orb a}/10^3r_g)^3/2 (M_smbh/10^8Msun)

    i.e. Orbit~3.15yr at 10^3r_g around a 10^8M_{sun} SMBH. Therefore in a timestep=1.e4yr, a binary at 10^3r_g orbits the SMBH N_orbit/timestep =3,000 times.

  4. Calculate binding energy of bins = [GM1M2/sep_bin1, GMiMi+1,sep_bin2, ….] where sep_bin1 is in meters and M1,M2 are binary mass components in kg.

  5. Find those single BH with e>e_crit and their associated semi-major axes a_ecc =[a_ecc1, a_ecc2, ..] and masses m_ecc =[m_ecc1,m_ecc2, ..] and calculate their average velocities v_ecc = [GM_smbh/a_ecc1, GM_smbh/a_ecc2,…]

  6. Where (1-ecc_i)*a_ecc_i < R_bin_j_com < (1+ecc_i)*a_ecc_i, interaction possible

  7. Among candidate encounters, calculate relative velocity of encounter.

    \(v_{peri,i}=\sqrt(Gm_{ecc,i}/a_{ecc,i}[1+ecc,i/1-ecc,i])\) \(v_{apo,i} =\sqrt(Gm_{ecc,i}/a_{ecc,i}[1-ecc,i/1+ecc,i])\) \(v_{ecc,i} =\sqrt(GM/a_{ecc_i})\) ..average Keplerian vel.

    \(v_{rel} = abs(v_{bin,i} - v_{ecc,i})\)

  8. Calculate relative K.E. of tertiary, (1/2)m_ecc_i*v_rel_^2

  9. Compare binding en of binary to K.E. of tertiary.
    Critical velocity for ionization of binary is v_crit, given by:

    :math:`v_{crit} = sqrt(GM_1M_2(M_1+M_2+M_3)/M_3(M_1+M_2)a_{bin})

    If binary is hard ie GM_1M_2/a_bin > m3v_rel^2 then:
    harden binary

    a_bin -> a_bin -da_bin and

    new binary eccentricity

    e_bin -> e_bin + de

    and give +da_bin worth of binding energy (GM_bin/(a_bin -da_bin) - GM_bin/a_bin) to extra eccentricity ecc_i and a_ecc,i of m_ecc,i. Say average en of encounter is de=0.1 (10%) then binary a_bin shrinks by 10%, ecc_bin is pumped by 10% And a_ecc_i shrinks by 10% and ecc_i also shrinks by 10%

    If binary is soft ie GM_bin/a_bin <m3v_rel^2 then:
    if v_rel (effectively v_infty) > v_crit
    ionize binary

    update singleton array with 2 new BH with orbital eccentricity e_crit+de remove binary from binary array

    else if v_rel < v_crit
    soften binary

    a_bin -> a_bin + da_bin and

    new binary eccentricity

    e_bin -> e_bin + de

    and remove -da_bin worth of binary energy from eccentricity of m3.

Note1: Will need to test binary eccentricity each timestep.

If bin_ecc> some value (0.9), check for da_bin due to GW bremsstrahlung at pericenter.

  1. As 4, except now include interactions between binaries and circularized BH. This should give us primarily

    hardening encounters as in Leigh+2018, since the v_rel is likely to be small for more binaries.

Given array of binaries at locations [a_bbh1,a_bbh2] with binary semi-major axes [a_bin1,a_bin2,…] and binary eccentricities [e_bin1,e_bin2,…], find all the single BH at locations a_i that within timestep

either pass between a_i(1-e_i)< a_bbh1 <a_i(1+e_i)

Calculate velocity of encounter compared to a_bin. If binary is hard ie GM1M2/a_bin > m3v_rel^2 then:

harden binary to a_bin = a_bin -da_bin and new binary eccentricity e_bin = e_bin + de around com and new binary orb eccentricity e_orb_com = e_orb_com + de and now give da_bin worth of binding energy to extra eccentricity of m3.

If binary is soft ie GM_bin/a_bin <m3v_rel^2 then:

soften binary to a_bin = a_bin + da_bin and new binary eccentricity e_bin = e_bin + de and take da_bin worth of binary energy from eccentricity of m3.

If binary is unbound ie GM_bin/a_bin << m3v_rel^2 then:

remove binary from binary array add binary components m1,m2 back to singleton arrays with new orbital eccentricities e_1,e_2 from energy of encounter. Equipartition energy so m1v1^2 =m2 v_2^2 and generate new individual orbital eccentricities e1=v1/v_kep_circ and e_2=v_2/v_kep_circ Take energy put into destroying binary from orb. eccentricity of m3.

mcfacts.physics.dynamics.circular_singles_encounters_prograde(smbh_mass, disk_bh_pro_orbs_a, disk_bh_pro_masses, disk_bh_pro_orbs_ecc, timestep_duration_yr, disk_bh_pro_orb_ecc_crit, delta_energy_strong, disk_radius_outer, rng_here=RandomState(MT19937) at 0x7602873A3C40)

“Adjust orb ecc due to encounters between 2 single circ pro BH

Parameters:

  • smbh_mass (float) – Mass [M_sun] of supermassive black hole

  • disk_bh_pro_orbs_a (numpy.ndarray) – Orbital semi-major axes [r_{g,SMBH}] of prograde singleton BH at start of a timestep (math:r_g=GM_{SMBH}/c^2) with float type

  • disk_bh_pro_masses (numpy.ndarray) – Masses [M_sun] of prograde singleton BH at start of timestep with float type

  • disk_bh_pro_orbs_ecc (numpy.ndarray) – Orbital eccentricity [unitless] of singleton prograde BH with float type

  • timestep_duration_yr (float) – Length of timestep [yr]

  • disk_bh_pro_orb_ecc_crit (float) – Critical orbital eccentricity [unitless] below which orbit is close enough to circularize

  • delta_energy_strong (float) – Average energy change [units??] per strong encounter

  • disk_radius_outer (float) – Outer radius of the inner disk (Rg)

Returns:

  • disk_bh_pro_orbs_a (numpy.ndarray) – Updated BH semi-major axis [r_{g,SMBH}] perturbed by dynamics with float type

  • disk_bh_pro_orbs_ecc (numpy.ndarray) – Updated BH orbital eccentricities [unitless] perturbed by dynamics with float type

Notes

Return array of modified singleton BH orbital eccentricities perturbed by encounters within \(f*R_{Hill}\), where f is some fraction/multiple of Hill sphere radius R_H

Assume encounters between damped BH (e<e_crit) and undamped BH (e>e_crit) are the only important ones for now. Since the e<e_crit population is the most likely BBH merger source.

1, find those orbiters with e<e_crit and their

associated semi-major axes a_circ =[a_circ1, a_circ2, ..] and masses m_circ =[m_circ1,m_circ2, ..].

2, calculate orbital timescales for a_circ1 and a_i and N_orbits/timestep.

For example, since \(T_orb =2\pi \sqrt(a^3/GM_{smbh})\) and .. math:: a^3/GM_{smbh} = (10^3r_g)^3/GM_{smbh} = 10^9 (a/10^3r_g)^3 (GM_{smbh}/c^2)^3/GM_{smbh}

= 10^9 (a/10^3r_g)^3 (G M_{smbh}/c^3)^2

So .. math:

T_orb   = 2\pi 10^{4.5} (a/10^3r_g)^{3/2} GM_{smbh}/c^3 \
        = 2\pi 10^{4.5} (a/10^3r_g)^{3/2} (6.7e-11*2e38/(3e8)^3) \
        = 2\pi 10^{4.5} (a/10^3r_g)^{3/2} (13.6e27/27e24) \
        = \pi 10^{7.5}  (a/10^3r_g)^{3/2} \
        ~ 3yr (a/10^3r_g)^3/2 (M_{smbh}/10^8M_{sun}) \

i.e. Orbit~3yr at 10^3r_g around a 10^8M_{sun} SMBH. Therefore in a timestep=1.e4yr, a BH at 10^3r_g orbits the SMBH N_orbit/timestep =3,000 times.

3, among population of orbiters with e>e_crit,

find those orbiters (a_i,e_i) where a_i*(1-e_i)< a_circ1,j <a_i*(1-e_i) for all members a_circ1,j of the circularized population so we can test for possible interactions.

4, calculate mutual Hill sphere R_H of candidate binary (a_circ1,j ,a_i).

5, calculate ratio of 2R_H of binary to size of circular orbit, or (2R_H/2pi a_circ1,j)

Hill sphere possible on both crossing inwards and outwards once per orbit, so 2xHill sphere =4R_H worth of circular orbit will have possible encounter. Thus, (4R_H/2pi a_circ1)= odds that a_circ1 is in the region of cross-over per orbit. For example, for BH at a_circ1 = 1e3r_g,

\[R_h = a_{circ1}*(m_{circ1} + m_i/3M_{smbh})^1/3\]
\[= 0.004a_{circ1} (m_{circ1}/10M_{sun})^1/3 (m_i/10M_{sun})^1/3 (M_{smbh}/1e8M_{sun})^-1/3\]
then

ratio (4R_H/2pi a_circ1) = 0.008/pi ~ 0.0026 (ie around 1/400 odds that BH at a_circ1 is in either area of crossing)

6, calculate number of orbits of a_i in 1 timestep.

If e.g. N_orb(a_i)/timestep = 200 orbits per timestep of 10kyr, then probability of encounter = (200orbits/timestep)*(4R_H/2pi a_circ1) ~ 0.5,

or 50% odds of an encounter on this timestep between (a_circ1,j , a_i).

If probability > 1, set probability = 1.

7, draw a random number from the uniform [0,1] distribution and

if rng < probability of encounter, there is an encounter during the timestep if rng > probability of encounter, there is no encounter during the timestep

8, if encounter:

Take energy (de) from high ecc. a_i and give energy (de) to a_circ1,j de is average fractional energy change per encounter.

So, a_circ1,j ->(1+de)a_circ1,j.

e_circ1,j ->(crit_ecc + de)

and

a_i ->(1-de)a_i e_i ->(1-de)e_i

Could be that average energy in gas-free cluster case is assume average energy transfer = 20% perturbation (from Sigurdsson & Phinney 1993).

Further notes for self: sigma_ecc = sqrt(ecc^2 + incl^2)v_kep so if incl=0 deg (for now) En of ecc. interloper = 1/2 m_i sigma_ecc^2.

Note: Can also use above logic for binary encounters except use binary binding energy instead.

or later could try

Deflection angle defl = tan (defl) = dV_perp/V = 2GM/bV^2 kg^-1 m^3 s^-2 kg / m (m s^-1)^2

so \(de/e =2GM/bV^2 = 2 G M_{bin}/0.5R_{hill}*\sigma^2\) and \(R_hill = a_{circ1}*(M_{bin}/3M_{smbh})^1/3 and \sigma^2 =ecc^2*v_{kep}^2\) So \(de/e = 4GM_{bin}/a_{circ1}(M_{bin}/3M_{smbh})^1/3 ecc^2 v_{kep}^2\) and \(v_{kep} = \sqrt(GM_{smbh}/a_i)\) So \(de/e = 4GM_{bin}^{2/3}M_{smbh}^1/3 a_i/a_{circ1} ecc^2 GM_{smbh} = 4(M_{bin}/M_{smbh})^{2/3} (a_i/a_{circ1})(1/ecc^2) where :math:`V_{rel} = \sigma\) say and \(b=R_H = a_{circ1} (q/3)^{1/3}\) So \(defl = 2GM/ a_{circ1}(q/3)^2/3 ecc^2 10^14 (m/s)^2 (R/10^3r_g)^-1\)

\(= 2 6.7e-11 2.e31/\)

!!Note: when doing this for binaries.

Calculate velocity of encounter compared to a_bin. If binary is hard ie GM_bin/a_bin > m3v_rel^2 then: harden binary

a_bin -> a_bin -da_bin and

new binary eccentricity

e_bin -> e_bin + de

and give da_bin worth of binding energy to extra eccentricity of m3. If binary is soft ie GM_bin/a_bin <m3v_rel^2 then: soften binary

a_bin -> a_bin + da_bin and

new binary eccentricity

e_bin -> e_bin + de

and remove da_bin worth of binary energy from eccentricity of m3.

mcfacts.physics.dynamics.circular_singles_encounters_prograde_star_bh(smbh_mass, disk_star_pro_orbs_a, disk_star_pro_masses, disk_star_pro_radius, disk_star_pro_orbs_ecc, disk_star_pro_id_nums, rstar_rhill_exponent, disk_bh_pro_orbs_a, disk_bh_pro_masses, disk_bh_pro_orbs_ecc, disk_bh_pro_id_nums, timestep_duration_yr, disk_bh_pro_orb_ecc_crit, delta_energy_strong_mu, delta_energy_strong_sigma, disk_radius_outer)

“Adjust orb ecc due to encounters between single circ star and single ecc black hole

Parameters:

  • smbh_mass (float) – Mass [M_sun] of supermassive black hole

  • disk_bh_pro_orbs_a (numpy.ndarray) – Orbital semi-major axes [r_{g,SMBH}] of prograde singleton star at start of a timestep (math:r_g=GM_{SMBH}/c^2) with float type

  • disk_bh_pro_masses (numpy.ndarray) – Masses [M_sun] of prograde singleton star at start of timestep with float type

  • disk_star_pro_radius (numpy.ndarray) – Radii [Rsun] of prograde singleton star at start of timestep with :obj: float type

  • disk_bh_pro_orbs_ecc (numpy.ndarray) – Orbital eccentricity [unitless] of singleton prograde star with float type

  • disk_star_pro_id_nums (numpy.ndarray) – ID numbers of singleton prograde stars

  • rstar_rhill_exponent (float) – Exponent for the ratio of R_star / R_Hill. Default is 2

  • timestep_duration_yr (float) – Length of timestep [yr]

  • disk_bh_pro_orb_ecc_crit (float) – Critical orbital eccentricity [unitless] below which orbit is close enough to circularize

  • delta_energy_strong (float) – Average energy change [units??] per strong encounter

Returns:

  • disk_star_pro_orbs_a (numpy.ndarray) – Updated stars semi-major axis [r_{g,SMBH}] perturbed by dynamics with float type

  • disk_star_pro_orbs_ecc (numpy.ndarray) – Updated stars orbital eccentricities [unitless] perturbed by dynamics with float type

  • disk_star_pro_id_nums_touch (numpy.ndarray) – ID numbers of stars that will touch each other

  • disk_bh_pro_orbs_a (numpy.ndarray) – Updated BH semi-major axis [r_{g,SMBH}] perturbed by dynamics with float type

  • disk_bh_pro_orbs_ecc (numpy.ndarray) – Updated BH orbital eccentricities [unitless] perturbed by dynamics with float type

Notes

Return array of modified singleton star orbital eccentricities perturbed by encounters within \(f*R_{Hill}\), where f is some fraction/multiple of Hill sphere radius R_H

Assume encounters between damped star (e<e_crit) and undamped star (e>e_crit) are the only important ones for now. Since the e<e_crit population is the most likely BBH merger source.

1, find those orbiters with e<e_crit and their

associated semi-major axes a_circ =[a_circ1, a_circ2, ..] and masses m_circ =[m_circ1,m_circ2, ..].

2, calculate orbital timescales for a_circ1 and a_i and N_orbits/timestep.

For example, since \(T_orb =2\pi \sqrt(a^3/GM_{smbh})\) and .. math:: a^3/GM_{smbh} = (10^3r_g)^3/GM_{smbh} = 10^9 (a/10^3r_g)^3 (GM_{smbh}/c^2)^3/GM_{smbh}

= 10^9 (a/10^3r_g)^3 (G M_{smbh}/c^3)^2

So .. math:

T_orb   = 2\pi 10^{4.5} (a/10^3r_g)^{3/2} GM_{smbh}/c^3 \
        = 2\pi 10^{4.5} (a/10^3r_g)^{3/2} (6.7e-11*2e38/(3e8)^3) \
        = 2\pi 10^{4.5} (a/10^3r_g)^{3/2} (13.6e27/27e24) \
        = \pi 10^{7.5}  (a/10^3r_g)^{3/2} \
        ~ 3yr (a/10^3r_g)^3/2 (M_{smbh}/10^8M_{sun}) \

i.e. Orbit~3yr at 10^3r_g around a 10^8M_{sun} SMBH. Therefore in a timestep=1.e4yr, a BH at 10^3r_g orbits the SMBH N_orbit/timestep =3,000 times.

3, among population of orbiters with e>e_crit,

find those orbiters (a_i,e_i) where a_i*(1-e_i)< a_circ1,j <a_i*(1-e_i) for all members a_circ1,j of the circularized population so we can test for possible interactions.

4, calculate mutual Hill sphere R_H of candidate binary (a_circ1,j ,a_i).

5, calculate ratio of 2R_H of binary to size of circular orbit, or (2R_H/2pi a_circ1,j)

Hill sphere possible on both crossing inwards and outwards once per orbit, so 2xHill sphere =4R_H worth of circular orbit will have possible encounter. Thus, (4R_H/2pi a_circ1)= odds that a_circ1 is in the region of cross-over per orbit. For example, for BH at a_circ1 = 1e3r_g,

\[R_h = a_{circ1}*(m_{circ1} + m_i/3M_{smbh})^1/3\]
\[= 0.004a_{circ1} (m_{circ1}/10M_{sun})^1/3 (m_i/10M_{sun})^1/3 (M_{smbh}/1e8M_{sun})^-1/3\]
then

ratio (4R_H/2pi a_circ1) = 0.008/pi ~ 0.0026 (ie around 1/400 odds that BH at a_circ1 is in either area of crossing)

6, calculate number of orbits of a_i in 1 timestep.

If e.g. N_orb(a_i)/timestep = 200 orbits per timestep of 10kyr, then probability of encounter = (200orbits/timestep)*(4R_H/2pi a_circ1) ~ 0.5,

or 50% odds of an encounter on this timestep between (a_circ1,j , a_i).

If probability > 1, set probability = 1.

7, draw a random number from the uniform [0,1] distribution and

if rng < probability of encounter, there is an encounter during the timestep if rng > probability of encounter, there is no encounter during the timestep

8, if encounter:

Take energy (de) from high ecc. a_i and give energy (de) to a_circ1,j de is average fractional energy change per encounter.

So, a_circ1,j ->(1+de)a_circ1,j.

e_circ1,j ->(crit_ecc + de)

and

a_i ->(1-de)a_i e_i ->(1-de)e_i

Could be that average energy in gas-free cluster case is assume average energy transfer = 20% perturbation (from Sigurdsson & Phinney 1993).

Further notes for self: sigma_ecc = sqrt(ecc^2 + incl^2)v_kep so if incl=0 deg (for now) En of ecc. interloper = 1/2 m_i sigma_ecc^2.

Note: Can also use above logic for binary encounters except use binary binding energy instead.

or later could try

Deflection angle defl = tan (defl) = dV_perp/V = 2GM/bV^2 kg^-1 m^3 s^-2 kg / m (m s^-1)^2

so \(de/e =2GM/bV^2 = 2 G M_{bin}/0.5R_{hill}*\sigma^2\) and \(R_hill = a_{circ1}*(M_{bin}/3M_{smbh})^1/3 and \sigma^2 =ecc^2*v_{kep}^2\) So \(de/e = 4GM_{bin}/a_{circ1}(M_{bin}/3M_{smbh})^1/3 ecc^2 v_{kep}^2\) and \(v_{kep} = \sqrt(GM_{smbh}/a_i)\) So \(de/e = 4GM_{bin}^{2/3}M_{smbh}^1/3 a_i/a_{circ1} ecc^2 GM_{smbh} = 4(M_{bin}/M_{smbh})^{2/3} (a_i/a_{circ1})(1/ecc^2) where :math:`V_{rel} = \sigma\) say and \(b=R_H = a_{circ1} (q/3)^{1/3}\) So \(defl = 2GM/ a_{circ1}(q/3)^2/3 ecc^2 10^14 (m/s)^2 (R/10^3r_g)^-1\)

\(= 2 6.7e-11 2.e31/\)

!!Note: when doing this for binaries.

Calculate velocity of encounter compared to a_bin. If binary is hard ie GM_bin/a_bin > m3v_rel^2 then: harden binary

a_bin -> a_bin -da_bin and

new binary eccentricity

e_bin -> e_bin + de

and give da_bin worth of binding energy to extra eccentricity of m3. If binary is soft ie GM_bin/a_bin <m3v_rel^2 then: soften binary

a_bin -> a_bin + da_bin and

new binary eccentricity

e_bin -> e_bin + de

and remove da_bin worth of binary energy from eccentricity of m3.

mcfacts.physics.dynamics.circular_singles_encounters_prograde_stars(smbh_mass, disk_star_pro_orbs_a, disk_star_pro_masses, disk_star_pro_radius, disk_star_pro_orbs_ecc, disk_star_pro_id_nums, rstar_rhill_exponent, timestep_duration_yr, disk_bh_pro_orb_ecc_crit, delta_energy_strong_mu, delta_energy_strong_sigma, disk_radius_outer, rng_here=RandomState(MT19937) at 0x7602873A3C40, fast_cube=False)

“Adjust orb ecc due to encounters between 2 single circ pro stars

Parameters:

  • smbh_mass (float) – Mass [M_sun] of supermassive black hole

  • disk_bh_pro_orbs_a (numpy.ndarray) – Orbital semi-major axes [r_{g,SMBH}] of prograde singleton star at start of a timestep (math:r_g=GM_{SMBH}/c^2) with float type

  • disk_bh_pro_masses (numpy.ndarray) – Masses [M_sun] of prograde singleton star at start of timestep with float type

  • disk_star_pro_radius (numpy.ndarray) – Radii [Rsun] of prograde singleton star at start of timestep with :obj: float type

  • disk_bh_pro_orbs_ecc (numpy.ndarray) – Orbital eccentricity [unitless] of singleton prograde star with float type

  • disk_star_pro_id_nums (numpy.ndarray) – ID numbers of singleton prograde stars

  • rstar_rhill_exponent (float) – Exponent for the ratio of R_star / R_Hill. Default is 2

  • timestep_duration_yr (float) – Length of timestep [yr]

  • disk_bh_pro_orb_ecc_crit (float) – Critical orbital eccentricity [unitless] below which orbit is close enough to circularize

  • delta_energy_strong_mu (float) – Average energy change [units??] per strong encounter

  • delta_energy_strong_sigma (float) – Standard deviation of average energy change per strong encounter

Returns:

  • disk_star_pro_orbs_a (numpy.ndarray) – Updated BH semi-major axis [r_{g,SMBH}] perturbed by dynamics with float type

  • disk_star_pro_orbs_ecc (numpy.ndarray) – Updated BH orbital eccentricities [unitless] perturbed by dynamics with float type

  • disk_star_pro_id_nums_touch (numpy.ndarray) – ID numbers of stars that will touch each other

Notes

Return array of modified singleton star orbital eccentricities perturbed by encounters within \(f*R_{Hill}\), where f is some fraction/multiple of Hill sphere radius R_H

Assume encounters between damped star (e<e_crit) and undamped star (e>e_crit) are the only important ones for now. Since the e<e_crit population is the most likely BBH merger source.

1, find those orbiters with e<e_crit and their

associated semi-major axes a_circ =[a_circ1, a_circ2, ..] and masses m_circ =[m_circ1,m_circ2, ..].

2, calculate orbital timescales for a_circ1 and a_i and N_orbits/timestep.

For example, since \(T_orb =2\pi \sqrt(a^3/GM_{smbh})\) and .. math:: a^3/GM_{smbh} = (10^3r_g)^3/GM_{smbh} = 10^9 (a/10^3r_g)^3 (GM_{smbh}/c^2)^3/GM_{smbh}

= 10^9 (a/10^3r_g)^3 (G M_{smbh}/c^3)^2

So .. math:

T_orb   = 2\pi 10^{4.5} (a/10^3r_g)^{3/2} GM_{smbh}/c^3 \
        = 2\pi 10^{4.5} (a/10^3r_g)^{3/2} (6.7e-11*2e38/(3e8)^3) \
        = 2\pi 10^{4.5} (a/10^3r_g)^{3/2} (13.6e27/27e24) \
        = \pi 10^{7.5}  (a/10^3r_g)^{3/2} \
        ~ 3yr (a/10^3r_g)^3/2 (M_{smbh}/10^8M_{sun}) \

i.e. Orbit~3yr at 10^3r_g around a 10^8M_{sun} SMBH. Therefore in a timestep=1.e4yr, a BH at 10^3r_g orbits the SMBH N_orbit/timestep =3,000 times.

3, among population of orbiters with e>e_crit,

find those orbiters (a_i,e_i) where a_i*(1-e_i)< a_circ1,j <a_i*(1-e_i) for all members a_circ1,j of the circularized population so we can test for possible interactions.

4, calculate mutual Hill sphere R_H of candidate binary (a_circ1,j ,a_i).

5, calculate ratio of 2R_H of binary to size of circular orbit, or (2R_H/2pi a_circ1,j)

Hill sphere possible on both crossing inwards and outwards once per orbit, so 2xHill sphere =4R_H worth of circular orbit will have possible encounter. Thus, (4R_H/2pi a_circ1)= odds that a_circ1 is in the region of cross-over per orbit. For example, for BH at a_circ1 = 1e3r_g,

\[R_h = a_{circ1}*(m_{circ1} + m_i/3M_{smbh})^1/3\]
\[= 0.004a_{circ1} (m_{circ1}/10M_{sun})^1/3 (m_i/10M_{sun})^1/3 (M_{smbh}/1e8M_{sun})^-1/3\]
then

ratio (4R_H/2pi a_circ1) = 0.008/pi ~ 0.0026 (ie around 1/400 odds that BH at a_circ1 is in either area of crossing)

6, calculate number of orbits of a_i in 1 timestep.

If e.g. N_orb(a_i)/timestep = 200 orbits per timestep of 10kyr, then probability of encounter = (200orbits/timestep)*(4R_H/2pi a_circ1) ~ 0.5,

or 50% odds of an encounter on this timestep between (a_circ1,j , a_i).

If probability > 1, set probability = 1.

7, draw a random number from the uniform [0,1] distribution and

if rng < probability of encounter, there is an encounter during the timestep if rng > probability of encounter, there is no encounter during the timestep

8, if encounter:

Take energy (de) from high ecc. a_i and give energy (de) to a_circ1,j de is average fractional energy change per encounter.

So, a_circ1,j ->(1+de)a_circ1,j.

e_circ1,j ->(crit_ecc + de)

and

a_i ->(1-de)a_i e_i ->(1-de)e_i

Could be that average energy in gas-free cluster case is assume average energy transfer = 20% perturbation (from Sigurdsson & Phinney 1993).

Further notes for self: sigma_ecc = sqrt(ecc^2 + incl^2)v_kep so if incl=0 deg (for now) En of ecc. interloper = 1/2 m_i sigma_ecc^2.

Note: Can also use above logic for binary encounters except use binary binding energy instead.

or later could try

Deflection angle defl = tan (defl) = dV_perp/V = 2GM/bV^2 kg^-1 m^3 s^-2 kg / m (m s^-1)^2

so \(de/e =2GM/bV^2 = 2 G M_{bin}/0.5R_{hill}*\sigma^2\) and \(R_hill = a_{circ1}*(M_{bin}/3M_{smbh})^1/3 and \sigma^2 =ecc^2*v_{kep}^2\) So \(de/e = 4GM_{bin}/a_{circ1}(M_{bin}/3M_{smbh})^1/3 ecc^2 v_{kep}^2\) and \(v_{kep} = \sqrt(GM_{smbh}/a_i)\) So \(de/e = 4GM_{bin}^{2/3}M_{smbh}^1/3 a_i/a_{circ1} ecc^2 GM_{smbh} = 4(M_{bin}/M_{smbh})^{2/3} (a_i/a_{circ1})(1/ecc^2) where :math:`V_{rel} = \sigma\) say and \(b=R_H = a_{circ1} (q/3)^{1/3}\) So \(defl = 2GM/ a_{circ1}(q/3)^2/3 ecc^2 10^14 (m/s)^2 (R/10^3r_g)^-1\)

\(= 2 6.7e-11 2.e31/\)

!!Note: when doing this for binaries.

Calculate velocity of encounter compared to a_bin. If binary is hard ie GM_bin/a_bin > m3v_rel^2 then: harden binary

a_bin -> a_bin -da_bin and

new binary eccentricity

e_bin -> e_bin + de

and give da_bin worth of binding energy to extra eccentricity of m3. If binary is soft ie GM_bin/a_bin <m3v_rel^2 then: soften binary

a_bin -> a_bin + da_bin and

new binary eccentricity

e_bin -> e_bin + de

and remove da_bin worth of binary energy from eccentricity of m3.

mcfacts.physics.dynamics.components_from_EL(E, L, units='geometric', smbh_mass=100000000.0)

Calculates new orb_a and eccentricity from specific energy and specific angular momentum

Parameters:

  • E (float) – specific energy (per unit mass)

  • L (float) – specific angular momentum (per unit mass)

  • units (str, optional) – whether to use geometric units, by default ‘geometric’

  • smbh_mass (float, optional) – Mass [Msun] of the SMBH, by default 1e8

Returns:

new orb_a [r_{g,SMBH}] and ecc values

Return type:

orb_a, ecc

mcfacts.physics.dynamics.cubic_finite_step_root(x0, y0, OmegaS, sanity=False)

Determine allowed finite step size

Parameters:

  • x0 (float) – dimensionless x0 value

  • y0 (float) – dimensionless y0 value

  • OmegaS (float) – Orbital frequency [???]

  • sanity (bool, optional) – Switch, turns on extra print statements, by default False

Returns:

roots_x,roots_y – Found roots

Return type:

np.array

mcfacts.physics.dynamics.cubic_finite_step_root_cardano(x0, y0, OmegaS, sanity=False)

Optimized version to determine allowed finite step size using an analytic solution for the depressed cubic equation.

mcfacts.physics.dynamics.cubic_y_root(x0, y0, sanity=False)

Calculate the root of cubic function f(y) = x0*y^3 + 1.5*y - y0

Parameters:

  • x0 (float) – dimensionless x0 variable

  • y0 (float) – dimensionless y0 variable

  • sanity (bool, optional) – Switch, turns on extra print statements, by default False

Returns:

roots of the function

Return type:

roots

mcfacts.physics.dynamics.cubic_y_root_cardano(x0, y0, sanity=False)

Optimized version of cubic_y_root using an analytic solver. Solves the equation: x0*y^3 + 1.5*y - y0 = 0

mcfacts.physics.dynamics.encounters_new_orba_ecc(smbh_mass, orb_a_give, orb_a_take, mass_give, mass_take, ecc_give, ecc_take, radius_give, radius_take, id_num_give, id_num_take, delta_energy_strong, flag_obj_types, fast_cube=False)

Calculate new orb_a and ecc values for two objects that dynamically interact

Parameters:

  • smbh_mass (float) – Mass [M_sun] of the SMBH

  • orb_a_give (float) – Semi-major axis [r_{g,SMBH}] of the object donating energy (typically the eccentric)

  • orb_a_take (float) – Semi-major axis [r_{g,SMBH}] of the object accreting energy (typically the circular)

  • mass_give (float) – Mass [M_sun] of the object donating energy

  • mass_take (float) – Mass [M_sun] of the object accreting energy

  • ecc_give (float) – Eccentricity of the object donating energy

  • ecc_take (float) – Eccentricity of the object accreting energy

  • radius_give (float) – Radius [r_{g,SMBH}] of the object donating energy

  • radius_take (float) – Radius [r_{g,SMBH}] of the object accreting energy

  • id_num_give (int) – ID number of the object donating energy

  • id_num_take (int) – ID number of the object accreting energy

  • delta_energy_strong (float) – Average energy change per strong encounter

  • flag_obj_types (int) – Switch determining the type of interaction 0 : eccentric star - circular star 1 : eccentric black hole - circular star 2 : eccentric black hole - circular black hole 3 : eccentric star - eccentric star

Returns:

  • orb_a_give_final (float) – New semi-major axis [r_{g,SMBH}] of the object donating energy

  • orb_a_take_final (float) – New semi-major axis [r_{g,SMBH}] of the object accreting energy

  • ecc_give_final (float) – New eccentricity of the object donating energy

  • ecc_take_final (float) – New eccentricity of the object accreting energy

  • id_num_unbound (int) – ID number of object unbound from the disk (if any, otherwise None)

  • id_num_flipped_rotation (int) – ID number of object flipped from prograde to retrograde (if any, otherwise None)

mcfacts.physics.dynamics.transition_physical_as_EL(E1, L1, E2, L2, DeltaE, m1, m2, units='geometric', smbh_mass=100000000.0, sanity=False, fast_cube=False)

Calculates final energy and angular momentum states

Parameters:

  • E1 (float) – energy of object 1

  • L1 (float) – angular momentum of object 1

  • E2 (float) – energy of object 2

  • L2 (float) – angular momentum of object 2

  • DeltaE (float) – change in energy

  • m1 (float) – mass of object 1

  • m2 (float) – mass of object 2

  • units (str, optional) – Switch to control type of units, by default ‘geometric’

  • smbh_mass (float, optional) – SMBH mass, by default 1e8

  • sanity (bool, optional) – Switch to turn on extra print statements, by default True