02. Main Belt Asteroids#
Properties of Main Belt Asteroids in DP0.3
Contact author(s): Jeff Carlin
Last verified to run: 2024-05-01
LSST Science Pipelines version: Weekly 2024_16
Container Size: medium
Targeted learning level: beginner
Description: An introduction to Main Belt asteroids in the DP0.3 dataset, and illustration of some simple analyses of their orbital families.
Skills: Use TAP queries to retrieve Solar System objects. Plot properties and orbits of Main Belt asteroids.
LSST Data Products: TAP DP0.3 MPCORB, SSObject, and SSSource tables.
Packages: lsst.rsp.get_tap_service
Credit: Originally developed by Jeff Carlin and the Rubin Community Science Team in the context of the Rubin DP0.
Get Support: Find DP0-related documentation and resources at dp0.lsst.io. Questions are welcome as new topics in the Support - Data Preview 0 Category of the Rubin Community Forum. Rubin staff will respond to all questions posted there.
1. Introduction¶
This tutorial demonstrates how to use TAP to query the MPCORB and SSObject tables. The tables are queried simultaneously and joined on their SSObjectId, which allows one to retrieve the measured magnitudes from SSObject and orbital properties from MPCORB. We examine the properties of a large sample of Solar System objects from the DP0.3 catalogs, then focus on asteroids in the Main Belt (between about 1.6-4.2 au). We examine the orbital properties of Main Belt asteroids, and plot orbits of some of them to look for differences between the orbital families.
A detailed introduction to the DP0.3 data products can be found in Notebook DP03_01_Introduction_to_DP03.
1.1 Package Imports¶
The matplotlib
and numpy
libraries are widely used Python libraries for plotting and scientific computing. We will use these packages below, including the matplotlib.pyplot
plotting sublibrary.
We also use the lsst.rsp
package to access the TAP service and query the DP0 catalogs.
# general python packages
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('tableau-colorblind10')
# LSST package for TAP queries
from lsst.rsp import get_tap_service
service = get_tap_service("ssotap")
1.2 Define Functions and Parameters¶
Define a function to convert a given perihelion distance (q
) and eccentricity (e
) to an orbital semi-major axis (a
). Their relationship is defined by q = a(1-e).
def calc_semimajor_axis(q, e):
"""
Given a perihelion distance and orbital eccentricity,
calculate the semi-major axis of the orbit.
Parameters
----------
q: float
Distance at perihelion, in au.
e: float
Orbital eccentricity.
Returns
-------
a: float
Semi-major axis of the orbit, in au.
q = a(1-e), so a = q/(1-e)
"""
return q / (1.0 - e)
Define a function to calculate the semi-major axis (a) at which a mean-motion resonance with a planet occurs via Kepler's Third Law: (P/P_planet)^2 = (a/a_planet)^3, where a
is the semi-major axis of a small body in au, a_planet
is the planet's semi-major axis in au, P
is the small body's orbital period in yr, and P_planet
is the planet's orbital period in yr.
def calc_mean_motion_res_a(P_over_P_planet, a_planet):
"""
Implements Kepler's Third law,
`(P/P_planet)^2 = (a/a_planet)^3`, where
`a` is the semi-major axis in au,
`a_planet` is the planet's semi-major axis in au,
`P` is the small body's orbital period in yr, and
`P_planet` is the planet's orbital period in yr,
to find the semimajor axis at which a mean-motion
resonance with a planet occurs, via
`a = a_planet * (P/P_planet)^(2/3)`.
Parameters
----------
P_over_P_planet: ndarray
Integer ratio of a small body's orbital period
to a planet's orbital period.
a_planet: ndarray
Semimajor axis of a planet, in au.
Returns
-------
a_res: ndarray
Semimajor axis of a mean motion resonance, in au.
a_res = (P/P_planet)^(2/3) * a_planet
"""
return (float(P_over_P_planet))**0.66667 * a_planet
Define a function to take an ssObjectId
as an input, and return all measurements of that object from the SSSource
table.
def get_obj_matched(objid):
"""
Given an objectId, retrieve all measurements from
SSSource for that object.
Parameters
----------
objid: int64
ObjectId of the object of interest.
Returns
-------
allsrc: dataframe
Measurements for the object from the DP0.3 SSSource table.
"""
allsrc = service.search("SELECT * FROM dp03_catalogs_10yr.DiaSource as diasrc \
JOIN dp03_catalogs_10yr.SSSource as sssrc \
ON diasrc.diaSourceId = sssrc.diaSourceId \
WHERE diasrc.ssObjectId = \
" + str(objid)).to_table().to_pandas()
return allsrc
Define a function to plot orbits of multiple Solar System objects on the same figure. The orbits are plotted in XYZ Cartesian coordinates, either in a heliocentric (Sun-centered) or topocentric (relative to the Earth's surface) frame.
def xyz_orbit_plot_manyobjects(tobs_list, x_list, y_list, z_list, title=None):
"""
Given lists (of lists) of observation times and XYZ coordinates,
plot the orbits of multiple SS objects.
Parameters
----------
tobs_list: list of lists
MidpointTAI times of each observation of each object.
{xyz}_list: list of lists
Solar system positions of each observation of each object.
Can be either heliocentric or topocentric.
title: string
(Optional) title to add to the top of the plot.
"""
params = {
'axes.labelsize': 12,
'font.size': 12,
'legend.fontsize': 8,
'xtick.major.width': 2,
'xtick.minor.width': 1,
'xtick.major.size': 6,
'xtick.minor.size': 4,
'xtick.direction': 'in',
'xtick.top': True,
'lines.linewidth': 2,
'axes.linewidth': 2,
'axes.labelweight': 2,
'axes.titleweight': 3,
'ytick.major.width': 2,
'ytick.minor.width': 1,
'ytick.major.size': 6,
'ytick.minor.size': 4,
'ytick.direction': 'in',
'ytick.right': True,
'text.usetex': False,
'figure.figsize': [6, 6],
'figure.facecolor': 'white'
}
plt.rcParams.update(params)
fig = plt.figure()
fig.subplots_adjust(hspace=0, wspace=0)
for i in range(len(tobs_list)):
col = (np.random.random(), np.random.random(), np.random.random())
ax1 = plt.subplot(221)
ax1.plot(x_list[i], z_list[i], '.', color=col, ms=5)
ax1.set_ylabel('Z (au)')
ax1.set_xlabel('X (au)')
if title:
ax1.set_title(title)
ax1.minorticks_on()
ax2 = plt.subplot(222, sharey=ax1)
ax2.plot(y_list[i], z_list[i], '.', color=col, ms=5)
ax2.set_xlabel('Y (au)')
ax2.minorticks_on()
plt.setp(ax2.get_yticklabels(), visible=False)
ax3 = plt.subplot(223, sharex=ax1)
ax3.plot(x_list[i], y_list[i], '.', color=col, ms=5)
ax3.set_ylabel('Y (au)')
ax3.set_xlabel('X (au)')
ax3.minorticks_on()
plt.show()
Set some plotting defaults to make plots look nice.
plt.style.use('tableau-colorblind10')
%matplotlib inline
params = {'axes.labelsize': 20,
'font.size': 20,
'legend.fontsize': 14,
'xtick.major.width': 3,
'xtick.minor.width': 2,
'xtick.major.size': 8,
'xtick.minor.size': 4,
'xtick.direction': 'in',
'xtick.top': True,
'lines.linewidth': 3,
'axes.linewidth': 3,
'axes.labelweight': 3,
'axes.titleweight': 3,
'ytick.major.width': 3,
'ytick.minor.width': 2,
'ytick.major.size': 8,
'ytick.minor.size': 4,
'ytick.direction': 'in',
'ytick.right': True,
'figure.figsize': [6, 6],
'figure.facecolor': 'White'
}
plt.rcParams.update(params)
2. Orbital parameters of Solar System objects¶
Query the DP0.3 catalogs, joining the MPCORB and SSObject tables on their ssObjectId. One could retrieve the entire catalog (containing millions of rows), but that query may take a long time. To minimize query time but still extract a large sample, we use a range of ssObjectId
values.
df_mpc_sso = service.search("SELECT mpc.ssObjectId, mpc.e, mpc.incl, mpc.q, mpc.peri, \
sso.ssObjectId, sso.g_H, sso.r_H, sso.i_H, sso.z_H \
FROM dp03_catalogs_10yr.MPCORB as mpc \
JOIN dp03_catalogs_10yr.SSObject as sso \
ON mpc.ssObjectId = sso.ssObjectId \
WHERE mpc.ssObjectId < 9223370875126069107 \
AND mpc.ssObjectId > 7331137166374808576 \
AND sso.numObs > 50").to_table().to_pandas()
print('The query returned ', len(df_mpc_sso), ' results.')
The query returned 339864 results.
2.1 Overall orbital properties¶
First, calculate the semi-major axis of each object's orbit, using the function we defined above.
a = calc_semimajor_axis(df_mpc_sso.q, df_mpc_sso.e)
Plot the orbital inclination and the eccentricity vs. semi-major axis. Note that we plot the sine of the inclination so its range is limited to values between 0 and 1.
fig, axs = plt.subplots(2, 1, figsize=(7, 5), sharex=True)
fig.subplots_adjust(hspace=0)
axs[0].hexbin(a, np.sin(np.deg2rad(df_mpc_sso.incl)),
gridsize=(int(103/0.2), 50), cmap='Spectral_r', bins='log',
extent=(0, 103, 0, 1))
axs[0].set_ylabel('sin(incl)')
axs[0].minorticks_on()
axs[0].set_xlim(0, 103)
axs[1].hexbin(a, df_mpc_sso.e,
gridsize=(int(103/0.2), 50), cmap='Spectral_r', bins='log',
extent=(0, 103, 0, 1))
axs[1].set_ylabel('eccentricity')
axs[1].set_xlabel('semi-major axis (au)')
axs[1].set_xlim(0, 103)
axs[1].set_ylim(0, 1.09)
axs[1].minorticks_on()
plt.show()
plt.close()
This plot shows Solar System objects extending out to at least 100 au, but the majority of objects are in the inner Solar System (at less than 6 au). For the rest of this notebook, we will focus on the inner Solar System -- specifically on the so-called "Main Belt" asteroids. (The outer Solar System bodies will be the subject of another notebook.)
2.2 Main Belt asteroids¶
2.2.1 Semi-major axes of Main Belt asteroids¶
This plot shows Solar System objects extending out to at least 100 au, but the majority of objects shown are at less than 10 au. For the rest of this notebook, we will focus on the so-called "Main Belt" asteroids, which roughly speaking inhabit the space between the orbits of Mars and Jupiter. The exact Main Belt definition is not uniform in the literature; for DP0.3 we adopt the definition used by the JPL small-body database, where the definitions of the orbital families are visible as the popup notes associated with the "orbital class" info. This definition places the Main Belt between 1.6 au < a < 4.6 au, with the constraint q > 1.6 au to ensure that objects don't cross Mars' orbit. However, in this notebook, we alter the definition a bit and consider objects with 1.6 au < a < 4.2 au, q > 1.6, to allow comparison to figures from other papers.
Before we get into the orbital elements of the Main Belt population, let's just see where they are located in the Solar System (via their semi-major axis (a)
):
main_belt = (a > 1.6) & (a < 4.2)
fig = plt.figure(figsize=(6, 4))
plt.hist(a[main_belt], bins=np.arange(1.6, 4.3, 0.01), color='Black')
plt.xlabel('semi-major axis (au)')
plt.ylabel('number of objects')
plt.minorticks_on()
plt.show()
Clearly it's not just a uniform "belt" between Mars and Jupiter -- it has some "gaps" where there are very few (or no) objects. These are known as Kirkwood gaps, which arise due to resonances between the asteroid's orbital period and that of Jupiter. For example, at a=2.5 au, an asteroid orbits 3 times for each single orbit of Jupiter, and is thus in a "3:1 resonance". Here's a figure illustrating the positions of some of the resonances, and the names of some of the prominent families (source).
You can see that the Hildas, Cybeles, and Hungarias are clearly visible in the DP0.3 dataset.
We can also add the locations of the mean-motion resonances present in the main belt to our histogram above. We use Kepler's Third Law, (P/P_planet)^2 = (a/a_planet)^3, with the ratio of the orbital periods of the asteroid (P
) and planet (P_planet
) along with the semi-major axis of the planet (a_planet
) to caclulate the semi-major axis at which the corresponding mean-motion resonance occurs: a_res = (P/P_planet)^(2/3) * a_planet. For example, the 3:1 mean-motion resonance with Jupiter, which corresponds to an asteroid orbital period that is 1/3 as long as Jupiter's orbital period, occurs at a_res = ((1/3)^0.66667 * (5.2 au) = 2.5 au.
a_res_J_3_1 = calc_mean_motion_res_a((1/3), 5.2)
a_res_J_5_2 = calc_mean_motion_res_a((2/5), 5.2)
a_res_J_7_3 = calc_mean_motion_res_a((3/7), 5.2)
a_res_J_2_1 = calc_mean_motion_res_a((1/2), 5.2)
a_res_J_4_1 = calc_mean_motion_res_a((1/4), 5.2)
a_res_J_3_2 = calc_mean_motion_res_a((2/3), 5.2)
fig = plt.figure(figsize=(6, 4))
plt.hist(a[main_belt], bins=np.arange(1.6, 4.3, 0.01), color='Black')
plt.axvline(x=a_res_J_3_1, color='black', linestyle='dashed', linewidth=0.5)
plt.axvline(x=a_res_J_5_2, color='black', linestyle='dashed', linewidth=0.5)
plt.axvline(x=a_res_J_7_3, color='black', linestyle='dashed', linewidth=0.5)
plt.axvline(x=a_res_J_2_1, color='black', linestyle='dashed', linewidth=0.5)
plt.axvline(x=a_res_J_4_1, color='black', linestyle='dashed', linewidth=0.5)
plt.axvline(x=a_res_J_3_2, color='black', linestyle='dashed', linewidth=0.5)
plt.annotate('J 3:1', xy=(a_res_J_3_1-0.08, 5000.), xytext=(a_res_J_3_1-0.08, 4800.), rotation=90.0, fontsize=10)
plt.annotate('J 5:2', xy=(a_res_J_5_2-0.08, 5000.), xytext=(a_res_J_5_2-0.08, 4800.), rotation=90.0, fontsize=10)
plt.annotate('J 7:3', xy=(a_res_J_7_3-0.08, 5000.), xytext=(a_res_J_7_3-0.08, 4800.), rotation=90.0, fontsize=10)
plt.annotate('J 2:1', xy=(a_res_J_2_1-0.08, 5000.), xytext=(a_res_J_2_1-0.08, 4800.), rotation=90.0, fontsize=10)
plt.annotate('J 4:1', xy=(a_res_J_4_1-0.08, 5000.), xytext=(a_res_J_4_1-0.08, 4800.), rotation=90.0, fontsize=10)
plt.annotate('J 3:2', xy=(a_res_J_3_2-0.08, 5000.), xytext=(a_res_J_3_2-0.08, 4800.), rotation=90.0, fontsize=10)
plt.xlabel('semi-major axis (au)')
plt.ylabel('number of objects')
plt.minorticks_on()
plt.show()
Now we look at the main-belt asteroid's closest approach to the Sun in their orbits, known as the perihelion distance (q)
:
fig = plt.figure(figsize=(6, 4))
plt.hist(df_mpc_sso[main_belt].q, bins=np.arange(0.0, 4.3, 0.02),
color='Black', log=True, histtype='step')
plt.xlabel('perihelion distance (au)')
plt.ylabel('number of objects')
plt.minorticks_on()
plt.show()
Hmmm... Should we be concerned that some of those cross the Earth's orbit at 1 au(!)? That's left as an exercise for the learner to explore. (Note that Main Belt asteroids are by definition not going to approach the Earth, as they are typically defined to have perihelion distances beyond the orbit of Mars.)
2.2.2 Orbital properties of Main Belt asteroids¶
The above illustration (from the Wikipedia article on Kirkwood gaps) shows planets as large orange circles, and asteroids as small dots of various colors. The Main Belt is clear as a large ring outside Mars' orbit (and inside Jupiter's), and you can clearly see the Kirkwood gaps in the Main Belt. But we can go even further by classifying objects based on their orbital properties.
The following figure is from the Grav et al. 2011 paper describing the PanSTARRS Synthetic Solar System Model, which (in part) forms the basis of DP0.3's simulated dataset. This figure shows the orbital inclinations and eccentricities of "Main Belt asteroids" (according to that paper's definition of the Main Belt) as a function of their semi-major axis.
Let's make a similar plot for DP0.3 data and see how it compares to the earlier simulation. (First, we'll make the cut to remove objects with perihelia < 1.6 au.)
not_near_mars_peri = df_mpc_sso.q > 1.6
Now plot the properties of these selected objects. Note that because DP0.3 (and, ultimately, LSST) has far more detected asteroids than previous SDSS-based studies, we will need to plot density of points instead of a scatter plot. Thus, in each of the following plots, the color scale represents the number of points, from purple/blue as the fewest to red regions containing the most objects.
fig, axs = plt.subplots(2, 1, figsize=(7, 5), sharex=True)
fig.subplots_adjust(hspace=0)
axs[0].hexbin(a[not_near_mars_peri], np.sin(np.deg2rad(df_mpc_sso[not_near_mars_peri].incl)),
gridsize=(int(4.7/0.01), 200), cmap='Spectral_r', bins='log',
extent=(0, 4.7, 0, 1))
axs[0].set_ylabel('sin(incl)')
axs[0].minorticks_on()
axs[0].set_xlim(0, 103)
axs[1].hexbin(a[not_near_mars_peri], df_mpc_sso[not_near_mars_peri].e,
gridsize=(int(4.7/0.01), 200), cmap='Spectral_r', bins='log',
extent=(0, 4.7, 0, 1))
axs[1].set_ylabel('eccentricity')
axs[1].set_xlabel('semi-major axis (au)')
axs[1].set_xlim(1.2, 4.7)
axs[1].set_ylim(0, 1.09)
axs[1].minorticks_on()
plt.show()
plt.close()
3. Orbit families among Main Belt asteroids¶
Besides the Kirkwood gaps, the distribution of asteroids in the above plot clearly shows concentrations ("blobs") of objects that share similar orbital properties. Let's examine those in more detail.
We will follow the analysis from this Parker et al. 2008 paper, which examined orbit families among Main Belt asteroids observed by SDSS. (See also Ivezic et al. 2002 and Chapter 5 of the LSST Science Book.)
First, we will divide the Main Belt asteroids into the following subsets, as in the above figure:
- inner belt: 2.0 < a < 2.5 au
- mid belt: 2.5 < a < 2.82 au
- outer belt: 2.82 < a < 3.25 au
Then, we will plot sin(incl)
(inclination) vs. eccentricity (e)
to compare to the following figure from the Parker et al. paper (apologies for the low resolution -- the x-axis is "e" and the y-axis "sin(i)"):
inner_belt = (a > 2.0) & (a < 2.5) & not_near_mars_peri
mid_belt = (a >= 2.5) & (a < 2.82) & not_near_mars_peri
outer_belt = (a > 2.82) & (a < 3.25) & not_near_mars_peri
Now plot the properties of these groups. Because there are so many points, we will plot the density of points; purple/blue colors represent low-density regions, while red areas have the highest densities.
fig, axs = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
axs[0].hexbin(df_mpc_sso[inner_belt].e, np.sin(np.deg2rad(df_mpc_sso[inner_belt].incl)),
gridsize=300, bins='log', cmap='Spectral_r')
axs[0].set_xlabel('e')
axs[0].set_ylabel('sin(i)')
axs[0].set_xlim(0, 1)
axs[0].set_ylim(0, 1)
axs[0].minorticks_on()
axs[0].set_title(r'inner belt (2.0 < a < 2.5)', fontsize=16)
axs[1].hexbin(df_mpc_sso[mid_belt].e, np.sin(np.deg2rad(df_mpc_sso[mid_belt].incl)),
gridsize=300, bins='log', cmap='Spectral_r')
axs[1].set_xlabel('e')
axs[1].set_xlim(0, 1)
axs[1].set_ylim(0, 1)
axs[1].minorticks_on()
axs[1].set_title(r'mid belt (2.5 < a < 2.82)', fontsize=16)
axs[2].hexbin(df_mpc_sso[outer_belt].e, np.sin(np.deg2rad(df_mpc_sso[outer_belt].incl)),
gridsize=300, bins='log', cmap='Spectral_r')
axs[2].set_xlabel('e')
axs[2].set_xlim(0, 1)
axs[2].set_ylim(0, 1)
axs[2].minorticks_on()
axs[2].set_title(r'outer belt (2.82 < a < 3.25)', fontsize=16)
plt.show()