In [1]:
# Code to evaluate the analytic toy model for divergence rate in Section 5.2
In [2]:
%matplotlib inline
In [3]:
import matplotlib as mpl; mpl.rcParams['savefig.dpi'] = 144
import matplotlib.pyplot as plt
import numpy as np

Support Functions

In [4]:
# Compute Local Number Density
def compute_number_density(R,w,Sigma0,R0,Npart):
    """
    @param  R, w   [AU]
    @param  Sigma0 [kg/AU2]
    @param  R0     [AU]
    @return nbin   [1/AU3]
    """
    Mbin = Sigma0 * R0 * twopi * w
    Mpart = 5.0 / Npart # M_Earth
    Nbin = Mbin / Mpart
    #@debug
    #print Nbin * len(R)
    Rlo = R - 0.5 * w
    Rhi = R + 0.5 * w
    Abin = twopi * ( Rhi**2.0 - Rlo**2.0 )
    Vbin = aspect * R * Abin
    #@debug
    #print Nbin
    #print Vbin
    n3d = Nbin / Vbin
    n2d = Nbin / Abin
    #@debug
    #print np.sum(nbin * Vbin)
    #print nbin/Vbin
    return n3d, n2d

# Relative Velocity between two Objects on Circular(!) Keplerian Orbits
def compute_vrel(R,w,GM):
    """
    @param  GM   [km3/s2]
    @param  R    [km]
    @param  w    [km]
    @return vrel [km/s]
    """
    vrel = np.sqrt(GM/R) - np.sqrt(GM/(R+w))
    return vrel

# Hill Radius
def compute_hill_radius(R,m):
    """
    @param  R      [SomeUnit]
    @return r_hill [SomeUnit]
    """
    r_hill = R * (2.0 * m / 3.0 / Msun)**(1.0/3.0)
    return r_hill

Constants

In [5]:
# Masses [kg]
m_sun = 1.9891e30
m_earth = 5.972e24
m_jup = 1.898e27
m_sat = 568.3e24

# Constants
au = 149597871.0 # km
r2d = 360.0 / np.pi / 2.0

# More Constants
twopi = 2.0 * np.pi
au2km = 149597871.0 # km
Msun = 1.99e30 # kg
Mearth = 5.97e24 # kg
GM = 132712440018.0 # km3/s2
G = GM/Msun
s2yr = 1.0 / 3600.0 / 24.0 / 365.25 # yr

Disk Parameters

In [6]:
# Surface Density Profile
Sigma0 = 6.1 # g/cm2
Sigma0 *= 1.0e5 * 1.0e5 # g/km2
Sigma0 *= 1.0 / 1000.0  # kg/km2
Sigma0 *= au2km * au2km # kg/AU2
Sigma0 *= 1.0 / Mearth  # M_Earth/AU2
R0 = 1.0 # AU

# Aspect Ratio
aspect = 0.002

# Particle Number, Particle Mass
Npart = 2048.0
m = (5.0 * Mearth / Npart) # kg

Evaluate Analytical Estimate

In [7]:
R1d, w = np.linspace(0.5, 4.0, 64, retstep=True)
p1d = np.logspace(-4.0, 0.0, 4096)
R, p = np.meshgrid(R1d, p1d)

# Make Disk [Number Density - 1/AU3]
n3d, n2d = compute_number_density(R, w, Sigma0, R0, Npart) 

# Hill Radius
rhill = compute_hill_radius(R, m)

# Relative Velocity, Shear/Dispersion Driven, Depending Which is Larger
vrel = compute_vrel(R*au2km, p*au2km, GM) # km/s
vrel = np.maximum(vrel,0.18) 

# Mean Free Path, Collision Time [Eq. 9]
mfp = 1.0 / (n3d * np.pi * p**2.0) # AU
tau = mfp * au2km / vrel           # s
tau = tau * s2yr                   # yr

# Displacement Factor [Part of Eq. 10]
ddp = (G * m * tau/s2yr) / ((p*au2km)**2.0 * vrel) # [-]

# Exponential Timescale  [Eq. 11]
te = (tau/s2yr) / np.log(1.0 + G * m * n3d/au2km**3.0 * (tau/s2yr)**2.0) # s
te = te * s2yr

# Displacement After 40 Years [Eq. 10]
d40 = (1.0+ddp)**(40.0/tau)

