The Solar System as an Exosystem: Planet Confusion

Future exoplanet direct-imaging missions seek to both discover and spectrally characterize Earth-like exoplanets (HabEx Study Team 2019). Spectral characterizations of Earth-like exoplanets with coronagraphs or starshades take substantially more time than detection observations. Detection observations are designed to maximize detection probability while minimizing integration time (Keithly et al. 2020). HabEx Study Team (2019) considers a planet to be detected when the planet has signal-to-noise ratio greater than 7. These observations are planned to be made with coronagraphs, over a portion of the visible spectrum (565 ± 56 nm; Keithly et al. 2020) where most planet–star flux ratios are largest (Madden & Kaltenegger 2018). Between three and four detection observations are required to fit a single detected planet's semimajor axis, eccentricity, and inclination with sufficient certainty to identify the planet as within the Habitable Zone and justify a follow-up spectral characterization (Guimond & Cowan 2019; Horning et al. 2019). HabEx Study Team (2019) optimizes and plans observations of Earth analogs in isolation, assuming no other planets exist in the system or that the Earth-like planet can be distinctly identified from the other planets when it is first detected. Considering other non-Earth-like planets like this introduces a planet classification confusion problem when other planets are included as exemplified by the Neptune–Earth false-positive Monte Carlo study in Guimond & Cowan (2018).

In this Letter, we focus on the first direct image of solar system planets where only estimates of the planet–star separation (s) and planet–star difference in magnitude (Δmag) for detected planets are measured. We assume each of the exoplanets detected are spatially resolved. Even if the limited information collected indicates an exoplanet has an (s, Δmag) characteristic of an Earth-like exoplanet, it is not possible to discern with certitude that the exoplanet is Earth-like. While orbit fitting of multiple simultaneously detected and resolved exoplanets could preclude this, it is possible for multiple exoplanets to be detected and the Earth-like exoplanet be indiscernible from non-Earth-like exoplanets.

While the set of potential planets around a host star has a broad diversity and the habitability classification is broadly defined, we choose to study the classification confusion problem of our own solar system treated as an exosystem. We are motivated to do this for multiple reasons. The decadal survey seeks insight and answers as to how our solar system was formed and how it fits into the vast collection of other planetary systems (National Research Council 2011a). The abundance or rarity of planetary systems similar to the solar system is currently unknown (National Research Council 2011b; Committee on Exoplanet Science Strategy 2019). Future exoplanet direct-imaging missions will have the capability to discover solar-system-like star systems (HabEx Study Team 2019). Engineers are motivated to design future telescopes to the most strict requirements stemming from the challenges of detecting Earth-like exoplanets, meaning many other types of exoplanets in the star system will also be detected (HabEx Study Team 2019). Third, we know more about the solar system planets than any other planets in the galaxy thanks to the multitude of missions and studies of these bodies.

For this analysis, we make several simplifying assumptions about the solar system planets. First, we assume the planets are spherical, allowing us to use their volumetric mean radius (R). Second, we assume the combination of the geometric albedo (p) and phase function to be sufficient to describe the fraction of incident light reflected. We use the high-order polynomial model fit planetary phase functions from Mallama & Hilton (2018) to account for the unique reflective properties of each body. Third, we assume the planets have circular orbits. As we can see from Table 1, with the exception of Mercury, the eccentricities (e) are small (≲0.05). An eccentricity of 0.05 would change apoastron and periastron by ∼5%, thus maximally changing the net planetary flux by ∼10%. Assuming circular orbits allows substantial simplification of the s and Δmag functions. Fourth, we assume all planets lie in a common system angular momentum plane. Table 1 contains the inclinations of solar system planets from the ecliptic, each varying by a few degrees. This would affect some (s, Δmag)-coincidences, particularly between interior and exterior planets. We include a summary of the planet parameters used in this work in Table 1.

Using this model for the solar system, we will show that an external observer would find that many of our planetary bodies have multiple points of (s, Δmag)-coincidence along their orbits. In Section 2 we present the underlying planet photometric and astrometric models and how we combine them to get continuous phase curves spanning the entire range of phase angles. Appendix contains the melded planet phase curves derived from Mallama & Hilton (2018). In Section 3 we show our process for finding the locations of (s, Δmag)-coincidences and inclination deviations from edge-on where intersections still occur. Finally, in Section 4, we calculate the fraction of solar systems where a given planet has (s, Δmag)-coincidence with Earth.

Our goal is to find the fraction of inclined solar systems that could have (s, Δmag)-coincidence between any two planets. To do this, we first need a method for finding the (s, Δmag)-coincidence points between any two solar system planets.

We start with the general equation for Δmag given in Equation (3) of Brown (2005). By assuming the orbits are circular (e = 0), this simplifies to

$$\begin{eqnarray}&&{\rm{\Delta }}\mathrm{mag}=-2.5{\mathrm{log}}_{10}\left(p{\left(\displaystyle \frac{R}{a}\right)}^{2}{\rm{\Phi }}(\beta )\right).\end{eqnarray} \tag{ 1 }$$

The geometric albedo (p) and planetary volumetric mean radius (R) can be substituted in from Table 1 for each planet. This results in Δmag as a function of the planet phase function (Φ) and phase angle (β).

The phase angles are limited by the common system inclination (i). The global β extrema are in the edge-on (i = 90°) system, but the β extrema for any given inclination are

$$\begin{eqnarray}&&{\beta }_{\mathrm{extrema}}=90^\circ \pm i.\end{eqnarray} \tag{ 2 }$$

In this work, we use the high-order, parametric, polynomial fit phase functions from Mallama & Hilton (2018) to cover a large portion of the phase angle space for most planets. Where the polynomial fit would be extrapolated beyond measured data and is therefore unreliable, we substitute in the Lambert phase function from Equation (4) of Brown (2005) first presented in Sobolev (1975). To make these parametric phase functions usable in a continuous optimization method, we "meld" the phase functions together using parameterized hyperbolic tangent functions. The smallest phase angle where melding between the planet's phase function and the Lambert phase function occurs is 130°, for Jupiter. The phase functions of inner solar system planets span nearly the entire range of phase angles. The limits of the model fit planet phase functions are included in the Appendix.

We selected the hyperbolic tangent function to join parametric phase function components together because it allows us to create a melded phase function with smooth transitions and is differentiable over the entire range. The $\tanh$ is a curve ranging from −∞ < x < ∞ where ${\mathrm{lim}}_{x\to \infty }\tanh (x)\,=1$, ${\mathrm{lim}}_{x\to -\infty }\tanh (x)=-1$, and $\tanh (x=0)=0$. We convert this into two separate equations used to transition between the "start" and "end" of an individual parametric model. The first is modified to ensure ${\mathrm{lim}}_{x\to -\infty }f(x)=0$ and ${{\rm{lim}}}_{x\to \infty }f(x)=1,$ used to "start" a model. The second is modified to ensure ${\mathrm{lim}}_{x\to -\infty }f(x)=1$ and ${{\rm{lim}}}_{x\to \infty }f\left(x\right)=0,$ used to "end" a model. We also add two constants unique to each planet in order shift the model transition midpoint with A and adjust the transition slope with B to get

$$\begin{eqnarray}&&{f}_{\mathrm{start}}(x,A,B)=0.5+0.5\tanh \left(\displaystyle \frac{x-A}{B}\right),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{f}_{\mathrm{end}}(x,A,B)=0.5-0.5\tanh \left(\displaystyle \frac{x-A}{B}\right).\end{eqnarray} \tag{ 4 }$$

The melding of phase functions introduces small, negligible errors at β = 0° and β = 180°. The melded phase functions are shown in Figure 1 along with the typically used Lambert phase function.

Zoom In Zoom Out Reset image size

Brown (2005) took the general planet–star separation equation and reduced it into a function of phase angle, by adding our circular orbit simplification we get

$$\begin{eqnarray}&&s=a\sin (\beta ).\end{eqnarray} \tag{ 5 }$$

