04. Phase curve fit analysis (Intermediate)

Contact authors: Yumi Choi and Melissa Graham

Last verified to run: 2023-12-22

Targeted learning level: Intermediate

Credits: This tutorial incorporates material from the DP0.3 tutorial notebook on the introduction to phase curves by Christina Williams and Yumi Choi.

Introduction

This portal tutorial is the same demonstration used in the tutorial notebook Introduction to Phase Curves to illustrate the phase curves of solar system objects, but focuses on Main Belt Asteroids (MBAs).

Phase curve fits and absolute magnitudes

Solar system objects in the DP0.3 catalogs change position and brightness between each Rubin image as they orbit about the Sun over time. In the DP0.3 catalogs, their intrinsic properties and orbital parameters are known, and are used to estimate what the measurements would be in a given image, and how they change between images. From these simulated observations, it is possible to reconstruct their intrinsic properties and orbital parameters in the same way as will be done using the real LSST data.

An important way to characterize intrinsic properties of a solar system object is by measuring its “phase curve”, which is the object brightness as a function of its “solar phase angle” (the angle made between the line of sight from the object to the Sun, and the line of sight from the object to Earth). For Solar System objects, absolute magnitudes (H) are defined to be for an object 1 au from the Sun and 1 au from the observer, and at a phase angle \(\alpha\) (the angle Sun-object-Earth) of 0 degrees. Absolute magnitudes are derived by fitting a function to the relationship between reduced magnitude \(H(\alpha)\) and phase angle \(\alpha\) (i.e., the phase curve), and evaluating the function at a phase angle of 0 degrees. The reduced magnitude, which is corrected for distance, \(H(\alpha)\), as:

\[H(\alpha) = m - 5 \log_{10}(d_{t} d_{h}),\]

where \(\alpha\) is the phase angle, \(d_{t}\) is the topocentric distance, \(d_{h}\) is the heliocentric distance, and m is the apparent magnitude.

The absolute magnitude H can be derived by fitting a function, where the choice of form for this function has several options (see Muinonen et al. 2010). The simple two-parameter model (HG model) for asteroids provides a best fit for the absolute magnitude H and the slope parameter G. The HG model parameters were used to predict the observed parameters for each object. These truth values are defined in the MPCORB table as mpcH (intrinsic absolute magnitude in V band) and mpcG (intrinsic slope, which has a constant value of 0.15 for all DP0.3 objects). The HG model has the form:

\[H = H(\alpha) + 2.5 \log_{10}((1-G)\phi_1(\alpha) + G\phi_2(\alpha)).\]

To better accommodate various observational effects (e.g., photometric quality, incomplete phase angle sampling) a more sophisticated HG1G2 model (a linear three-parameter function) and its nonlinear two-parameter version HG12 model were developed (Muinonen et al. 2010). The two-parameter HG12 model is generally very effective for deriving reliable values of absolute magnitude when the phase angle sampling is not optimal (e.g., poor phase angle coverage at a range of phase angle). Thus, the LSST data products will compute estimated parameters of the HG12 model and this will be the focus of this tutorial. The HG12 model expresses the G1 and G2 parameters as a piecewise linear function of a single parameter, G12:

\[H(\alpha) = H − 2.5 \log_{10}[G1\phi_1(\alpha)+G2\phi_2(\alpha) + (1-G1-G2)\phi_3(\alpha)],\]

where:

\(G1 = 0.9529 \times G12 + 0.02162\), \(G2 = -0.6125 \times G12 + 0.5572\) for \(G12 \ge 0.2\), and \(G1 = 0.7527 \times G12 + 0.06164\), \(G2 = -0.9612 \times G12 + 0.6270\) for \(G12 < 0.2\).

The results of phase curve fits (H and G12) in the four LSST filters, griz, are stored in the SSObject table. The explanation for the absence of u and y bands in DP0.3 catalogs can be found in the DP0.3 documentation. Note that rotation curves or complex geometry of solar system objects are not included in DP0.3 simulations. Thus, any changes over time in an object’s apparent magnitude are due only to changes in its distance and phase angle.

Step 1. Query the DP0.3 tables for the Main Belt Asteroids

1.1. Log into the Rubin Science Platform at data.lsst.cloud and select the Portal Aspect. At upper right, next to “TAP Services” choose to “Show”, and then select “LSST DP0.3 SSO” from the drop-down menu that appears at the top.

1.2. At upper right, next to “View” choose “Edit ADQL”. Enter the query statement below into the ADQL Query box and execute the query to select a good number of MBAs with a fair number of total observations (numObs > 100) to explore the distribution of their properties. Following the population definitions used by the JPL Horizons small body database query tool, we select MBAs as objects with semi-major axes a between 2.0 au and 3.25 au and perihelion distance q > 1.666 au. Note that semi-major axes are not directly available in the MPCORB table, so the constraint on a is derived from perihelion q and eccentricity e. You might want to increase the “Row limit” to 100,000 to have an appreciable sample of objects by entering this number in the box on the lower left. In order to have the query execution not to take too long, we restrict the number of returned objects to have their mpc.ssObjectId in the limited range.

SELECT mpc.ssObjectId, mpc.e, mpc.q, mpc.mpcG, mpcH,
sso.numObs,
sso.g_H, sso.g_Herr, sso.g_G12, sso.g_G12err, sso.g_Ndata,
sso.r_H, sso.r_Herr, sso.r_G12, sso.r_G12err, sso.r_Ndata,
sso.i_H, sso.i_Herr, sso.i_G12, sso.i_G12err, sso.i_Ndata,
sso.z_H, sso.z_Herr, sso.z_G12, sso.z_G12err, sso.z_Ndata
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 (mpc.q / (1-mpc.e)) > 2.0
AND (mpc.q / (1-mpc.e)) < 3.25
AND (mpc.q > 1.666)
AND sso.numObs > 100

1.3. Plot the distribution of semi-major axes a of orbits of the objects in your query. The execution of the query will result in a blank panel for the plot, with a comment “Cannot display the requested data.” To plot the distribution of a you need to click on the “Chart options and tools” icon (two gears), click on “Add New Chart”, select “Histogram” for “Plot Type”, enter “q / (1-e)” as the “Column or expression” and “100” for number of bins as on the screenshot below.

Figure 1: The “Plot Parameters” pop-up window to set parameters for making a histogram of semi-major axes for MBAs.

1.4. Click “Ok” and close the chart stating “cannot display requested data” by clicking the blue “X” mark in its upper right hand corner. It will result in the following plot and table below. Note that the distribution of asteroids as a function of semi-major axis is not uniform, but it reveals a number of peaks and gaps where there are very few (or no) objects. These are known as Kirkwood gaps, which arise due to resonances between the asteroid’s and Jupiter’s orbital periods.

Figure 2: The distribution of semi-major axes for MBAs. The prominent Kirkwood gaps in this plot are located at 2.065 au (4:1 resonance), 2.502 au (3:1 resonance), 2.825 au (5:2 resonance), and 2.958 au (7:3 resonance).

Step 2. Select a well-observed MBA, and plot its phase curve

2.1. Unique solar system objects in the SSObject and MPCORB tables will be observed many times over the full LSST survey. Individual observations of each unique object in each filter are recorded in the SSSource and diaSource tables. Below, we query these tables to obtain all of the individual observed time series data (we call indivObs) for an MBA that has more than 2000 observations. First, select MBAs with 2000 or more observations by entering “>2000” in the box underneath the table heading numObs as shown as below and hitting the return key. This will leave only a small fraction of queried 100,000 MBAs above, 23 MBAs in this tutorial. To go back to the originally retreived table by removing the applied filter, click the remove filter icon, which is the first icon on the top right of the table.

Figure 3: The resulting table of 23 MBAs with 2000 or more observations out of the retrieved 100,000 MBAs in Step 1.2.

2.2. Pick and copy one ssObjectId. Hovering over a table cell shows you a triple-dot box. Right-click on that box, two options will pop up: “Copy to clipboard” and “View as plain text”. Here, copy ssObjectId = 7470575696289418102 to clipboard and click “RSP TAP Search” button on the top left to go back to the ADQL Query page.

Figure 4: How to copy a selected ssObjectId to clipboard.

2.3 Execute the following ADQL query to retrieve the apparent magnitudes, magnitude errors, filters, phase angles, topocentric and heliocentric distances of the individual observations for a well-observed MBA.

SELECT
dia.ssObjectId, dia.mag, dia.magErr, dia.band,
sss.phaseAngle, sss.topocentricDist, sss.heliocentricDist
FROM dp03_catalogs_10yr.DiaSource as dia
INNER JOIN dp03_catalogs_10yr.SSSource as sss ON dia.diaSourceId = sss.diaSourceId
WHERE dia.ssObjectId = 7470575696289418102

2.4. To plot the phase curve in the g-band (i.e, reduced magnitude versus phase angle), first select the g-band data by entering “=’g’” in the box underneath the table heading band and hitting the return key. Then open the “Plot Parameters” pop-up window (click on the two-gear icon), click on “Modify Trace”, set the “X” to phaseAngle and “Y” to mag - 5 * log10(topocentricDist * heliocentricDist). Check the “Error” box for the y-axis and select “Symm”, and put magErr. Click on the “Chart Options” arrow, and set the “X Label” to be “Phase angle [deg]” and the “Y Label” to be “Reduced magnitude”. Check the “reverse” box for the y-axis option.

Figure 5: The “Plot Parameters” pop-up window to plot the phase curve in g-band.

2.5. Click on the “Apply” button. This will result in the g-band phase curve plot with error bars for the MBA with ssObjectId = 7470575696289418102 as shown below.

Figure 6: The g-band phase curve for the MBA with ssObjectId = 7470575696289418102.

2.6. In order to plot a phase curve in a different band, for example in z-band, enter “=’z’” in the box underneath the table heading band and hitting the return key. The g-band phase curve plot will be replaced with the z-band phase curve plot as shown below. It is clear that the phase curves of the source are offset from each other in these two filters, reflecting the variation in brightness of asteroids in different filters. Also the reduced magnitude qualities (i.e., photometric uncertainties) are significantly different.

Figure 7: Same as Figure 6, but in z-band.

Step 3. Exploring phase curve data products from the DP0.3 catalogs

3.1. This section explores the distribution of G12 slope parameter values as a function of absolute magnitudes H for MBAs in griz bands. Return to the originally retrieved table in Step 1.2 by clicking the first table tab. Remove the numObs > 2000 condition either by clicking the remove filter icon on the top right or by deleting the condition and hitting the return key. Then, open the “Plot Parameters” pop-up window (click on the two-gear icon), choose “Add New Chart”, opt for “Heatmap” as the “Plot Type”, and create a new plot for the G12 vs. H in g band, adhering to the specified plot settings below.

Figure 8: The “Plot Parameters” pop-up window to plot the G12 vs. H in g band.

3.2. Once creating the G12 vs. H plot for g-band, close the histogram of semi-major axes of MBAs we made in Step 1.4, and add three more new plots for riz bands by repeating the creation of the G12 vs. H plot in Step 3.1, but going through the riz bands. This will generate four panels as shown below.

Figure 9: The G12 vs. H plots in griz from top left to bottom right clockwise.

3.3. Recall that the input (truth) G value using the HG_model that was used to generate the DP0.3 simulated object’s observed properties was fixed across the population to a constant value of G = 0.15 (refer The DP0.3 Simulation). The DP0.3 automated phase curve fitting (which uses HG12_model) produces a nearly constant value for G12 with a relatively small spread at bright magnitudes, and the scatter in measured G12 starts to deviate more substantially at fainter magnitudes where it is likely harder to recover the intrinsic value.

