Basic usage
Computing solid surface densities
You can compute the solid surface densities of the annular disks required to form each planet from their masses and semi-major axes:
sigmas = solid_surface_density_CL2013(planet_masses, planet_semi_a) # g/cm^2
where planet_masses
is an array of the planet masses (in Earth masses) and planet_semi_a
is an array of the planet semi-major axes (in AU).
Here, we’ve used the Chiang & Laughlin (2013) prescription for the feeding zone of each planet, which simply sets them equal to their semi-major axes; the output sigmas
is an array of solid surface densities computed for each planet.
Tip
Several other prescriptions are available; see solid_surface_density_S2014
, solid_surface_density_nHill
, and solid_surface_density_system_RC2014
.
Some of these may require additional parameters, such as the planet radii (for “S2014”) and the stellar mass (for “S2014” and “nHill”).
You can also define and use your own custom prescription for the feeding zone width (delta_a
) by using the core function that all of the above functions call, solid_surface_density
.
Fitting a power-law
It is easy to fit a power-law to the solid surface densities as a function of semi-major axis:
a0 = 0.3 # AU; separation for the normalization constant
sigma0, beta = fit_power_law_MMEN(planet_semi_a, sigmas, a0=a0)
Here, sigma0
is the normalization constant for the solid surface density (g/cm^2) at the chosen separation a0
, and beta
is the slope or power-law index.
An example using a SysSim catalog
A full working example using a simulated catalog is shown below (make sure to replace load_dir
with the path to where you downloaded the simulated catalog):
from syssimpyplots.general import *
from syssimpyplots.load_sims import *
from syssimpymmen.mmen import *
# Load a simulated catalog:
load_dir = '/path/to/a/simulated/catalog/' # replace with your path where you downloaded a simulated catalog!
sss_per_sys, sss = compute_summary_stats_from_cat_obs(file_name_path=load_dir)
# Compute solid surface densities for all planets in the catalog:
sigma_obs_CL2013, mass_obs, a_obs = solid_surface_density_CL2013_given_observed_catalog(sss_per_sys, max_core_mass=np.inf)
# Fit a power-law profile for the MMEN:
a0 = 0.3 # AU; choice for the normalization point
sigma0_obs_CL2013, beta_obs_CL2013 = fit_power_law_MMEN(a_obs, sigma_obs_CL2013, a0=a0)
# Plot solid surface density vs semi-major axis:
fig = plt.figure(figsize=(8,8))
plot = GridSpec(1,1,left=0.15,bottom=0.15,right=0.95,top=0.95,wspace=0,hspace=0)
ax = plt.subplot(plot[0,0])
plt.scatter(a_obs, np.log10(sigma_obs_CL2013), marker='o', s=10, alpha=0.2, color='b')
a_array = np.linspace(1e-3,2,1001)
plt.plot(a_array, np.log10(MMEN_power_law(a_array, sigma0_obs_CL2013, beta_obs_CL2013, a0=a0)), lw=2, ls='--', color='k', label=r'Power-law fit ($\Sigma_0 = {:0.0f}$, $\beta = {:0.2f}$)'.format(sigma0_obs_CL2013, beta_obs_CL2013))
ax.tick_params(axis='both', labelsize=20)
plt.gca().set_xscale("log")
ax.get_xaxis().set_major_formatter(ticker.ScalarFormatter())
plt.xticks([0.05, 0.1, 0.2, 0.4, 0.8])
plt.xlim([0.04,0.9])
plt.ylim([0.5,5.5])
plt.xlabel(r'Semimajor axis, $a$ (AU)', fontsize=20)
plt.ylabel(r'Surface density, $\log_{10}(\Sigma/{\rm g cm}^{-2})$', fontsize=20)
plt.legend(loc='upper right', bbox_to_anchor=(1.,1.), ncol=1, frameon=False, fontsize=16)
plt.show()

The solid surface densities and MMEN fit to a simulated observed catalog.
Note
In this example, the function solid_surface_density_CL2013_given_observed_catalog
takes the observed catalog, uses a mass-radius relation to draw a set of planet masses from the planet radii, and enforces a limit on the maximum mass via max_core_mass
before calling the function solid_surface_density_CL2013
we started with at the top of this page.
We’ve removed the maximum core mass limit (by setting max_core_mass=np.inf
) to show you the broad range of surface densities arising from the broad range of planet masses. By default, it is set to 10 Earth masses.