We plot the Δmag versus s curves for each solar system planet in Figure 2 at varying inclinations with HabEx angular measurement uncertainty of σ_WA = 0.05 mas at 10 pc and Δmag uncertainty derived from HabEx Study Team (2019) to be σ_Δmag = 0.145. This photometric uncertainty assumes a signal-to-noise ratio of 7 will be achieved on every planet and is achievable across the entire Δmag range of the planets. We plot the measurement uncertainty bounds by sampling Δmag and s over β and plotting the resulting 1σ bounds in Figure 2(a) and 3σ in Figures 2(b)–(f). The Δmag versus s curves in Figure 2(a) show 21 different locations where planet pairs have coincidence in the edge-on system. As the system inclination changes from edge-on to face-on in Figures 2(b)–(f), we see the number of intersections and visible range of the planets decrease. These Δmag versus s curves give us the core components required to find the (s, Δmag)-coincidence for any give planet pair.

Zoom In Zoom Out Reset image size

We cannot solve for (s, Δmag)-coincidence using conventional root-finding techniques because the underlying equations are nonlinear, their combination is underconstrained, and we do not know good initial guesses of the (s, Δmag)-coincidence points. This problem has three degrees of freedom: the phase angle of each planet (β_s and β_l for the interior and exterior planet, respectively) and the inclination of the common orbital plane. We only have the Δmag and s constraint equations making this problem underconstrained. To circumvent this issue, we formulate the problem as a constrained minimization problem and include i as a constraint on β_s and β_l .

In our optimization formulation, we minimize the absolute difference between the Δmag values of each planet, Δmag error. We additionally require both planets have the same planet–star separation. For each planet–star separation, there are two associated Δmag values. The larger Δmag, and therefore the dimmer side, occurs where β > 90°. The smaller Δmag, and therefore the brighter side, occurs where β < 90°. When optimizing, it is not possible for β_s or β_l to cross the infinite slope point of the Δmag versus s curve. We therefore formulate four separate optimization initial conditions and phase angle constraints associated with the portions of the phase curve the interior and exterior planets could have coincidence on: where the interior–exterior planets are brighter–dimmer, brighter–brighter, dimmer–dimmer, and dimmer–brighter, which we label as Q ∈ {0, 1, 2, 3}. Therefore, there are four sets of initial guesses of (β_s,0,β_l,0) to test the optimization process, which we differentiate with Q ∈ {0, 1, 2, 3}:

$$\begin{eqnarray}({\beta }_{s,0},{\beta }_{l,0})=\left\{\begin{array}{ll}\left(\displaystyle \frac{90+i}{2},\left({\sin }^{-1}\left(\displaystyle \frac{{a}_{s}}{{a}_{l}}\right)\displaystyle \frac{180}{\pi }+90\right)\times 0.3+(180-i)\times 0.7\right), & Q=0\\ \left(\displaystyle \frac{90+i}{2},\displaystyle \frac{1}{2}\left(i+{\sin }^{-1}\left(\displaystyle \frac{{a}_{s}}{{a}_{l}}\right)\displaystyle \frac{180}{\pi }\right)\right), & Q=1\\ \left(90\times 0.3+(180-i)\times 0.7,\left({\sin }^{-1}\left(\displaystyle \frac{{a}_{s}}{{a}_{l}}\right)\displaystyle \frac{180}{\pi }+90\right)\times 0.3+(180-i)\times 0.7\right), & Q=2\\ \left(90\times 0.3+(180-i)\times 0.7,\displaystyle \frac{1}{2}\left(i+{\sin }^{-1}\left(\displaystyle \frac{{a}_{s}}{{a}_{l}}\right)\displaystyle \frac{180}{\pi }\right)\right), & Q=3.\end{array}\right.\end{eqnarray} \tag{ 6 }$$

The 0.3 and 0.7 values were selected to create a midpoint value between the minimum and maximum possible β for the interior and exterior planet.

Algorithm 1. Minimum Δmag error producing phase angles

Input: i, as , al , Rs , Rl , ps , pl , ${\beta }_{s,0}$, ${\beta }_{l,0}$, and Q

Output: ${\beta }_{s}^{* }$ and ${\beta }_{l}^{* }$, the optimal planet phase angles of the small as and large al planet respectively

$\mathop{\arg \,min}\limits_{{\beta }_{s},{\beta }_{l}}$ $\left|2.5{\mathrm{log}}_{10}\left({p}_{s}{\left(\displaystyle \frac{{R}_{s}}{{a}_{s}}\right)}^{2}{{\rm{\Phi }}}_{s}({\beta }_{s})\right)-2.5{\mathrm{log}}_{10}\left({p}_{l}{\left(\displaystyle \frac{{R}_{l}}{{a}_{l}}\right)}^{2}{{\rm{\Phi }}}_{l}({\beta }_{l})\right)\right|$

${\beta }_{s}^{\ast },{\beta }_{l}^{\ast }\,=$

${\rm{s}}.{\rm{t}}.$

${a}_{s}\sin ({\beta }_{s})-{a}_{l}\sin ({\beta }_{l})=0$

$\,-{\beta }_{s}\leqslant -i\,\mathrm{If}\ Q=0,1,$

$\,{\beta }_{s}\leqslant 180-i\,\mathrm{If}\ Q=2,3,$

$\,-{\beta }_{l}\leqslant -{\sin }^{-1}\left(\displaystyle \frac{{a}_{s}}{{a}_{l}}\right)\displaystyle \frac{180}{\pi }-90\,\mathrm{If}\ Q=0,2,$

$\,{\beta }_{l}\leqslant 180-i\,\mathrm{If}\ Q=0,2,$

$\,-{\beta }_{l}\leqslant -i\,\mathrm{If}\ Q=1,3,$

$\,{\beta }_{l}\leqslant {\sin }^{-1}\left(\displaystyle \frac{{a}_{s}}{{a}_{l}}\right)\displaystyle \frac{180}{\pi }\,\mathrm{If}\ Q=1,3\quad$

Download table as:  ASCIITypeset image

We run Algorithm 1 over each unique pair of solar system planets and each constraint associated with Q ∈ {0, 1, 2, 3}. Some of these optimization processes do not successfully terminate. This occurs when the inclination constraints do not allow the separation constraint to be satisfied. While other optimization formulations may successfully terminate via convergence to a minimum error solution, not all minimum error solutions are locations of (s, Δmag)-coincidence. We apply a threshold, defining planet pairs with ∣Δmag_s − Δmag_l ∣ < 10⁻⁵ as coincident. This selection forces the omission of the Mars–Uranus intersection despite the two planets having a substantial region of overlapping uncertainty, but still includes the Earth–Saturn intersection. Table 2 contains the s and Δmag of coincidence as well as the phase angles of the interior and exterior planet. Only the Mars–Jupiter pair has two (s, Δmag)-coincidence points. There are only five instances where (s, Δmag)-coincidence occurs over a phase angle region using the Lambert phase function. They occur for Jupiter at Jupiter–Neptune, Jupiter–Uranus, Mars–Jupiter (1), and Mars–Jupiter (2) as well as Saturn–Neptune.

With these phase angles of coincidence, we can calculate the deviations from edge-on inclinations where the solar system no longer has (s, Δmag)-coincidence for any given planet pair. We define these critical inclination deviations as (δ i_crit,±). They can be found by

$$\begin{eqnarray}&&\begin{array}{l}\delta {i}_{\mathrm{crit},\pm }=90^\circ \pm \min \left[\min \left({\beta }_{s},180^\circ -{\beta }_{s}\right)\right.,\\ \left.\min \left({\beta }_{l},180^\circ -{\beta }_{l}\right)\right].\end{array}\end{eqnarray} \tag{ 7 }$$

We consider the Earth's Δmag versus s curve in Figure 2, which we take as representative of the highest scientific priority exoplanet type. The Earth's curve crosses those of Mercury, Neptune, Uranus, Mars, Venus, and Saturn (within the 3σ uncertainty region). Furthermore, the critical inclinations in Table 3 indicate (s, Δmag)-coincidence persists across a broad range of inclinations.

When simulating a multitude of inclined star systems, we randomly sample inclinations such that the probability density function is ${f}_{\overline{i}}(i)=\sin (i)/2$ following after Brown (2005), Savransky et al. (2009), and Keithly et al. (2020). Since we are assuming circular orbits, we can calculate the percentage of randomly generated solar systems where any solar system planet will have (s, Δmag)-coincidence with the Earth by integrating over the inclination probability density function,

$$\begin{eqnarray}&&P(\delta {i}_{\mathrm{crit},-}\leqslant i\leqslant \delta {i}_{\mathrm{crit},+})=\displaystyle \frac{1}{2}{\int }_{\delta {i}_{\mathrm{crit},-}}^{\delta {i}_{\mathrm{crit},+}}\sin (i)\delta i.\end{eqnarray} \tag{ 8 }$$

This distribution of inclinations peaks at i = 90° (edge-on) and has minimums at i = 0° and i = 180°. We compute these probabilities and compile them in Table 3.

Table 3 shows the fraction of solar systems with (s, Δmag)-coincidence independent of instrument capabilities. We can determine which intersecting planet pairs are visible to an instrument by comparing the Δmag and s of intersection with the instrument limited Δmag (${\rm{\Delta }}{\mathrm{mag}}_{\mathrm{lim}}$) and working angle limits at a star distance. For the Earth's (s, Δmag)-coincidence referenced in Table 3, an instrument with a contrast of 10⁻¹⁰ and inner working angle of 45 mas at 10 pc would only be able to see Earth's coincidence with Venus. As we can see from Table 2, increasing the limiting Δmag of a telescope exacerbates the potential for planet-type confusion by including more instances of (s, Δmag)-coincidence. Assuming the same instrument as in Table 3, only the Earth–Venus, Venus–Saturn, and Jupiter–Saturn coincidences are detectable. If, at some point in the future, a contrast of 10⁻¹¹ with an inner working angle of 30 mas were achievable; then 15 intersections could be observable in solar system analogs.

Assuming a system produces (s, Δmag)-coincidence between two planets, we can compute the probability a planet randomly located along its orbit is within n σ of the coincidence point. We randomly sample inclinations from ${f}_{\overline{i}}(i)$ between the critical inclination limits and randomly distribute these planets uniformly in time along their orbit. s and Δmag of each randomly sampled planet can be computed and we can find the fraction of these planets within the n σ uncertainty region of the (s, Δmag)-coincidence point. This fraction is the fraction of an orbit the planet spends within this uncertainty region and is converted into the percentages in Table 2. In all cases, the probability of the larger semimajor axis planet being within the uncertainty bounds is less than the probability of the smaller semimajor axis planet being within the measurement uncertainty bounds (P_{n σ l} < P_{n σ s} ). The maximum probability of coincidence is between Venus–Earth followed by Mercury–Venus, Mercury–Earth, and Earth–Mars. If a planet from a solar-system-like star system at the (s, Δmag)-coincidence of Earth with either Uranus or Neptune is detected, it is ∼150× more likely that the planet is an Earth than a Uranus or Neptune. In general, the probabilities of either planet being within the instrument uncertainty bounds are within an order of magnitude of one another. Only the instances of Mercury–Uranus, Mercury–Neptune, Venus–Saturn, Earth–Uranus, and Earth–Neptune have occurrence disparities greater than an order of magnitude.

Future exoplanet direct-imaging missions must make multiple observations to differentiate between Earth-like exoplanets and the myriad of other planets in the population. We took phase functions derived from a variety of deep space missions to create melded phase functions of solar system planets. We showed that up to 21 cases of (s, Δmag)-coincidence between planet pairs exist in the solar system. We additionally showed how an Earth can have the same (s, Δmag)-coincidence with up to six other solar system planets. We found the inclination range where each solar system planet could still have coincidence with another planet. We found 36%–69% of inner solar system planets and 1%–4% of outer solar system planets share (s, Δmag)-coincidence with Earth. While the Nancy Grace Roman Space Telescope would only be capable of seeing coincidences between Earth and Venus, further improvement in instrument contrast and inner working angles will exacerbate the planet confusion problem.

This work was funded by the Science Investigation Team of the Nancy Grace Roman Space Telescope under NASA grant NNX15AB40G.

Here, we combine the phase functions extracted from the model fit visual magnitude functions (V_mag) of each solar system planet in Mallama & Hilton (2018) using the parameterized hyperbolic tangent functions to meld parametric phase functions into a continuous phase function. In β ranges where a phase curve model for the planet is unavailable, we "fill in the gaps" by substituting in the Lambert phase function (Equation (4) from Brown 2005).

In general, we have knowledge of the phase function over most of the β range for Earth and planets interior to Earth. The phase date we have for planets exterior to Earth are limited to the phase angles observable by Earth and the ranges imaged via various flybys. The regions generally requiring melding with the Lambert phase function to span the β range are in excess of 130° and are where the planet is dimmest and least detectable.

Mercury's phase function is given by

$$\begin{eqnarray}\begin{array}{rcl}-2.5{\mathrm{log}}_{10}\left({{\rm{\Phi }}}_{\unicode{x0263F}}(\beta )\right) & = & 6.3280\times {10}^{-2}\beta -1.6336\times {10}^{-3}{\beta }^{2}\\ & & +\,3.3644\times {10}^{-5}{\beta }^{3}-3.4265\times {10}^{-7}{\beta }^{4}\\ & & +\,1.6893\,\times \,{10}^{-9}{\beta }^{5}-3.0334\times {10}^{-12}{\beta }^{6}.\end{array}\end{eqnarray} \tag{ A1 }$$

It spans the entire β range and does not require melding.

Venus's phase function is separated over two regions. The first,

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02640},1}(\beta )\\ =\,{10}^{-0.4(-1.044\times {10}^{-3}\beta +3.687\times {10}^{-4}{\beta }^{2}-2.814\times {10}^{-6}{\beta }^{3}+8.938\times {10}^{-9}{\beta }^{4})},\end{array}\end{eqnarray} \tag{ A2 }$$

is defined from 0° ≤ β ≤ 1637. The second,

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02640},2}(\beta )\\ \quad =\,{10}^{-0.4(-2.81914\beta +8.39034\times {10}^{-3}{\beta }^{2}))},\end{array}\end{eqnarray} \tag{ A3 }$$

is defined over the region 1637 ≤ β ≤ 179°. The parametric phase functions have minor discontinuities that we account for by adding a constant to and scaling the second parametric phase function into

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02640},3}(\beta ) & = & {{\rm{\Phi }}}_{\mathrm{int},\unicode{x02640},1}(163\buildrel{\circ}\over{.} 7)\\ & & \times \,\left(1+\displaystyle \frac{{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02640},2}(\beta ){{\rm{\Phi }}}_{\mathrm{int},\unicode{x02640},1}(163\buildrel{\circ}\over{.} 7)}{{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02640},2}(163\buildrel{\circ}\over{.} 7)-{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02640},2}(179^\circ )}\right.\\ & & -\left.\displaystyle \frac{{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02640},2}(163\buildrel{\circ}\over{.} 7){{\rm{\Phi }}}_{\mathrm{int},\unicode{x02640},2}(163\buildrel{\circ}\over{.} 7)}{{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02640},2}(163\buildrel{\circ}\over{.} 7)-{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02640},2}(179^\circ )}\right).\end{array}\end{eqnarray} \tag{ A4 }$$

Finally, we meld Φ_int,♀,1(β) and Φ_int,♀,3(β) together using the hyperbolic tangents to arrive at Venus's melded phase function

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{\unicode{x02640}}(\beta ) & = & {f}_{\mathrm{end}}(\beta ,A=163\buildrel{\circ}\over{.} 7,B=5^\circ ){{\rm{\Phi }}}_{\mathrm{int},\unicode{x02640},1}(\beta )\\ & & +\,{f}_{\mathrm{start}}(\beta ,A=163\buildrel{\circ}\over{.} 7,B=5^\circ )\\ & & \times \,{f}_{\mathrm{end}}(\beta ,A=179^\circ ,B=0\buildrel{\circ}\over{.} 5){{\rm{\Phi }}}_{\mathrm{int},\unicode{x02640},3}(\beta )\\ & & +\,{f}_{\mathrm{start}}(\beta ,A=179^\circ ,\\ B & = & 0\buildrel{\circ}\over{.} 5){{\rm{\Phi }}}_{L}(\beta )+2.766\times {10}^{-4}.\end{array}\end{eqnarray} \tag{ A5 }$$

The phase function for Earth spans the range of β and is

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\oplus }(\beta )={10}^{-0.4(-1.060\times {10}^{-3}\beta +2.054\times {10}^{-4}{\beta }^{2})}.\end{eqnarray} \tag{ A6 }$$

Mars's phase function is separated over two regions. The first,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02642},1}(\beta )={10}^{(-0.4(0.02267\beta -0.0001302{\beta }^{2}))},\end{eqnarray} \tag{ A7 }$$

is valid from 0° ≤ β ≤ 50°. The second,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02642},2}(\beta )={10}^{(-0.4(-0.02573\beta +0.0003445{\beta }^{2}))},\end{eqnarray} \tag{ A8 }$$

is valid over 50° < β ≤ 180°. We need to normalize Φ_int,♂,2(β) to account for discontinuities between the two phase functions and arrive at the corrected second phase function of

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02642},3}(\beta )={{\rm{\Phi }}}_{\mathrm{int},\unicode{x02642},1}(50^\circ )/{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02642},2}(50^\circ ){{\rm{\Phi }}}_{\mathrm{int},\unicode{x02642},2}(\beta ).\end{eqnarray} \tag{ A9 }$$

Combining these two phase functions, we get Mars's melded phase function

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{\unicode{x02642}}(\beta ) & = & {f}_{\mathrm{end}}(\beta ,A=50^\circ ,B=5^\circ ){{\rm{\Phi }}}_{\mathrm{int},\unicode{x02642},1}(\beta )\\ & & +\,{f}_{\mathrm{start}}(\beta ,A=50^\circ ,B=5^\circ ){{\rm{\Phi }}}_{\mathrm{int},\unicode{x02642},3}(\beta ).\end{array}\end{eqnarray} \tag{ A10 }$$

In the original formulation from Mallama & Hilton (2018), both parametric V_mag functions of Mars vary depending upon the rotational and orbital longitude of the planet. We assume the average of these correction terms, which are both 0.

The phase function for Jupiter is defined over two separate regions. The first is

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02643},1}(\beta )\,=\,{10}^{(-0.4(-3.7\times {10}^{-4}\beta +6.16\times {10}^{-4}{\beta }^{2}))}.\end{array}\end{eqnarray} \tag{ A11 }$$

The second is

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02643},2}(\beta )=1-1.507(\beta /180^\circ )\\ \quad -\,0.363{\left(\beta /180^\circ \right)}^{2}-0.062{\left(\beta /180^\circ \right)}^{3}\\ \quad +\,2.809{\left(\beta /180^\circ \right)}^{4}-1.876{\left(\beta /180^\circ \right)}^{5}.\end{array}\end{eqnarray} \tag{ A12 }$$

We offset Φ_int,♃,2(β) accounting for the the last valid value of Φ_int,♃,1(β) to arrive at the modified second phase function of

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02643},3}(\beta )={{\rm{\Phi }}}_{\mathrm{int},\unicode{x02643},1}(12^\circ )\\ \quad -\,{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02643},2}(12^\circ )+{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02643},2}(\beta ).\end{array}\end{eqnarray} \tag{ A13 }$$

We then combine the phase functions for Jupiter to get

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{\unicode{x02643}}(\beta ) & = & {f}_{\mathrm{end}}(\beta ,A=12^\circ ,B=5^\circ ){{\rm{\Phi }}}_{\mathrm{int},\unicode{x02643},1}(\beta )\\ & & +\,{f}_{\mathrm{start}}(\beta ,A=12^\circ ,B=5^\circ )\\ & & \times \,{f}_{\mathrm{end}}(\beta ,A=130^\circ ,B=5^\circ ){{\rm{\Phi }}}_{\mathrm{int},\unicode{x02643},3}(\beta )\\ & & +\,{f}_{\mathrm{start}}(\beta ,A=130^\circ ,B=5^\circ ){{\rm{\Phi }}}_{L}(\beta ).\end{array}\end{eqnarray} \tag{ A14 }$$

The phase function for Saturn is complicated by measurements of the planet obfuscated and augmented by the rings that have a unique phase function from the planet. The phase function of Saturn without the rings is defined over two separate regions. The first,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02644},1}(\beta )={10}^{(-0.4(-3.7\times {10}^{-4}\beta +6.16\times {10}^{-4}{\beta }^{2}))},\end{eqnarray} \tag{ A15 }$$

is based on Earth observations valid from 0° ≤ β ≤ 65. The second,

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02644},2}(\beta )\\ \,={10}^{(-0.4(2.446\times {10}^{-4}\beta +2.672\times {10}^{-4}{\beta }^{2}-1.505\times {10}^{-6}{\beta }^{3}+4.767\times {10}^{-9}{\beta }^{4}))}\end{array}\end{eqnarray} \tag{ A16 }$$

is defined over the range 6° ≤ β ≤ 150°. By properly adjusting the discontinuity between the first and second phase functions, we arrive at

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02644},3}(\beta ) & = & {{\rm{\Phi }}}_{\mathrm{int},\unicode{x02644},1}(6\buildrel{\circ}\over{.} 5)\\ & & -\,{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02644},2}(6\buildrel{\circ}\over{.} 5)+{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02644},2}(\beta ).\end{array}\end{eqnarray} \tag{ A17 }$$

We combine the phase functions for Saturn to get

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{\unicode{x02644}}(\beta ) & = & {f}_{\mathrm{end}}(\beta ,A=6\buildrel{\circ}\over{.} 5,B=5^\circ ){{\rm{\Phi }}}_{\mathrm{int},\unicode{x02644},1}(\beta )\\ & & +\,{f}_{\mathrm{start}}(\beta ,A=6\buildrel{\circ}\over{.} 5,B=5^\circ )\\ & & \times \,{f}_{\mathrm{end}}(\beta ,A=150^\circ ,B=5^\circ ){{\rm{\Phi }}}_{\mathrm{int},\unicode{x02644},3}(\beta )\\ & & +\,{f}_{\mathrm{start}}(\beta ,A=150^\circ ,B=5^\circ ){{\rm{\Phi }}}_{L}(\beta ).\end{array}\end{eqnarray} \tag{ A18 }$$

Uranus has a unique phase function due to the inclination of the pole's rotational axis, but only has one phase function. We assume a subsolar latitude of − 82°, which makes Uranus as bright as possible by adding a 0.0689 correction to the phase function (the minimum correction occurs at 82° and adds −0.0689). The phase function for Uranus is

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02645}}(\beta )={10}^{(-0.4(0.0689+6.587e-3\beta +1.045e-4{\beta }^{2}))},\end{eqnarray} \tag{ A19 }$$

valid over the range 0° ≤ β ≤ 154°. The combined phase function model for Uranus is

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{\unicode{x02645}}(\beta ) & = & {f}_{\mathrm{end}}(\beta ,A=154^\circ ,B=5^\circ ){{\rm{\Phi }}}_{\mathrm{ind},\unicode{x02645}}(\beta )\\ & & +\,{f}_{\mathrm{start}}(\beta ,A=154^\circ ,B=5^\circ ){{\rm{\Phi }}}_{L}(\beta ).\end{array}\end{eqnarray} \tag{ A20 }$$

Neptune has two phase functions, one based on Earth measurements that cover such a small range of phase angles that we ignore it and use the phase curve from Voyager 2 radiometer measurements. Neptune's one phase function is

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{int},\unicode{x02646}}(\beta )={10}^{-0.4(7.944\times {10}^{-3}\beta +9.617\times {10}^{-5}{\beta }^{2})},\end{eqnarray} \tag{ A21 }$$

and is valid over the range 0° ≤ β ≤ 13314. The combined phase function for Neptune is therefore

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{\unicode{x02646}}(\beta ) & = & {f}_{\mathrm{end}}(\beta ,A=133\buildrel{\circ}\over{.} 14,B=5^\circ ){{\rm{\Phi }}}_{\mathrm{int},\unicode{x02646}}(\beta )\\ & & +\,{f}_{\mathrm{start}}(\beta ,A=133\buildrel{\circ}\over{.} 14,B=5^\circ ){{\rm{\Phi }}}_{L}(\beta ).\end{array}\end{eqnarray} \tag{ A22 }$$