3.4. This section explores the impact of the total number of observations for a given source (numObs) and the perihelion distance (q) on the quality of phase curve fitting in i-band as an example. First close any open plots except for one heatmap, and then click on “Chart options and tools” icon to make a new plot. Select “Modify Trace”, set the “X” to numObs, “Y” to i_Herr, the number of “X”- and “Y”-bins to 200. Lastly, set the min and max for the y-axis under the “Chart Options” to be 0 and 0.05 as follows. Make another plot by repeating the same paramter setting, but selecting “Add New Chart” and entering q on the x-axis.

Figure 10: The “Plot Parameters” pop-up window to plot the i_Herr vs. numObs.

3.5. Make two new plots by repeating the above, but setting the “Y” to i_G12err. This will generate four panels showing how the H and G12 parameter uncertainties vary with the total number of observations and the perihelion distance for MBAs. In left panels, it is clear that the phase curve fit uncertainties decrease with the number of observations of each source. So as LSST accumulates data over time, precision in the phase curve modeling will improve. The right panels show that uncertainties in the phase curve parameters modestly increase for objects with larger perihelion distances.

Figure 11: Uncertainties in i_Herr and i_G12err as a function of the total number of observations and the perihelion distance.

3.6. The above plots compare numObs (total in all bands) with model fits, which may not be the ideal metric since the quality of phase curves can vary quite a bit between filters. Instead, one can look at the number of datapoints included in the phase curve modeling on a per filter basis (e.g., r_Ndata for the r-band in the SSObject table). To make a plot showing the distribution of the number of observations in each filter, again first close any open plots except for one, and then click on the “Chart options and tools” icon. Select “Add New Chart”, set the “Plot Type” to “Histogram”, the “Column or expression” to g_Ndata. Select the “Uniform binning” algorithm, set the number of bins to 100 with the min and max to be 0 and 1300, respectively. Under the “Chart Options”, check the “log” box for the y-axis. It will plot the histogram of the g-band number of observations. Once creating the g_Ndata histogram, close the remaining plot from Step 3.5. To overplot the histogram for r_Ndata, select “Overplot New Trace” on the “Plot Parameters” pop-up window, and use the same plot parameters, but change the “Column or expression” to r_Ndata. Now the “Name” box under the “Trace Options” will appear, where you can set legend for each histogram. Once overplotting the r_Ndata histogram, the “Choose Trace” field with a drop-down menu becomes available when you reopen the “Plot Parameters” pop-up window. Choose “Trace 0” and enter a label for the g_Ndata histogram in the “Name” box under the “Trace Options”. Repeat this process for the i and z bands as well. For the z band plot, set the “X Label” to “Number of observations”. Note that r- and i-bands produce the most data points for recovering phase curves, while g- and z-band produce much less. Phase curves measured in r- and i-bands will thus be better sampled. Clicking the labels in the legend makes it possible to show and hide each histogram.

Figure 12: Histograms of the number of observations in each filter.

3.7. To confirm whether phase curve fitting in i band is indeed more precise than in z band, let’s compare the uncertainty in H values for i and z bands by adding a new plot. Click on the two-gear icon, select “Add New Chart” and “Heatamp”. Set the “X” to i_Herr and “Y” to z_Herr with the X and Y MIN/MAX of 0 and 0.1. To make the plot with a more proper display ratio, set the “X/Y ratio” under the “Chart Options” to be 1, check the “width” box, and hit the apply button. The right panel in the figure below, one can see that poorer sampling drives higher uncertainty in the derived absolute magnitude H using z band compared to i band for MBAs.

Figure 13: Comparison of the uncertainty in the measured H values in i and z bands.

Step 4. Exercises for the learner

4.1 Explore phase curves for objects with less phase angle coverage and compare them with those for MBAs. For example, Trans-Neptunian Objecs (refer the portal tutorial 03. Explore Trans-Neptunian Objects (TNOs) in DP0.3 (Intermediate)) or Jupiter Trojans (refer the notebook tutorial DP03_04a_Introduction_to_Phase_Curves).