Generate Figure

Displacement Map

In [8]:
fig = plt.figure()
fig.set_size_inches(8.27*0.5,8.27*(6./8.)*0.5)
ax = fig.add_subplot(1,1,1)

# Contour Shading
ct1 = ax.contourf(R, p, d40, \
                  levels=[1,1e1,1e2,1e3,1e4,1e5,1e6], \
                  cmap='gray_r', \
                  locator=mpl.ticker.LogLocator())

# Contour Lines
ctx = ax.contour(R, p, d40, \
                 levels=[1,1e1,1e2,1e3,1e4,1e5,1e6], \
                 linewidths=1.0, colors='k')

# Hill Radii
ax.plot(R1d, 10.0 * compute_hill_radius(R1d, m), \
        'k--', lw=2, label='10 Hill Radii')
ax.plot(R1d, 3.0 * compute_hill_radius(R1d, m), \
        'k--', lw=2, label='3 Hill Radii')
ax.plot(R1d, 1.0 * compute_hill_radius(R1d, m), \
        'k--', lw=2, label='1 Hill Radii')

# Hill Radius Markers
ax.text(0.9, 0.725, r"10 R$_\mathrm{Hill}$", \
         horizontalalignment='right', color='black', \
         transform=ax.transAxes, size='x-small')

ax.text(0.9, 0.60, r"3 R$_\mathrm{Hill}$", \
         horizontalalignment='right', color='black', \
         transform=ax.transAxes, size='x-small')

ax.text(0.9, 0.36, r"1 R$_\mathrm{Hill}$", \
         horizontalalignment='right', color='black', \
         transform=ax.transAxes, size='x-small')

# Colors
fig.colorbar(ct1)

# Scaling, Labels
ax.set_yscale('log')
ax.set_xlabel('Semi-Major Axis (AU)')
ax.set_ylabel('Impact Parameter (AU)')

# Title
ax.set_title("Separation Ratio @ 40 Years")

pass

Separation ratios grow to $10^4$ within 40 years in the most extreme regions. This conincides with the empirical evidence to within an order of magnitude. Not bad for a toy model that only consider the effect of one particle at a time.

Bonus: E-Folding Time

In [9]:
fig = plt.figure()
fig.set_size_inches(8.27*0.5,8.27*(6./8.)*0.5)
ax = fig.add_subplot(1,1,1)

# Contour Shading
ct1 = ax.contourf(R, p, te, \
                  levels=[1,1e1,1e2,1e3,1e4,1e5,1e6], \
                  cmap='gray_r', \
                  locator=mpl.ticker.LogLocator())

# Contour Lines
ctx = ax.contour(R, p, te, \
                 levels=[1,1e1,1e2,1e3,1e4,1e5,1e6], \
                 linewidths=1.0, colors='k')

# Hill Radii
ax.plot(R1d, 10.0 * compute_hill_radius(R1d, m), \
        'k--', lw=2, label='10 Hill Radii')
ax.plot(R1d, 3.0 * compute_hill_radius(R1d, m), \
        'k--', lw=2, label='3 Hill Radii')
ax.plot(R1d, 1.0 * compute_hill_radius(R1d, m), \
        'k--', lw=2, label='1 Hill Radii')

# Hill Radius Markers
ax.text(0.9, 0.725, r"10 R$_\mathrm{Hill}$", \
         horizontalalignment='right', color='black', \
         transform=ax.transAxes, size='x-small')

ax.text(0.9, 0.60, r"3 R$_\mathrm{Hill}$", \
         horizontalalignment='right', color='black', \
         transform=ax.transAxes, size='x-small')

ax.text(0.9, 0.36, r"1 R$_\mathrm{Hill}$", \
         horizontalalignment='right', color='black', \
         transform=ax.transAxes, size='x-small')

# Colors
fig.colorbar(ct1)

# Scaling, Labels
ax.set_yscale('log')
ax.set_xlabel('Semi-Major Axis (AU)')
ax.set_ylabel('Impact Parameter (AU)')

# Title
ax.set_title("E-Folding Time (Years)")

pass

E-folding times for encounters with impact parameters ~ Hill radius are on the order of 1-100 years. Compared the empirical tests, this is an overestimate, but remains within an order of magnitude.

In [ ]: