Analytic Light Curves in Reflected Light: Phase Curves, Occultations, and Non-Lambertian Scattering for Spherical Planets and Moons

Figure 1 shows a sample application of the algorithm developed in this paper: a reflected light phase curve of the Earth over the course of one year. The model is computed using the methodology in Appendix B from an l = 25 spherical harmonic expansion of the cloudless Earth, where the oceans are given an albedo of zero and the continents an albedo of unity (note, however, that since the light curve is normalized, the model does not depend on the value of the latter). The Earth is assumed to be a perfect Lambertian scatterer, so effects like the phase dependence of Rayleigh scattering and specular reflection (glint) from the oceans are neglected (but see Section 4 for an extension of the model to non-Lambertian scatterers). The observer is assumed to be along the ecliptic, so the illumination source is along the x–z plane of a right-handed Cartesian coordinate system, with $\hat{z}$ pointing toward the observer and $\hat{x}$ pointing to the right on the sky. The axis of rotation of the Earth is therefore tilted clockwise away from $\hat{y}$ by 235. The images at the top show snapshots of the disk of the Earth throughout the observation; below each one, we plot in blue the normalized phase curve at that phase over a single rotation. The orange dots correspond to a brute force numerical solution, obtained by discretizing the disk on a grid of ∼10⁵ points and summing over the dayside. The models agree to within the numerical precision of the brute force solution (about 100 ppm of the planetary flux in this case).

Zoom In Zoom Out Reset image size

While the dominant signal in the phase curve is the sine-like envelope due to the changing phases of the Earth, the local behavior of the light curve at each phase is complex and varies significantly over the course of the year. Unlike phase curves in thermal light, which primarily encode low-order spatial information (since the region of integration is always the full disk), phase curves in reflected light encode information at different scales depending on the phase. At crescent phase, the region of the disk contributing to the total flux is a narrow lune; these measurements therefore encode information primarily about high-l modes. At full phase, the region of integration is the full disk, so these measurements encode information about low-l modes. Furthermore, because of the obliquity of the Earth, the orientation of the crescent lune changes relative to features on the surface over the course of one orbit, changing the relative contribution of different portions of the surface to the flux and increasing the overall information content of the observation. As we will show in Section 5.2, the information content of reflected light phase curves is overwhelmingly higher than that of phase curves in thermal light, particularly for planets with significant obliquity.

Figure 2 shows another light curve of the rotating Earth, but this time taken during an occultation by the Moon. The map of the Earth is the same as before, but the observer is now along the equatorial plane of the Earth. The top panel is a reproduction of Figure 7 in Luger et al. (2019a) for the case of thermal light, where the Moon is seen to travel across the disk of the Earth from southwest to northeast, progressively occulting South America (dip), the Atlantic (peak), and Africa (dip). As before, the blue curve is the analytic solution and the orange dots correspond to the numerical solution.

Zoom In Zoom Out Reset image size

The bottom panel of the figure shows a light curve for the same occultation geometry, but seen instead in reflected light, with the Sun to the top left and slightly out of the page, corresponding to some point during northern summer. Note the same dip-peak-dip pattern, albeit with significantly different amplitudes. In particular, the transit across South America is deeper, because it occurs close to local noon, when the illumination is highest; conversely, the transit across Africa occurs close to local dusk, when the illumination is close to zero. As before, the light curves computed using starry agree to within the numerical precision of the brute force solution.

Our last sample light curve is Figure 3, which shows a secondary eclipse light curve of the Earth as it is occulted by the Sun. The model for the Earth is the same as above, and the observer is now close to the ecliptic, but slightly misaligned so that the Earth is occulted behind a solar latitude of 30° (i.e., at a solar impact parameter of 0.5). The observation takes place at the June solstice, so the Earth is tilted by 235 out of the page. As before, the top panel shows the light curve in thermal light; this is similar to the top panel of Figure 13 in Luger et al. (2019a). The orange dots again correspond to the numerical solution. The center panel shows the same light curve in reflected light. Because the observation occurs very close to full phase, the normalized light curves look very similar to each other. The bottom panel shows the difference between the two (reflected minus thermal), which is only on the order of a few percent. In fact, because the illumination profile is proportional to the cosine of the viewing angle, μみゅー, and the reflection is assumed to be isotropic, the (normalized) secondary eclipse light curve in reflected light is to good approximation equal to a limb-darkened thermal occultation light curve with linear limb darkening coefficient u₁ = 1. As we will see later, for very close-in planets, this approximation breaks down, since the illumination phases at secondary eclipse ingress and egress are sufficiently different from full phase.

Zoom In Zoom Out Reset image size

As we discuss in the Appendix, the model for phase curves and occultation light curves in reflected light may be expressed analytically in terms of purely algebraic and trigonometric functions and in some cases incomplete elliptic integrals of the first, second, and third kinds. We have derived efficient and numerically stable recursion relations to compute the relevant expressions and their derivatives. At times, these involve the evaluation of certain expressions numerically, especially when doing so leads to either a speed-up or a significant gain in numerical precision. In particular, as we discuss in Appendix C.1, the integration boundaries during an occultation sometimes depend on the solution to a quartic equation. While this can be solved in closed form, the analytic solution can often be very unstable. We therefore solve the quartic numerically, attaining a precision for the roots within a few orders of magnitude of machine (double) precision.

Figures 4 and 5 summarize the precision and computation time of the starry algorithm for two typical scenarios: a phase curve evaluation (Figure 4) and an occultation evaluation (Figure 5). Blue points correspond to the reflected light algorithm developed in this paper, while purple points correspond to the thermal light algorithm from Luger et al. (2019a) for the same occultation geometry, but without an illumination source. Solid and dashed lines correspond to evaluations without and with gradient propagation, respectively (see Section 3.3 for details). The orange and red dots correspond to numerical evaluation of the flux: brute force integration by summation on a grid of ∼10⁶ points (orange) and two-dimensional adaptive Gaussian quadrature using the dblquad function in scipy (Jones et al. 2001) with both absolute and relative error tolerances set to 10⁻³ (red). In both figures, the vertical axis corresponds to the evaluation time in seconds for a single flux computation, while the size of the points is proportional to the base-10 log of the relative error. For the starry solutions, the latter is estimated as the max-min difference in the flux over one thousand evaluations in which the input parameters are perturbed within an order of magnitude of machine epsilon; this is therefore a probe of the condition number of the starry algorithm and captures only error due to numerical instabilities. It is worth emphasizing that this is a measurement of the precision of the algorithm, rather than the accuracy, because it would be computationally intractable to compute a solution more accurate than this using a different algorithm. We argue that this measurement can be interpreted to mean that the algorithm is also accurate, but detailed quantification of this difference is beyond the scope of this paper. For the numerical solutions, the error is estimated as the difference between the numerical flux and the starry flux.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For both phase curves and occultations, the starry reflected light algorithm is 1–2 orders of magnitude slower than the emitted light algorithm, owing primarily to the increased complexity of the reflected light model. For phase curves (Figure 4), the thermal solution vector r ^⊤ (Equation (A7)) is a constant that can be precomputed, while the analogous vector in the reflected light case, ${{\mathbb{r}}}^{\top }$ (Equation (B2)), must be evaluated recursively each time. For occultations (Figure 5), the slower evaluation in the reflected light case is primarily due to the time spent solving the quartic equation for the points of intersection between the occultor and the day/night terminator of the illuminated body (Appendix C). This contributes the same overhead at all map degrees l, resulting in a gentler scaling in l than for the thermal case; for large l, the evaluation time for the two algorithms is within a factor of 2–3. In terms of precision, the algorithms are comparable, particularly for occultations. For both phase curves and occultations, the numerical error up to l = 10 is less than one part per trillion (10⁻¹²) for both thermal and reflected light curves.

Compared to either numerical evaluation method, the starry reflected light solutions are 4–5 orders of magnitude faster and about 10 orders of magnitude more precise. While different grid sizes and different error settings for the numerical integration change the balance slightly between these numbers, the starry solution is always many orders of magnitude faster and more precise than either method. In particular, because of the complicated integration boundaries (see, for example, Figure 15), two-dimensional Gaussian quadrature struggles to reach adequate accuracy in a reasonable amount of time, while integration on a grid fails to capture the curvature of the integration boundaries. Moreover, neither method yields the gradient of the solution with respect to the input parameters, which can be extremely useful for optimization and inference problems (see Section 3.3) below.

Note, importantly, that as we mentioned above, the reported error of the starry solution is only the numerical error of the algorithm: it does not capture any systematic error due to, say, an error in the derivation of the method. To this end, we rely on the Jupyter notebooks containing derivations and validations of the main equations in the Appendix, whose links (as footnote) appear below to the equation.We have also developed an extensive suite of unit tests comparing the starry solution to the numerical solution over a large grid of input parameter values, and verified that the solutions agree to within the precision of the numerical method. That said, there are specific cases in which the algorithm presented in the Appendix suffers from numerical instabilities. These generally happen due to division by small numbers or catastrophic cancellation in the recursion, and often occur near configurations involving grazing occultations, near-total occultations, terminator semiminor axis b ≈ 0 or b ≈ 1, etc. To mitigate these, we introduce various tolerance parameters in the code to either nudge the inputs away from these singular points or switch to a different evaluation method. These parameters are outlined in Table 5 at the end. In the vicinity of the cases described in that table, the precision of the starry algorithm will be reduced to (roughly) the value of the tolerance parameter, which in extremely rare cases can be as high as 10⁻⁵. ¹⁴

Table 5. Tolerance Parameters Used in the Code

Symbol	Description	Value
0	If ξくしー is this close to $\tfrac{n\pi }{2}$, compute ${{\mathbb{t}}}_{2}$ in the limit $\xi =\tfrac{n\pi }{2}$	10−13
1	If $| \sin \theta |$ or $| \cos \theta |$ are less than this value, set to this value	10−12
2	If k2 is within this value of unity, nudge it away	10−12
3	If ∣bo − ro∣ is less than this value, nudge bo away from ro	10−8
4	If bo is within this amount of ro − 1, nudge it away	10−8
5	If bo is within this amount of ro + 1, nudge it away	10−8
6	If $| \sin \alpha |$ is less than this value, set to this value	10−8
7	If b is within this value of zero, nudge it away	10−8
8	If bo is within this amount of 1 − ro, nudge it away	10−7
9	If two quartic roots are this close, eliminate one of them	10−7
10	If b is within this value of unity, set it to unity	10−6
11	If θしーた is within this amount of $\tfrac{\pi }{2}$ when ro = 1, nudge it away	10−5

Download table as:  ASCIITypeset image

The algorithm presented in this paper has been implemented in the Python package starry, which can be installed from GitHub or via the Python package manager pip. The algorithm is coded in a mixture of C++ with forward automatic differentiation using the Eigen library (Guennebaud & Jacob 2010) and Python with backward differentiation using just-in-time compiled theano operations (Theano Development Team 2016). The user interface, however, is purely in Python. The theano backend facilitates integration with the exoplanet modeling package (Foreman-Mackey et al. 2020) and in particular with pymc3 (Salvatier et al. 2016) for inference with gradient-based Markov Chain Monte Carlo (MCMC) schemes such as Hamiltonian Monte Carlo (HMC; Duane et al. 1987) and No-U-Turn Sampling (NUTS; Hoffman & Gelman 2011). Complete documentation and an extensive library of tutorials is available online. The links next to each of the figures (inside parenthesis) point to the Python scripts used to generate them and may also help in learning how to use starry. ^{15 , 16}

In the limit that the angular size of the star as seen from the planet is small, the illumination profile on the surface of the planet will decrease as the cosine of the angle between the surface normal and the star, reaching zero at the day/night terminator, an angle πぱい/2 away from the substellar point. However, if the star is sufficiently large and the planet is sufficiently close-in, rays originating from near the limb of the star will reach points on the planet surface beyond this angle. If the stellar radius R_⋆ is larger than the planet radius R_p, the angular extent of the true day/night terminator past πぱい/2 is given by

$$\begin{eqnarray}&&\tau =\arcsin \left(\displaystyle \frac{1-{R}_{{\rm{p}}}/{R}_{\star }}{a/{R}_{\star }}\right),\end{eqnarray} \tag{ 4 }$$

where a the semimajor axis of the orbit (where we implictly assume the eccentricity is zero). For planets like the Earth, this quantity is only about 026, resulting in a negligible effect on the planet's light curve. However, for very close-in planets, the effect can be significant (Knuth et al. 2017; Carter 2019). For instance, the hot Jupiter Kelt-9b has R_p/R_⋆ = 0.083 and a/R_⋆ = 3.16 planet (Wong et al. 2019). Assuming zero eccentricity and ignoring any stellar oblateness (see Ahlers et al. 2020), the day/night terminator extends τたう ≈ 17° past the limb of that planet. Figure 6 shows the illumination profile of this planet in a Mollweide projection, with the substellar point at the center, for the point-source approximation (left) and accounting for the finite size of the star (right). In addition to the displaced day/night terminator, the main difference between the two profiles is the substellar intensity, which is significantly higher in the extended source case. This is due to the simple fact that, in the point-source case, the illumination source is placed at the center of the star, which is one stellar radius farther from the planet than the point closest to the planet (the subplanetary point) in the extended source case. Once accounting for this difference, the fractional change in the intensity on the planet away from the substellar point is similar in both cases, and the intensity anywhere beyond πぱい/2 is less than one-tenth the peak value. ¹⁷

Zoom In Zoom Out Reset image size

The illumination profile in the extended source case may be computed as the two-dimensional integral of the point-source illumination profile over the visible portion of the stellar disk. While this integral may in theory be computed analytically (see, for instance, Kopal 1954, who derived series solutions to this problem), the resulting profile on the planet surface will not in general be exactly expressible in terms of spherical harmonics, a necessary condition for the starry algorithm. For simplicity, we therefore compute the illumination profile for extended sources by averaging the contribution of source_npts point sources uniformly distributed across the projected disk corresponding to the portion of the stellar surface visible from the planet, where source_npts is a user-supplied value. In the limit R_p ≪ R_⋆, this is a spherical cap centered at the subplanet point with radius ${R}_{\star }\cos \tau$. In Figure 6 we set source_npts = 300, but in practice we find that ∼30 points are sufficient for even the most extreme cases such as Kelt-9b. Note, importantly, that while this method allows one to account for the effect of stellar limb darkening on the illumination profile of the planet, this has not been implemented in starry.

Figure 7 shows the practical implications of the finite stellar size for the phase curve of Kelt-9b. The left panel shows the reflected light phase curve of the planet in parts per million, assuming a spherical albedo of 0.2, for the point-source approximation (blue) and including the effect of the extended source with source_npts = 300 (purple). The increased illumination at the substellar point results in a ∼30 ppm increase in the value of the phase curve close to full phase (±1/2). Close to a phase of zero, the extended source results in decreased flux, since the portion of the star illuminating the planet (the region close to the limb) is slightly farther away, by a factor of $\sqrt{1+{(a/{R}_{\star })}^{-2}}$. This results in a steeper phase curve, which can be seen in the right panel, where the light curves have been normalized to their maximum value.

Zoom In Zoom Out Reset image size

While modeling the extended size of the star is essential to getting the shape of the phase curve correct, the same is not true for secondary eclipse. Figure 8 shows the normalized secondary eclipse model for Kelt-9b under the point-source approximation (blue) and the extended source model (purple). Neglecting the fact that the depth of secondary eclipse is significantly different between the two models (see the left panel of Figure 7), the difference in shape between the two curves is almost negligible. For reference, the figure shows two additional models one might consider using to fit a secondary eclipse light curve: a uniform (unilluminated) disk (solid orange) and a disk whose intensity falls as μみゅー, the cosine of the viewing angle (dashed orange). Both can be computed using the classical Mandel & Agol (2002) transit model; the latter corresponds to a linearly limb-darkened sphere and is functionally equivalent to a Lambertian sphere seen at full phase. However, neither approximation is particularly good, since Kelt-9b changes illumination phase significantly from ingress to egress owing to its proximity to the star.

Zoom In Zoom Out Reset image size

The second assumption we now seek to relax is that of Lambertian scattering. A perfect Lambert sphere reflects light isotropically, so the measured intensity at a point on the surface is strictly proportional to the product of the cosine of the angle of incidence and the cosine of the viewing angle. While this is convenient from a modeling standpoint, it is hardly ever true in practice. For planets and moons in particular, there is often a strong phase dependence in the scattering. Rayleigh scattering in planetary atmospheres is preferentially in the forward/backward direction, while clouds and oceans can contribute strong specular reflection. Moreover, rough surfaces can have complex scattering behavior due to changes in the orientation of the surface normal on small scales and effects such as multiple reflections and self-shading.

In principle, any of these processes can be accounted for in the starry algorithm by modifying the linear operator I (Equation (A17)), which in the Lambertian case simply weights the spherical harmonic expansion of the albedo by the cosine-like illumination profile to obtain the observed intensity at a point on the surface (see Equations (A18) and (A20)). For non-Lambertian scattering, this matrix must also account for the phase dependence of the reflection: in particular, it will depend not only on the angle between the surface normal and the incident radiation, ϑ_i, but also on the angle between the surface normal and the reflected radiation (i.e., the direction toward the observer), ϑ_r. It may also depend on the azimuthal angles of the incident and reflected rays, ϕ_i and ϕ_r, respectively. These four angles are shown in Figure 9, showing the incoming radiation (source) vector s and the outgoing radiation (viewer) vector v in a frame in which the z-axis points along the surface normal.

Zoom In Zoom Out Reset image size

Treatment of a generalized, flexible scattering model is beyond the scope of this paper; see Heng et al. (2021) for recent results on this front. However, as an example of how a scattering model may be incorporated into the starry algorithm, we consider in detail the case of the rough surface scattering model of Oren & Nayar (1994), commonly used in computer graphics applications and solar system body modeling (e.g., Morgado et al. 2019). In this model, the surface is treated as a collection of a large number of Lambertian facets oriented at random angles relative to the average surface normal, whose net contribution to the total intensity can depart significantly from the Lambertian case. While the general model accounts for interreflections, shadowing, and an arbitrary distribution of facet orientations, in its simplest form the intensity observed at a point (x, y) on the (projected) surface of a body of unit spherical albedo may be approximated as (see Equation (30) in Oren & Nayar 1994)

$$\begin{eqnarray}&&{ \mathcal I }={{ \mathcal I }}_{\mathrm{Lamb}}\left\{{c}_{0}+{c}_{1}\,\max \left(0,\cos ({\phi }_{{\rm{r}}}-{\phi }_{{\rm{i}}})\right)\sin \alpha \tan \beta \right\},\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{c}_{0} & = & 1-0.5\left(\displaystyle \frac{{\sigma }^{2}}{{\sigma }^{2}+0.33}\right)\quad \alpha =\max \left({\vartheta }_{{\rm{r}}},{\vartheta }_{{\rm{i}}}\right)\\ {c}_{1} & = & 0.45\left(\displaystyle \frac{{\sigma }^{2}}{{\sigma }^{2}+0.09}\right)\quad \beta =\min \left({\vartheta }_{{\rm{r}}},{\vartheta }_{{\rm{i}}}\right),\end{array}\end{eqnarray} \tag{ 6 }$$

and the angles ϑ_i, ϑ_r, ϕ_i, and ϕ_r are all implicit functions of x, y, and the illumination source position. The term ${{ \mathcal I }}_{\mathrm{Lamb}}$ is the Lambertian illumination profile, given by Equation (2). At a given point on the surface, and for a given source position, the intensity ${ \mathcal I }$ is therefore a function of a single parameter, σしぐま, defined as the standard deviation in radians of the distribution of facet angles (which is assumed to be a zero-mean Gaussian). ^{21 , 22}

In order to incorporate this scattering model into starry, we must weight the spherical harmonic expansion of the albedo, y , by Equation (5) instead of Equation (2). Weighting by Equation (2) is (relatively) straightforward, since ${{ \mathcal I }}_{\mathrm{Lamb}}$ is a piecewise function of the l = 1 spherical harmonics (see Appendix A.2). The function we must integrate when computing fluxes is therefore exactly expressible in terms of spherical harmonics and thus starry-integrable. However, Equation (5) cannot be expressed exactly in terms of spherical harmonics, so we must instead approximate it. To this end, we evaluate Equation (5) on a grid of x and y spanning the unit disk, as well as the illumination phase, parametrized by b, the semiminor axis of the elliptical segment defining the day/night terminator (see Appendix A.2). We then fit to this a polynomial of total degree 5 in x, y, and $z\equiv \sqrt{1-{x}^{2}-{y}^{2}}$ and degree 5 in b and degree 4 in ${b}_{{\rm{c}}}\equiv \sqrt{1-{b}^{2}}$. Then, for a given value of b and b_c, we construct the operator I out of the polynomial coefficients in the same way as we constructed the Lambertian operator in Appendix A.2. More details about our approximation can be found in the Jupyter notebook accompanying Equation (5).

Figure 10 shows spheres of unit albedo with different surface roughness coefficients σしぐま and their corresponding phase curves. The sphere in the top row is perfectly Lambertian; its phase curve (blue) peaks at a value of 2/3, equal to the geometric albedo of a Lambert sphere. Increasing the surface roughness results in a greater relative contribution of flux from the limb of the object near full phase, since there now exist facets reflecting light directly back toward the observer (remaining rows and curves in the phase curve plot). Conversely, less light is scattered back to the observer at the subillumination point. These competing effects lead to phase curves that peak at a super-Lambertian value for small roughness coefficients (σしぐま = 15°, orange) and at a sub-Lambertian value for large roughness coefficients (σしぐま = 45°, red). We validate our calculations by computing the phase curves by numerically integrating the Oren & Nayar (1994) model over the visible disk; these are shown as the small dots in the figure, which agree to within 350 ppm of the body's flux for σしぐま = 45°.

Zoom In Zoom Out Reset image size

Our implementation of the scattering model extends just as easily to occultations and to cases where the surface does not have uniform albedo. Figure 11 shows a visible-light observation of the occultation of Io by Europa on 2009 December 4 taken by the PHEMU09 campaign (Arlot et al. 2014). The trajectory of Europa relative to Io, computed from the JPL Horizons database using ephemerides from Folkner et al. (2014) is shown at the top. We fit to this data a starry occultation model with the scattering law discussed above. For simplicity, we set the surface map equal too an l = 15 spherical harmonic expansion of the Galileo global color mosaic of Io (Becker & Geissler 2005) and use ephemerides from the JPL Horizons database, allowing for a small static x–y offset between Europa and Io as in Arlot et al. (2014) due to the uncertainty in the database. In total, we fit for five parameters: the two Cartesian offset terms, the flux contribution from Europa, the average albedo of Io, and the average surface roughness of Io, σしぐま. The model is displayed in blue and closely matches the data. Note, however, that this is meant simply as a demonstration of the starry algorithm, as our model for the surface is approximate at best, given differences in the wavelength band between the Galileo observations and those of the PHEMU09 campaign, changes in the albedo of Io since the Galileo measurements, and the fact that the orientation and extent of any shadows due to volcanoes on the surface are likely different between the Galileo and PHEMU09 observations. Furthermore, proper modeling would entail the joint analysis of all light curves of Io taken in a given season, for which we can afford to simultaneously fit for the surface map without risk of overfitting (Bartolić et al. 2022). ²⁵

Zoom In Zoom Out Reset image size

In the Appendix we derive closed-form expressions for the flux as a function of the spherical harmonic expansion of the albedo, y : Equation (A18) for occultations and Equation (A7) for phase curves. Inspection of those equations reveals that they are both linear in y : the flux is simply the dot product of several matrices and the vector of spherical harmonic coefficients. We may therefore write both expressions in the form

$$\begin{eqnarray}&&f={{\boldsymbol{x}}}^{\top }{\boldsymbol{y}},\end{eqnarray} \tag{ 7 }$$

where f is a scalar representing the model for the flux at a particular point in time and x ^⊤ is a row vector equal to ${{\mathbb{s}}}^{\top }{{\boldsymbol{A}}}_{{\bf{2}}}{\boldsymbol{I}}{{\boldsymbol{A}}}_{{\bf{1}}}{\boldsymbol{R}}^{\prime} {\boldsymbol{R}}$ (in the case of an occultation) or ${{\mathbb{r}}}^{\top }{\boldsymbol{I}}{{\boldsymbol{A}}}_{{\bf{1}}}{\boldsymbol{R}}^{\prime\prime} {\boldsymbol{R}}$ (in the case of a phase curve; see Appendix A.2 for details on what each of the terms represent). Now, if we let f be the vector of values of f for each point in the timeseries and construct the matrix X out of the stacked row vectors x ^⊤, we may write our model for the entire timeseries as the dot product

$$\begin{eqnarray}&&{\boldsymbol{f}}={\boldsymbol{X}}{\boldsymbol{y}}.\end{eqnarray} \tag{ 8 }$$

The linearity of the starry model is useful in several ways. For one, it can be exploited to cheaply compute the same model for different input vectors y . This is useful for multiband light curves, where the same matrix X dots into several vectors y , one for each observation band, or for time-dependent models, in which the model for the flux might be the Taylor series

$$\begin{eqnarray}&&{\boldsymbol{f}}={\boldsymbol{X}}\left({\left.{\boldsymbol{y}}(t)\right|}_{t={t}_{0}}+{\left.\displaystyle \frac{d{\boldsymbol{y}}^{\prime} (t)}{{dt}}\right|}_{t={t}_{0}}t+{\left.\displaystyle \frac{1}{2}\displaystyle \frac{{d}^{2}{\boldsymbol{y}}^{\prime} (t)}{{{dt}}^{2}}\right|}_{t={t}_{0}}{t}^{2}+\cdots \right),\end{eqnarray} \tag{ 9 }$$

where y (t) is a time-dependent representation of the surface map, which we expand about t = t₀ (Luger et al. 2019b); see also Kawahara & Masuda (2020) for an alternative linear model for time-dependent maps. But perhaps even more importantly, linear models are particularly useful for inference, since under Gaussian noise properties the posterior is analytic. In particular, if our light-curve measurements are given by the data vector d whose noise model is specified by the covariance matrix Σしぐま, and we place a Gaussian prior on y with mean μみゅー and covariance Λらむだ, the posterior mean may be written

$$\begin{eqnarray}&&\hat{{\boldsymbol{y}}}={\boldsymbol{C}}\left({{\boldsymbol{X}}}^{\top }{{\boldsymbol{\Sigma }}}^{-1}{\boldsymbol{d}}+{{\boldsymbol{\Lambda }}}^{-1}{\boldsymbol{\mu }}\right),\end{eqnarray} \tag{ 10 }$$

where C is the posterior covariance, given by

$$\begin{eqnarray}&&{\boldsymbol{C}}={\left({{\boldsymbol{X}}}^{\top }{{\boldsymbol{\Sigma }}}^{-1}{\boldsymbol{X}}+{{\boldsymbol{\Lambda }}}^{-1}\right)}^{-1}.\end{eqnarray} \tag{ 11 }$$

Because of this linearity, and the analyticity of the starry model, inference on data sets comprising thousands of points and spherical harmonic degree l ≤ 20 takes a fraction of a second on a typical computer. Full posterior inference with starry can thus be faster than the numerical evaluation of a single forward model (see Figure 5).

It is important to note, however, that the starry model is linear only in the map coefficients y . In any real application, there will be uncertainty in the inputs of X , such as the orbital parameters, the occultor parameters, the scattering law, etc. These parameters must typically be sampled over, since the model is a nonlinear function of them. However, the analyticity—and in particular, the differentiability—of the starry model makes sampling via gradient-based MCMC easy. Moreover, the linearity of the model with respect to y allows one to efficiently marginalize over those parameters when sampling over the nonlinear parameters. Tutorials on how to do this can be found in the starry documentation. ²⁷

One of the fundamental difficulties with the mapping problem is that the process of inferring a two-dimensional map from a light curve is almost always ill-posed. This has been known since at least the turn of the last century, when Russell (1906) discussed how, because of symmetry, all odd harmonics above l = 1 are in the null space for thermal phase curves of spherical bodies, meaning those terms do not contribute at all to the disk-integrated flux. As discussed in Luger et al. (2019a, 2021c), the problem is even more ill-posed than that: for any even degree l > 0, there are 2l + 1 modes on the surface (one for each value of m), but only 2 Fourier modes in the light curve (i.e., a sine and a cosine). Thus, for every mode that can be constrained from the light curve, there are far more modes that cannot, a problem that only gets worse as l increases.

The left panel of Figure 12 shows this issue in practice. We generate a mock thermal light curve (center left) from an l = 20 expansion of the Earth (top) with 1000 points over the course of one year with an exquisite photometric precision of 1 ppm. The Earth is given an obliquity of 235 on the plane of the sky but is viewed along the ecliptic, rotating edge-on with an inclination of 90°. The data are shown in black, and in blue is the posterior mean model (Equation (10)), in which we assume a prior variance of 10⁻³ for all spherical harmonic coefficients (and zero covariance). The corresponding surface map is shown at the bottom. As expected, this looks nothing like the true map of the Earth. For a body seen rotating edge-on, the information content of the light curve is strictly longitudinal. While the inferred map captures the average brightness of the Earth at each longitude fairly well, it is missing all latitudinal information. This is independent of the signal-to-noise or the cadence of the data set—it is a fundamental limitation of phase curves in thermal light. ²⁸

Zoom In Zoom Out Reset image size

The same is not true for the case of reflected light phase curves. In the right panel of Figure 12 we show the same mock light curve and perform the same inference step, but this time for observations in reflected light; the light curve is similar to that in Figure 1. Because of the presence of a day-night terminator beyond which features on the surface contribute zero flux, none of the symmetry arguments above apply. In particular, the facts that (1) the Earth is seen at different phases and (2) the terminator is inclined relative to the rotational axis mean that the region of the surface contributing to the phase curve is always changing, resulting in a complex light curve that encodes significantly more information than its thermal counterpart. The result is an inferred map that is largely faithful to the true map (bottom). While the continental outlines are somewhat fuzzy and some artifacts are present at high latitudes, it is clear that the inference problem is much less ill-posed in this case.

Unlike thermal light curves, reflected light phase curves have the potential to robustly constrain two-dimensional maps of exoplanets. This result is not new, and has been discussed at length in the literature (e.g., Fujii & Kawahara 2012; Berdyugina & Kuhn 2017; Luger et al. 2019b; Aizawa et al. 2020; Kawahara 2020). In particular, Kawahara & Fujii (2010) demonstrated the uniqueness of their inferred map from mock reflected light curves of the Earth. As we argued above, for specific geometrical configurations, the mapping problem in reflected light can actually be well-posed, meaning it has no null space up to a certain degree l. Reflected light phase curves of terrestrial planets with James Webb Space Telescope (JWST) and future direct imaging missions thus have the potential to reveal detailed information about their surfaces.

Nevertheless, at finite signal to noise and for limited observation duration or cadence, there may still be significant degeneracies in the reflected phase curve problem. In Luger et al. (2019a), we argued that occultations can be used to break many of these degeneracies, since they directly probe the surface at scales inaccessible to phase curves. The same is true in reflected light. Figure 13 demonstrates this for mock observations of the same Earth-like planet as in Figure 12, but this time taken over ten days near quadrature, when the disk of the planet is seen at half phase. The top rows show the thermal and reflected phase curves and the inferred maps, which look similar to those in the previous figure. At the bottom we show the same light curve, but this time including seven equatorial occultations by a moon one-quarter the size of the planet. The moon's period and occultation duration are unrealistically short, but the inferred maps at the bottom show the exquisite constraining power of these occultations. In the thermal case, the presence of the moon allows us to infer the two-dimensional distribution of surface features along its equatorial occultation path; however, there is little information in the light curve about features at higher latitudes, and the starry model overfits. Conversely, in the reflected light case, there is information about all latitudes and longitudes, and the inferred map is largely faithful to the true map.

Zoom In Zoom Out Reset image size

It should be kept in mind that the kinds of observations mentioned above will be very challenging for exoplanets, even with next-generation observatories such as HabEx or LUVOIR. Observations of occultations of planets by moons, in particular, are likely several decades away at least. There is some hope that planet–planet occultations may be detectable in the near future for specific planetary systems such as TRAPPIST-1 (Luger et al. 2019b), but at extremely limited signal-to-noise. Even phase curve observations will be difficult because of their limited signal-to-noise, and in practice many degeneracies will likely remain. Several studies have found that color information can greatly help in the interpretation of reflected phase curves (Cowan et al. 2009; Kawahara & Fujii 2011; Lustig-Yaeger et al. 2018; Kawahara 2020), while others have explored in detail the best kinds of priors to assume (e.g., Aizawa et al. 2020; Asensio Ramos & Pallé 2021). Future maps of exoplanets—particularly terrestrial ones in the habitable zone—will require every tool in the toolbox.

There are a few limitations to the starry model that are important to bear in mind. The primary limitation concerns the maximum spherical harmonic degree of the model. While the expressions derived here are valid at arbitrary degree, we find that their numerical stability quickly degrades above l ∼ 20–25 for occultations and l ∼ 35–40 for phase curves. The same is true for the model in thermal light (Luger et al. 2019a), and is due to (1) the large condition number of the change-of-basis matrix from spherical harmonics to polynomials and (2) instabilities in the many recursion relations used to evaluate the solution vectors ${{\mathbb{s}}}^{\top }$ and ${{\mathbb{r}}}^{\top }$. In principle, one could improve the numerical stability by evaluating all expressions at higher floating-point precision, but in practice the computational cost of this becomes quickly prohibitive. However, it is important to keep in mind that our current best image of an exoplanet is the l = 1 map of HD189733b (Knutson et al. 2007; de Wit et al. 2012; Majeau et al. 2012). Next-generation facilities such as JWST may allow us to probe surface modes as small as l = 5 for some planets (R. Luger et al. 2022, in preparation), but even with future telescopes such as LUVOIR, it is extremely unlikely we will do better than l = 20. If cases arise requiring a resolution smaller than about 180°/20 = 9° on the surface, the starry algorithm will have to be revisited.

The second limitation concerns the flexibility of the starry model. While we presented ways to capture non-Lambertian scattering in starry (Section 4.2), there are certain aspects of light curves in reflected light that cannot be captured by the model. One example of this is shadowing. Craters on the moon or volcanoes on Io can cast large shadows visible from space, particularly if viewed near crescent phase. Unfortunately, there is no way to model this within the starry framework. Another example is multiple scattering, as in optically thick atmospheres, for which a proper radiative transfer model must be used. There is also the case of specular reflection, or "glint," which is a pronounced signal for the Earth due its oceans (Robinson et al. 2014) and on Titan due to its hydrocarbon lakes (Barnes et al. 2011). In principle, glint could be modeled in the same way as non-Lambertian scattering, by constructing the linear operator I in such a way as to downweight portions of the projected disk where ϑ_i ≠ ϑ_r. In practice, however, if the size of the glint spot is small (which is typically the case), an expansion at extremely high l (l ∼ 360 in the case of the Earth) would be required, which would not work for the reasons above. Instead, it may be possible to combine the starry algorithm with the formalism of Haggard & Cowan (2018), who derived analytic expressions for phase curves of delta function maps (δでるた-map), to model glint. As discussed in Lustig-Yaeger et al. (2018), glint mapping is an extremely powerful way to not only map terrestrial planets but also to confirm their habitability via the presence of an ocean. Lastly, the presence of time-variable features, such as clouds, dust storms, or seasonal variations of vegetation are not accounted for in the model, although the documentation discusses how one may approach the modeling of temporal features. ³⁰

We would also like to emphasize that while spherical harmonics are a convenient basis for the purpose of computing light curves, they have a significant drawback when it comes to modeling real planetary surfaces: it can be difficult to strictly enforce physical values of the albedo everywhere on the surface map when doing inference. That is because there is no analytic way to determine whether a spherical harmonic representation is positive-valued (or restricted to a given range) everywhere on the sphere. Instead, this must be checked numerically, by evaluating the function on a discrete grid. This makes it somewhat cumbersome to implement positivity as a prior when doing inference; in particular, this prior cannot be expressed as a Gaussian, so the analytic expression for the posterior discussed in Section 5.1 will generally have nonzero support for negative albedo values. This is particularly problematic when the data is not very constraining and the posterior is prior-dominated, as positivity can be an extremely informative prior (e.g., Fienup 1978). We therefore recommend that in such cases users of the starry algorithm use HMC/MCMC to do inference, either (1) sampling over the spherical harmonic coefficients y and imposing a uniform prior in the range [0, 1] on the albedo values a evaluated on a discrete grid on the sphere or (2) sampling over the albedo values a with the same uniform prior, but using y to compute the actual light curve model (Bartolić et al. 2022). In both cases, there exists a linear operator that transforms between y and a :

$$\begin{eqnarray}&&{\boldsymbol{a}}={\boldsymbol{P}}{\boldsymbol{y}}\end{eqnarray} \tag{ 12 }$$

and

$$\begin{eqnarray}&&{\boldsymbol{y}}={{\boldsymbol{P}}}^{+}{\boldsymbol{a}},\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}&&{{\boldsymbol{P}}}^{+}={\left({{\boldsymbol{P}}}^{\top }{\boldsymbol{P}}+\lambda {\bf\mathsf{I}}\right)}^{-1}{{\boldsymbol{P}}}^{\top }\end{eqnarray} \tag{ 14 }$$

is the pseudoinverse of P , with (small) regularization parameter λらむだ and where I is the identity matrix. Each row of the matrix P is constructed from the value of each of the spherical harmonics at the corresponding point on the grid; both P and its inverse may be precomputed for efficiency. In both cases, the grid should be fine enough to ensure positivity over most of the sphere but not so fine as to throttle the computation; as a rule of thumb, we find that grids with ∼4 times as many pixels as spherical harmonic coefficients are sufficient. The documentation includes tutorials on how to implement this in practice. ³¹

Finally, while the algorithm presented here can be used to model planet–planet and planet–moon occultations in reflected light, we only account for the physical blocking of light rays from the planet by the occultor; we do not account for the attenuation of the reflected light due to the shadow of the occultor on the planet. This can be an important effect at certain orbital phases (e.g., Cabrera & Schneider 2007). However, since the shadow of the occultor is not necessarily circular, it is significantly harder to model in closed form. A proper treatment of this and other effects, such as occultations and shadows due to rings (e.g., Arnold & Schneider 2004), is deferred to future work.

Without loss of generality, assume the body whose flux we wish to compute has radius unity and sits at the origin of a right-handed Cartesian coordinate system in some frame ${{ \mathcal F }}_{0}$. In this frame, the surface (emitted) intensity field of the body is described by a vector y of coefficients in the spherical harmonic basis $\tilde{{\boldsymbol{y}}}$:

$$\begin{eqnarray}&&\begin{array}{l}\tilde{{\boldsymbol{y}}}(x,y)\\ \,=\,{\left(\begin{array}{cccccccccc}{Y}_{\mathrm{0,0}} & {Y}_{1,-1} & {Y}_{\mathrm{1,0}} & {Y}_{\mathrm{1,1}} & {Y}_{2,-2} & {Y}_{2,-1} & {Y}_{\mathrm{2,0}} & {Y}_{\mathrm{2,1}} & {Y}_{\mathrm{2,2}} & \cdot \cdot \cdot \end{array}\right)}^{\top },\end{array}\end{eqnarray} \tag{ A1 }$$

where the component at index n is the spherical harmonic Y_l,m (x, y) with

$$\begin{eqnarray}\begin{array}{rcl}l & = & \lfloor \sqrt{n}\rfloor \\ m & = & n-{l}^{2}-l.\end{array}\end{eqnarray} \tag{ A2 }$$

The spherical harmonics are traditionally expressed in spherical coordinates, but for our purposes it is more convenient to express them in Cartesian coordinates on the sky-projected disk, in which case they are simply polynomials in x, y, and z (see Appendix A in Luger et al. 2019a).

An observer views the body from a large distance in the sky frame ${ \mathcal F }$, in which the x-axis points to the right, the y-axis points up, and the z-axis points out of the sky toward the observer. Following Luger et al. (2019a), if an occultor of radius r_o is located at sky position (x_o, y_o), we compute the visible thermal flux f_T from

$$\begin{eqnarray}&&{f}_{{\rm{T}}}={{\boldsymbol{s}}}^{\top }{\boldsymbol{A}}{\boldsymbol{R}}^{\prime} {\boldsymbol{R}}{\boldsymbol{y}},\end{eqnarray} \tag{ A3 }$$

where, from right to left, R = R (I, Λらむだ, Θしーた) is a Wigner rotation matrix that rotates y from ${{ \mathcal F }}_{0}$ to the sky frame ${ \mathcal F }$ given the body's inclination I, obliquity Λらむだ, and rotational phase Θしーた (Appendix C in Luger et al. 2019a), ${\boldsymbol{R}}^{\prime} ={\boldsymbol{R}}^{\prime} ({x}_{{\rm{o}}},{y}_{{\rm{o}}})$ rotates the body on the plane of the sky into the integration frame ${ \mathcal F }^{\prime}$, in which the occultor lies along the $+y^{\prime}$-axis, A (Equation (B13) in Luger et al. 2019a) is the change-of-basis matrix from $\tilde{{\boldsymbol{y}}}$ to the Green's basis $\tilde{{\boldsymbol{g}}}$ in which the integrals are computed, whose component at index n is

$$\begin{eqnarray}&&\begin{array}{l}{\mathop{g}\limits^{}}_{n}(x,y)\\ =\,\left\{\begin{array}{ll}\displaystyle \frac{\mu +2}{2}{x}^{\tfrac{\mu }{2}}{y}^{\tfrac{\nu }{2}} & \mu ,\nu \,\mathrm{even}\\ z(x,y) & \mu =\nu =1\\ 3{x}^{l-2}{yz}(x,y) & \nu \,\mathrm{odd},\mu =1,\displaystyle \frac{\mu +\nu }{2}\,\mathrm{even}\\ z(x,y)\left(-{x}^{l-3}+{x}^{l-1}+4{x}^{l-3}{y}^{2}\right) & \nu \,\mathrm{odd},\mu =1,\,\mathrm{odd}\\ \begin{array}{l}z(x,y)\left(\displaystyle \frac{\mu -3}{2}{x}^{\tfrac{\mu -5}{2}}{y}^{\tfrac{\nu -1}{2}}-\displaystyle \frac{\mu -3}{2}{x}^{\tfrac{\mu -5}{2}}{y}^{\tfrac{\nu +3}{2}}\right.\\ \left.\qquad -\displaystyle \frac{\mu +3}{2}{x}^{\tfrac{\mu -1}{2}}{y}^{\tfrac{\nu -1}{2}}\right)\end{array} & \mathrm{otherwise},\end{array}\right.\end{array}\end{eqnarray} \tag{ A4 }$$

with

$$\begin{eqnarray}\begin{array}{rcl}\mu & \equiv & l-m\\ \nu & \equiv & l+m\end{array}\end{eqnarray} \tag{ A5 }$$

and

$$\begin{eqnarray}&&z(x,y)\equiv \sqrt{1-{x}^{2}-{y}^{2}},\end{eqnarray} \tag{ A6 }$$

and s ^⊤ = s ^⊤(b_o, r_o) is the vector of solutions to the integral over the projected visible disk of the body for each term in $\tilde{{\boldsymbol{g}}}$ (Equation (26) in Luger et al. 2019a), with ${b}_{{\rm{o}}}=\sqrt{{x}_{{\rm{o}}}^{2}+{y}_{{\rm{o}}}^{2}}$.

If instead no occultor is present, we compute the total visible thermal flux ${f}_{{{\rm{T}}}_{0}}$ from this body as

$$\begin{eqnarray}&&{f}_{{{\rm{T}}}_{0}}={{\boldsymbol{r}}}^{\top }{{\boldsymbol{A}}}_{{\bf{1}}}{\boldsymbol{R}}^{\prime\prime} {\boldsymbol{R}}{\boldsymbol{y}},\end{eqnarray} \tag{ A7 }$$

where, as before, R = R (I, Λらむだ, Θしーた) rotates the body from ${{ \mathcal F }}_{0}$ to the sky frame ${ \mathcal F }$, R '' rotates the body on the plane of the sky into the integration frame ${ \mathcal F }^{\prime\prime}$, A ₁ (Equation (B11) in Luger et al. 2019a) is the change-of-basis matrix from the spherical harmonic basis $\tilde{{\boldsymbol{y}}}$ to the polynomial basis $\tilde{{\boldsymbol{p}}}$ in which the integrals are computed, whose component at index n is

$$\begin{eqnarray}{\tilde{p}}_{n}(x,y)=\left\{\begin{array}{ll}{x}^{\tfrac{\mu }{2}}{y}^{\tfrac{\nu }{2}} & \mu ,\nu \,\mathrm{even}\\ {x}^{\tfrac{\mu -1}{2}}{y}^{\tfrac{\nu -1}{2}}z(x,y) & \mathrm{otherwise},\end{array}\right.\end{eqnarray} \tag{ A8 }$$

and r ^⊤ is the vector of solutions to the integral over the projected visible disk of the body for each term in $\tilde{{\boldsymbol{p}}}$ (Equation (19) in Luger et al. 2019a). ³²

In order to compute light curves in reflected light, we must make two modifications to the starry algorithm. First, the expressions above assume that the coefficient vector y describes the emissivity of the body, which (in the absence of limb darkening) is assumed to be Lambertian, i.e., all points on the surface emit equally in all directions. Here, we wish to derive the solution for the flux in the case of Lambertian reflectance, in which case the vector y is taken to describe the spherical albedo of the surface, A.

Second, we must explicitly model the illumination of the body. We assume the body is illuminated by a point-like source whose flux measured by the observer is unity. In this case, the observed intensity at any point on the surface is proportional to the cosine of the angle ϑ_i between the incident light and the surface normal. Points for which ϑ_i ≥ πぱい/2 are unilluminated and therefore have an intensity of zero. If the point-like illumination source is placed at sky coordinates (x_s, y_s, z_s) in units of the radius of the illuminated body, the day/night terminator on the body is a half-ellipse of semimajor axis unity that is fully described by its (signed) semiminor axis,

$$\begin{eqnarray}&&b=-\displaystyle \frac{{z}_{{\rm{s}}}}{{r}_{{\rm{s}}}},\end{eqnarray} \tag{ A9 }$$

where ${r}_{{\rm{s}}}=\sqrt{{x}_{{\rm{s}}}^{2}+{y}_{{\rm{s}}}^{2}+{z}_{{\rm{s}}}^{2}}$ is the distance to the source, and the angle by which its semimajor axis is rotated away from the + x-axis,

$$\begin{eqnarray}&&\theta =-\arctan 2({x}_{{\rm{s}}},{y}_{{\rm{s}}}),\end{eqnarray} \tag{ A10 }$$

where $\arctan 2(a,b)$ is the quadrant-aware arctangent of a/b. Given this formulation, and assuming that r_s ≫ 1, it is straightforward to show that the illumination ${ \mathcal I }$ at a point (x, y) on the projected disk of the body is given by the function

$$\begin{eqnarray}&&{ \mathcal I }(b,\theta ,{r}_{{\rm{s}}};x,y)=\max \left(0,I(b,\theta ,{r}_{{\rm{s}}};x,y)\right),\end{eqnarray} \tag{ A11 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}I(b,\theta ,{r}_{{\rm{s}}};x,y) & = & \displaystyle \frac{1}{\pi {r}_{{\rm{s}}}^{2}}\cos {\vartheta }_{{\rm{i}}}\\ & = & \displaystyle \frac{1}{\pi {r}_{{\rm{s}}}^{2}}\left(-{b}_{{\rm{c}}}\sin \theta x+{b}_{{\rm{c}}}\cos \theta y-{bz}(x,y)\right),\end{array}\end{eqnarray} \tag{ A12 }$$

with ${b}_{{\rm{c}}}\equiv \sqrt{1-{b}^{2}}$ and $z(x,y)=\sqrt{1-{x}^{2}-{y}^{2}}$. The illumination ${ \mathcal I }$ is a unitless quantity, normalized such that the integral of $A{ \mathcal I }$ over the unit disk is equal to the flux measured by the observer as a fraction of the flux of the illumination source.^, In particular, if we place the illumination source along the + z-axis at (0, 0, 1), the body is seen at full phase, so b = −1, b_c = 0, and

$$\begin{eqnarray}&&{{ \mathcal I }}_{\mathrm{full}}(x,y)=\displaystyle \frac{\sqrt{1-{x}^{2}-{y}^{2}}}{\pi }.\end{eqnarray} \tag{ A13 }$$

Multiplying this by the albedo and integrating over the unit disk, we obtain the reflected flux measured by the observer in units of the flux of the illumination source:

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\mathrm{full}} & = & {\displaystyle \int }_{-1}^{1}{\displaystyle \int }_{-\sqrt{1-{x}^{2}}}^{\sqrt{1-{x}^{2}}}A\displaystyle \frac{\sqrt{1-{x}^{2}-{y}^{2}}}{\pi }\,{dy}\,{dx}\\ & = & \displaystyle \frac{2}{3}A,\end{array}\end{eqnarray} \tag{ A14 }$$

which is precisely the geometric albedo of a Lambert sphere of spherical albedo A (see, e.g., Seager 2010). ^{33 , 34 , 35 , 36 , 37 , 38}

In principle, our task is now straightforward: weight each of the terms in the Green's basis (Equation (A4)) and integrate them over the visible portion of the body's disk to obtain the reflected light solution vector, ${{\mathbb{s}}}^{\top }$. Unfortunately, the piecewise nature of Equation (A11) makes direct evaluation of these integrals extremely difficult in practice. We find that it is more tractable to weight our basis terms by the function I (Equation (A12)) and to modify the limits of integration to exclude the nightside of the body, where I is (unphysically) negative. In particular, since I is just a polynomial in x, y, and z(x, y), we can express it as a vector i (b, θしーた) in the polynomial basis $\tilde{{\boldsymbol{p}}}$. Recalling the structure of the basis (Equation (A8)), we may write

$$\begin{eqnarray}&&{\boldsymbol{i}}(b,\theta ,{r}_{{\rm{s}}})=\displaystyle \frac{1}{\pi {r}_{{\rm{s}}}^{2}}\left(\begin{array}{c}0\\ -{b}_{{\rm{c}}}\sin \theta \\ -b\\ {b}_{{\rm{c}}}\cos \theta \end{array}\right).\end{eqnarray} \tag{ A15 }$$

This fact allows us to construct a linear operator I to weight a map vector in the polynomial basis by the illumination profile. If we think about how each of the terms in $\tilde{{\boldsymbol{p}}}$ transforms under I ,

$$\begin{eqnarray}&&\begin{array}{rcl}\left(\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right) & {\boldsymbol{\to }} & \left(\begin{array}{c}{{\boldsymbol{i}}}_{0}\\ {{\boldsymbol{i}}}_{1}\\ {{\boldsymbol{i}}}_{2}\\ {{\boldsymbol{i}}}_{3}\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right)\qquad \left(\begin{array}{c}0\\ 1\\ 0\\ 0\end{array}\right){\boldsymbol{\to }}\left(\begin{array}{c}0\\ {{\boldsymbol{i}}}_{0}\\ 0\\ 0\\ {{\boldsymbol{i}}}_{1}\\ {{\boldsymbol{i}}}_{2}\\ {{\boldsymbol{i}}}_{3}\\ 0\\ 0\end{array}\right)\\ \left(\begin{array}{c}0\\ 0\\ 1\\ 0\end{array}\right) & {\boldsymbol{\to }} & \left(\begin{array}{c}{{\boldsymbol{i}}}_{2}\\ 0\\ {{\boldsymbol{i}}}_{0}\\ 0\\ -{{\boldsymbol{i}}}_{2}\\ {{\boldsymbol{i}}}_{1}\\ 0\\ {{\boldsymbol{i}}}_{3}\\ -{{\boldsymbol{i}}}_{2}\end{array}\right)\qquad \left(\begin{array}{c}0\\ 0\\ 0\\ 1\end{array}\right){\boldsymbol{\to }}\left(\begin{array}{c}0\\ 0\\ 0\\ {{\boldsymbol{i}}}_{0}\\ 0\\ 0\\ {{\boldsymbol{i}}}_{1}\\ {{\boldsymbol{i}}}_{2}\\ {{\boldsymbol{i}}}_{3}\end{array}\right),\end{array}\end{eqnarray} \tag{ A16 }$$

we can compose I out of these column vectors:

$$\begin{eqnarray}{\boldsymbol{I}}(b,\theta ,{r}_{{\rm{s}}})=\displaystyle \frac{1}{\pi {r}_{{\rm{s}}}^{2}}\left(\begin{array}{ccccc}0 & 0 & -b & 0 & \cdots \\ -{b}_{{\rm{c}}}\sin \theta & 0 & 0 & 0 & \cdots \\ -b & 0 & 0 & 0 & \cdots \\ {b}_{{\rm{c}}}\cos \theta & 0 & 0 & 0 & \cdots \\ 0 & -{b}_{{\rm{c}}}\sin \theta & b & 0 & \cdots \\ 0 & -b & -{b}_{{\rm{c}}}\sin \theta & 0 & \cdots \\ 0 & {b}_{{\rm{c}}}\cos \theta & 0 & -{b}_{{\rm{c}}}\sin \theta & \cdots \\ 0 & 0 & {b}_{{\rm{c}}}\cos \theta & -b & \cdots \\ 0 & 0 & b & {b}_{{\rm{c}}}\cos \theta & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{array}\right),\end{eqnarray} \tag{ A17 }$$

where the dimensions of the matrix are $\left({(l+2)}^{2},{(l+1)}^{2}\right)$, where l is the spherical harmonic degree of the map (this operator raises the degree of the map by one). Note, again, that this weighting is valid only on the dayside hemisphere (see Equation (A11)), as the operator I weights points on the nightside by a negative amount, which is clearly unphysical. As we will see momentarily, we account for this by excluding the nightside from the integration region in our flux integrals. ^{40 , 41}

We may now re-write Equations (A3) and (A7) to account for this illumination transformation. The flux during an occultation is now given by

$$\begin{eqnarray}&&f={{\mathbb{s}}}^{\top }(b,\theta ^{\prime} ,{b}_{{\rm{o}}},{r}_{{\rm{o}}}){{\boldsymbol{A}}}_{{\bf{2}}}{\boldsymbol{I}}(b,\theta ^{\prime} ,{r}_{{\rm{s}}}){{\boldsymbol{A}}}_{{\bf{1}}}{\boldsymbol{R}}^{\prime} ({x}_{{\rm{o}}},{y}_{{\rm{o}}}){\boldsymbol{R}}({\rm{I}},{\rm{\Lambda }},{\rm{\Theta }}){\boldsymbol{y}},\end{eqnarray} \tag{ A18 }$$

where

$$\begin{eqnarray}&&\theta ^{\prime} =\arctan 2({x}_{{\rm{o}}},{y}_{{\rm{o}}})-\arctan 2({x}_{{\rm{s}}},{y}_{{\rm{s}}})\end{eqnarray} \tag{ A19 }$$

is the angle of the terminator in the frame ${ \mathcal F }^{\prime}$. Note that we made use of the fact that A = A ₂ A ₁ (Equation (14) in Luger et al. 2019a), where A ₁ transforms from the spherical harmonic basis $\tilde{{\boldsymbol{y}}}$ to the polynomial basis $\tilde{{\boldsymbol{p}}}$, and A ₂ transforms from $\tilde{{\boldsymbol{p}}}$ to the Green's basis $\tilde{{\boldsymbol{g}}}$.

Similarly, the flux when there is no occultation is now given by

$$\begin{eqnarray}&&{f}_{0}={{\mathbb{r}}}^{\top }(b){\boldsymbol{I}}(b,\theta ^{\prime\prime} ,{r}_{{\rm{s}}}){{\boldsymbol{A}}}_{{\bf{1}}}{\boldsymbol{R}}^{\prime\prime} ({x}_{{\rm{s}}},{y}_{{\rm{s}}}){\boldsymbol{R}}({\rm{I}},{\rm{\Lambda }},{\rm{\Theta }}){\boldsymbol{y}},\end{eqnarray} \tag{ A20 }$$

where

$$\begin{eqnarray}&&\theta ^{\prime\prime} =0\end{eqnarray} \tag{ A21 }$$

is the angle of the terminator in the frame ${ \mathcal F }^{\prime\prime}$, by construction. The transformation R '' = R ''(x_s, y_s) rotates the body through an angle $\arctan 2({x}_{{\rm{s}}},{y}_{{\rm{s}}})$ so the semimajor axis of the terminator is aligned with the x''-axis; as will become clear in Appendix B below, this greatly simplifies the integration step.

Note that in both equations we replaced the integral vectors r ^⊤ and s ^⊤(b_o, r_o) with the vectors ${{\mathbb{r}}}^{\top }(b)$ and ${{\mathbb{s}}}^{\top }(b,\theta ^{\prime} ,{b}_{{\rm{o}}},{r}_{{\rm{o}}})$, respectively. As we mentioned above, we must modify the integration limits to exclude the nightside, where the weighting by I is unphysical. The vectors ${{\mathbb{r}}}^{\top }$ and ${{\mathbb{s}}}^{\top }$ correspond to these modified integrals, which we devote the rest of this paper to computing.

Figure 14 summarizes the transformations involved in the two equations above. Starting on the right with a map vector y in the spherical harmonic basis $\tilde{{\boldsymbol{y}}}$, defined in some observer-independent frame ${{ \mathcal F }}_{0}$, we first rotate it via R to the sky frame ${ \mathcal F }$, in which the body is viewed by the observer. If an occultor is present (upper branch of the figure), we rotate the map from ${ \mathcal F }$ via ${\boldsymbol{R}}^{\prime}$ to the frame ${ \mathcal F }^{\prime}$, in which the occultor lies along the $+y^{\prime}$-axis. We then apply A ₁ to change basis to $\tilde{{\boldsymbol{p}}}$ and I to weight the map by the illumination. Finally, we change basis via A ₂ to the Green's basis, in which we compute and dot the integrals ${{\mathbb{s}}}^{\top }$. If, on the other hand, there is no occultation (lower branch of the figure), we instead rotate the map via R '' to the integration frame ${ \mathcal F }^{\prime\prime}$, in which the terminator is parallel to the x''-axis. We then apply A ₁ to change basis to $\tilde{{\boldsymbol{p}}}$, apply the illumination transform I , and finally dot in the solutions to the surface integrals ${{\mathbb{r}}}^{\top }$.

Zoom In Zoom Out Reset image size

The key to identifying the case corresponding to a given configuration is to determine whether or not the limb of the occultor intersects the terminator of the body, and if so, the points of intersection. While we perform the integration in frame ${ \mathcal F }^{\prime}$, finding the points of intersection with the terminator is easier if we temporarily switch to the frame ${ \mathcal F }^{\prime\prime}$, in which the terminator is parallel to the x''-axis. In this frame, the equations defining the terminator and the limb of the occultor are, respectively,

$$\begin{eqnarray}\begin{array}{rcl}{y}_{1}^{{\prime\prime} }(x^{\prime\prime} ) & = & b\sqrt{1-{x}^{{\prime\prime} 2}}\\ {y}_{2}^{{\prime\prime} }(x^{\prime\prime} ) & = & {y}_{{\rm{o}}}^{{\prime\prime} }\pm \sqrt{{r}_{{\rm{o}}}^{2}-{\left(x^{\prime\prime} -{x}_{{\rm{o}}}^{{\prime\prime} }\right)}^{2}},\end{array}\end{eqnarray} \tag{ C1 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{x}_{{\rm{o}}}^{{\prime\prime} } & = & {b}_{{\rm{o}}}\sin \theta ^{\prime} \\ {y}_{{\rm{o}}}^{{\prime\prime} } & = & {b}_{{\rm{o}}}\cos \theta ^{\prime} \end{array}\end{eqnarray} \tag{ C2 }$$

are the coordinates of the occultor in ${ \mathcal F }^{\prime\prime}$.^, We wish to find the vector of N points ${\boldsymbol{x}}^{\prime\prime} ={\left({x}_{0},{x}_{1},\cdot \cdot \cdot ,{x}_{N-1}\right)}^{\top }$ for which ${y}_{1}^{{\prime\prime} }({x}_{n}^{\prime\prime} )-{y}_{2}^{{\prime\prime} }({x}_{n}^{\prime\prime} )=0$. Following Luger et al. (2017), we may express this condition as the quartic equation

$$\begin{eqnarray}&&{Ax}{{\prime\prime} }^{4}+{Bx}{{\prime\prime} }^{3}+{Cx}{{\prime\prime} }^{2}+{Dx}^{\prime\prime} +E=0,\end{eqnarray} \tag{ C3 }$$

with coefficients

$$\begin{eqnarray}\begin{array}{rcl}A & = & {\left(1-{b}^{2}\right)}^{2}\\ B & = & -4{x}_{{\rm{o}}}^{{\prime\prime} }(1-{b}^{2})\\ C & = & -2\left({b}^{4}+{r}_{{\rm{o}}}^{2}-3{{x}_{{\rm{o}}}^{{\prime\prime} }}^{2}-{{y}_{{\rm{o}}}^{{\prime\prime} }}^{2}-{b}^{2}\left(1+{r}_{{\rm{o}}}^{2}-{{x}_{{\rm{o}}}^{{\prime\prime} }}^{2}+{{y}_{{\rm{o}}}^{{\prime\prime} }}^{2}\right)\right)\\ D & = & -4{x}_{{\rm{o}}}^{{\prime\prime} }({b}^{2}-{r}_{{\rm{o}}}^{2}+{{x}_{{\rm{o}}}^{{\prime\prime} }}^{2}+{{y}_{{\rm{o}}}^{{\prime\prime} }}^{2})\\ E & = & {b}^{4}-2{b}^{2}\left({r}_{{\rm{o}}}^{2}-{{x}_{{\rm{o}}}^{{\prime\prime} }}^{2}+{{y}_{{\rm{o}}}^{{\prime\prime} }}^{2}\right)+{\left({r}_{{\rm{o}}}^{2}-{{x}_{{\rm{o}}}^{{\prime\prime} }}^{2}-{{y}_{{\rm{o}}}^{{\prime\prime} }}^{2}\right)}^{2}.\end{array}\end{eqnarray} \tag{ C4 }$$

Although closed-form solutions to quartic equations exist (see, e.g., Hughes & Chraibi 2011, who solve for the area of overlap between two ellipses analytically), they are prone to significant numerical instabilities. Instead, we solve for the roots of the quartic numerically by casting it as an eigenvalue problem (e.g., Edelman & Murakami 1995) and polish the results with a few iterations of Newton's method. We find that this is reasonably computationally efficient and yields roots with precision within a couple orders of magnitude of machine epsilon (see Section 3.2). ^{47 , 48 , 49 , 50}

In general, the quartic defined by Equation (C3) has N = 4 (potentially degenerate) roots, some of which may be complex, and some of which correspond to intersections with the wrong half of the terminator ellipse (i.e., the section of the terminator on the far side of the body). After excluding the unphysical solutions, we are still left with anywhere between zero and four roots.

Cases with zero roots (case 1–case 5) are treated in Appendix C.2, while cases with one or two roots (case 6–case 10) are treated in Appendix C.3. Cases with three or four roots (case 11–case 14) are rarely encountered in practice, but are possible for some pathological configurations; these are treated in Appendix C.4.

Cases 1–5 (see Figure 15) involve configurations in which the occultor does not intersect with the terminator of the occulted body, and are therefore fairly straightforward to solve. In particular, we can use the original emitted light solution from Luger et al. (2019a), provided we weight the map by our polynomial illumination function:

$$\begin{eqnarray}&&{f}_{{\rm{I}}}={{\boldsymbol{s}}}^{\top }({b}_{{\rm{o}}},{r}_{{\rm{o}}}){{\boldsymbol{A}}}_{{\bf{2}}}{\boldsymbol{I}}(b,\theta ^{\prime} ,{r}_{{\rm{s}}}){{\boldsymbol{A}}}_{{\bf{1}}}{\boldsymbol{R}}^{\prime} ({x}_{{\rm{o}}},{y}_{{\rm{o}}}){\boldsymbol{R}}({\rm{I}},{\rm{\Lambda }},{\rm{\Theta }}){\boldsymbol{y}},\end{eqnarray} \tag{ C5 }$$

where s ^⊤(b_o, r_o) is the emitted light solution vector (Equation (26) in Luger et al. 2019a). The flux f_I is the flux one would measure from a body whose surface map is weighted by the illumination function ${\boldsymbol{I}}(b,\theta ^{\prime} ,{r}_{{\rm{s}}})$ during an occultation. Note that this is not necessarily the observed flux, since this may include the unphysical negative contribution from the nightside. We must compute the actual observed flux on a case-by-case basis.

Case 1 corresponds to any complete occultation of the body (b_o ≤ r_o − 1), so the solution for the flux is trivial:

$$\begin{eqnarray}&&{f}_{1}=0.\end{eqnarray} \tag{ C6 }$$

Case 2 corresponds to occultations in which the occultor blocks all of the dayside of the body and some of the nightside. In this configuration, the unocculted part of the disk consists only of nightside, so the solution is again trivial:

$$\begin{eqnarray}&&{f}_{2}=0.\end{eqnarray} \tag{ C7 }$$

Conversely, case 3 corresponds to occultations in which the occultor blocks all of the nightside of the body and some of the dayside. Since the visible portion of the disk consists only of dayside, we can simply use the weighted solution in emitted light (Equation (C5)):

$$\begin{eqnarray}&&{f}_{3}={f}_{{\rm{I}}}.\end{eqnarray} \tag{ C8 }$$

Case 4 involves any occultation in which the occultor blocks only the nightside of the body (regardless of whether or not it intersects with the limb of the body). Since the nightside intensity is zero everywhere, this case is also trivial, as the flux is equal to the flux in the no occultation case (Equation (B6)):

$$\begin{eqnarray}&&{f}_{4}={f}_{0}.\end{eqnarray} \tag{ C9 }$$

Finally, case 5 involves any occultation in which the occultor blocks only the dayside of the body (regardless of whether or not it intersects with the limb). We first compute the illumination-weighted flux f_I as above, then negate the unphysical nightside contribution using Equation (B7):

$$\begin{eqnarray}&&{f}_{5}={f}_{{\rm{I}}}-{\hat{f}}_{0}.\end{eqnarray} \tag{ C10 }$$

Cases 6–10 (see Figure 15) correspond to configurations in which the limb of the occultor intersects with the terminator at either one point (case 6) or two points (cases 7–10). Because of these intersections, we cannot simply reweight the emitted light solution, as the integration boundaries are now different. In general, we may compute the flux by integrating the components of the Green's basis $\tilde{{\boldsymbol{g}}}$ over the region S bounded by three curves, which we denote ${ \mathcal P }$, ${ \mathcal T }$, and ${ \mathcal Q }$. These are shown in orange in Figure 15 and presented in more detail in Figure 16. Curve ${ \mathcal P }$ is a segment of the limb of the occultor, parametrized by the angle ϕ ∈ [ϕ₀, ϕ₁]; curve ${ \mathcal T }$ is a segment of the terminator, parametrized by the angle ξくしー ∈ [ξくしー₀, ξくしー₁]; and curve ${ \mathcal Q }$ is a segment of the limb of the occulted body, parametrized by the angle λらむだ ∈ [λらむだ₀, λらむだ₁]. The endpoints ϕ₀, ϕ₁, ξくしー₀, ξくしー₁, λらむだ₀, and λらむだ₁ are functions of the solutions to the quartic from Appendix C.1 and will be presented in Appendix C.5.

Zoom In Zoom Out Reset image size

Let ${{\mathbb{s}}}^{\top }$ be the integral of ${\tilde{{\boldsymbol{g}}}}^{\top }$ over S:

$$\begin{eqnarray}&&{{\mathbb{s}}}^{\top }(b,\theta ^{\prime} ,{b}_{{\rm{o}}},{r}_{{\rm{o}}})=\mathop{\iint }\limits_{S(b,\theta ^{\prime} ,{b}_{{\rm{o}}},{r}_{{\rm{o}}})}{\tilde{{\boldsymbol{g}}}}^{\top }(x^{\prime} ,y^{\prime} )\ {dx}^{\prime} \ {dy}^{\prime} .\end{eqnarray} \tag{ C11 }$$

We defer the solution to Equation (C11) to Appendix C.5 below, as it is quite lengthy. Given ${{\mathbb{s}}}^{\top }$, the flux f_S over the integration region S is computed from Equation (A18):

$$\begin{eqnarray}&&{f}_{{\rm{S}}}={{\mathbb{s}}}^{\top }(b,\theta ^{\prime} ,{b}_{{\rm{o}}},{r}_{{\rm{o}}}){{\boldsymbol{A}}}_{{\bf{2}}}{\boldsymbol{I}}(b,\theta ^{\prime} ,{r}_{{\rm{s}}}){{\boldsymbol{A}}}_{{\bf{1}}}{\boldsymbol{R}}^{\prime} ({x}_{{\rm{o}}},{y}_{{\rm{o}}}){\boldsymbol{R}}({\rm{I}},{\rm{\Lambda }},{\rm{\Theta }}){\boldsymbol{y}}.\end{eqnarray} \tag{ C12 }$$

Note again that this is not necessarily the observed flux, which we must compute on a case-by-case basis below.

Case 6 corresponds to configurations in which the limb of the occultor intersects the terminator at a single point. The integration region (see Figures 15 and 16) is the occulted portion of the dayside, which is bounded by all three curves ${ \mathcal P }$, ${ \mathcal Q }$, and ${ \mathcal T }$. The total flux may be computed by subtracting the occulted flux f_S from the total dayside flux f₀:

$$\begin{eqnarray}&&{f}_{6}={f}_{0}-{f}_{{\rm{S}}}.\end{eqnarray} \tag{ C13 }$$

Cases 7–10 involve two points of intersection between the occultor limb and the terminator. Cases 7 and 8 correspond to occultors that block some of the nightside and some of the dayside, but neither of the extrema of the terminator ellipse. In case 7 a lens-shaped region is formed by the intersection of the occultor limb and the terminator on the nightside, while in case 8 this region is formed on the dayside. In case 7, we begin by computing the flux over the unocculted region, f_I, which includes the spurious nightside contribution. We then remove this contribution by noting that it is equal to the total nightside contribution, ${\hat{f}}_{0}$, minus the occulted nightside flux, f_S:

$$\begin{eqnarray}&&{f}_{7}={f}_{{\rm{I}}}-({\hat{f}}_{0}-{f}_{{\rm{S}}}).\end{eqnarray} \tag{ C14 }$$

Case 8, on the other hand, is equivalent to case 6, since the integration region consists of occulted dayside:

$$\begin{eqnarray}&&{f}_{8}={f}_{0}-{f}_{{\rm{S}}}.\end{eqnarray} \tag{ C15 }$$

Cases 9 and 10 correspond to occultors that also block some nightside and some dayside, along with both of the extrema of the ellipse; these are therefore exclusively for large occultors (r_o > 1). Case 9 involves occultations in which only a small lens-shaped region of the nightside is visible. The total flux is the visible dayside plus unphysical nightside contribution, f_I, minus the nightside contribution, which we compute from Equation (C12):

$$\begin{eqnarray}&&{f}_{9}={f}_{{\rm{I}}}-{f}_{{\rm{S}}}.\end{eqnarray} \tag{ C16 }$$

Conversely, case 10 involves occultations in which only a small lens-shaped region of the dayside is visible. In this case, we may compute the observed flux from Equation (C12) directly:

$$\begin{eqnarray}&&{f}_{10}={f}_{{\rm{S}}}.\end{eqnarray} \tag{ C17 }$$

Cases 11–14 correspond to (rare) configurations involving three or four roots to Equation (C3) and are illustrated in Figure 17. All four involve integration over two disjoint regions (see the figure). Cases 11 and 12 involve three points of intersection between the terminator and the occultor limb. In case 11, the regions of integration S₁ and S₂ are the occulted portion of the dayside, so the solution is similar to that of cases 6 and 8:

$$\begin{eqnarray}&&{f}_{11}={f}_{0}-({f}_{{{\rm{S}}}_{1}}+{f}_{{{\rm{S}}}_{2}}),\end{eqnarray} \tag{ C18 }$$

where ${f}_{{{\rm{S}}}_{1}}$ and ${f}_{{{\rm{S}}}_{2}}$ are computed from Equation (C12) for each of the integration regions. Conversely, in case 12 the two regions are the occulted portion of the nightside, so the solution is similar to that of case 7:

$$\begin{eqnarray}&&{f}_{12}={f}_{{\rm{I}}}-\left({\hat{f}}_{0}-({f}_{{{\rm{S}}}_{1}}+{f}_{{{\rm{S}}}_{2}})\right).\end{eqnarray} \tag{ C19 }$$

Finally, cases 13 and 14 involve four points of intersection between the terminator and the occultor limb. The regions of integration in case 13 are the visible portion of the nightside, so this case is equivalent to case 9:

$$\begin{eqnarray}&&{f}_{13}={f}_{{\rm{I}}}-({f}_{{{\rm{S}}}_{1}}+{f}_{{{\rm{S}}}_{2}}).\end{eqnarray} \tag{ C20 }$$

Conversely, the regions of integration in case 14 are the visible portion of the dayside, so this case is equivalent to case 10:

$$\begin{eqnarray}&&{f}_{14}={f}_{{{\rm{S}}}_{1}}+{f}_{{{\rm{S}}}_{2}}.\end{eqnarray} \tag{ C21 }$$

Zoom In Zoom Out Reset image size

In the previous sections, we discussed how to identify the case corresponding to a specific configuration of the occultor and the illumination source. We showed how in some cases (1–5; Appendix C.2) the total flux may be computed by exploiting the classical starry integrals (Equation (C5)). In all other cases (6–14; Appendices C.3 and C.4), however, the flux computation involves evaluation of Equation (C12), where the solution vector ${{\mathbb{s}}}^{\top }$ is the vector of integrals (in the frame ${ \mathcal F }^{\prime}$) of each of the terms in the Green's basis $\tilde{{\boldsymbol{g}}}$ over a region S of the projected disk of the occultor (Equation (C11)). As in Luger et al. (2019a), the approach to computing ${{\mathbb{s}}}^{\top }$ is to use Green's theorem to transform the surface integrals into line integrals along the curves ${ \mathcal P }$, ${ \mathcal T }$, and ${ \mathcal Q }$ (see Figure 16). Specifically, we write

$$\begin{eqnarray}\begin{array}{rcl}{{\mathbb{s}}}^{\top } & = & \mathop{\iint }\limits_{S}{\tilde{{\boldsymbol{g}}}}^{\top }(x^{\prime} ,y^{\prime} )\ {dx}^{\prime} \ {dy}^{\prime} \\ & = & \oint {{\boldsymbol{G}}}^{\top }(x^{\prime} ,y^{\prime} )\cdot d{\boldsymbol{r}}(x^{\prime} ,y^{\prime} ),\end{array}\end{eqnarray} \tag{ C22 }$$

where ${\boldsymbol{G}}(x^{\prime} ,y^{\prime} )$ is a vector of two-dimensional Cartesian vectors chosen such that its exterior derivative is $\tilde{{\boldsymbol{g}}}$,

$$\begin{eqnarray}&&\displaystyle \frac{d{{\boldsymbol{G}}}_{y^{\prime} }(x^{\prime} ,y^{\prime} )}{{dx}^{\prime} }-\displaystyle \frac{d{{\boldsymbol{G}}}_{x^{\prime} }(x^{\prime} ,y^{\prime} )}{{dy}^{\prime} }=\tilde{{\boldsymbol{g}}}(x^{\prime} ,y^{\prime} ),\end{eqnarray} \tag{ C23 }$$

and

$$\begin{eqnarray}&&d{\boldsymbol{r}}(x^{\prime} ,y^{\prime} )=\left(\displaystyle \frac{{dx}^{\prime} }{d\varphi }\right)d\varphi \,\hat{{\boldsymbol{x}}}^{\prime} +\left(\displaystyle \frac{{dy}^{\prime} }{d\varphi }\right)d\varphi \,\hat{{\boldsymbol{y}}}^{\prime} ,\end{eqnarray} \tag{ C24 }$$

where φふぁい is the parametrized angle along the integration path and the integral is taken in a counter-clockwise direction relative to the center of the integration region. Luger et al. (2019a) showed that one possible solution to Equation (C23) consists of the vector whose nth component is given by

$$\begin{eqnarray}&&\begin{array}{l}{{\boldsymbol{G}}}_{n}(x^{\prime} ,y^{\prime} )\,=\\ \left\{\begin{array}{ll}x{{\prime} }^{\tfrac{\mu +2}{2}}y{{\prime} }^{\tfrac{\nu }{2}}\,\hat{{\boldsymbol{y}}}^{\prime} & \mu ,\nu \,\mathrm{even}\\ \displaystyle \frac{1-z{\left(x^{\prime} ,y^{\prime} \right)}^{3}}{3{\left(1-z(x^{\prime} ,y^{\prime} \right)}^{2})}\left(-y^{\prime} \,\hat{{\boldsymbol{x}}}^{\prime} +x^{\prime} \,\hat{{\boldsymbol{y}}}^{\prime} \right) & \qquad \mu =\nu =1\\ x{{\prime} }^{l-2}z{\left(x^{\prime} ,y^{\prime} \right)}^{3}\,\hat{{\boldsymbol{x}}} & \nu \,\mathrm{odd},\mu =1,l\,\mathrm{even}\\ x{{\prime} }^{l-3}y^{\prime} z{\left(x^{\prime} ,y^{\prime} \right)}^{3}\,\hat{{\boldsymbol{x}}} & \nu \,\mathrm{odd},\mu =1,l\,\mathrm{odd}\\ x{{\prime} }^{\tfrac{\mu -3}{2}}y{{\prime} }^{\tfrac{\nu -1}{2}}z{\left(x^{\prime} ,y^{\prime} \right)}^{3}\,\hat{{\boldsymbol{y}}} & \mathrm{otherwise},\end{array}\right.,\end{array}\end{eqnarray} \tag{ C25 }$$

where the indices l, m, μみゅー, νにゅー are given by Equations (A2) and (A5). ⁵³

We showed in the previous sections that there are at most three curves ${ \mathcal P }$, ${ \mathcal T }$, and ${ \mathcal Q }$ bounding a given closed surface of integration (see Figure 16). We may therefore express Equation (C22) as

$$\begin{eqnarray}&&{{\mathbb{s}}}^{\top }={{\mathbb{p}}}^{\top }+{{\mathbb{t}}}^{\top }+{{\mathbb{q}}}^{\top },\end{eqnarray} \tag{ C26 }$$

where we define the primitive integrals

$$\begin{eqnarray}&&{{\mathbb{p}}}^{\top }={\int }_{{\boldsymbol{\phi }}}{{\boldsymbol{G}}}^{\top }({x}_{{\rm{p}}}^{{\prime} },{y}_{{\rm{p}}}^{{\prime} })\cdot d{\boldsymbol{r}}({x}_{{\rm{p}}}^{{\prime} },{y}_{{\rm{p}}}^{{\prime} }),\end{eqnarray} \tag{ C27 }$$

$$\begin{eqnarray}&&{{\mathbb{t}}}^{\top }={\int }_{{\boldsymbol{\xi }}}{{\boldsymbol{G}}}^{\top }({x}_{{\rm{t}}}^{{\prime} },{y}_{{\rm{t}}}^{{\prime} })\cdot d{\boldsymbol{r}}({x}_{{\rm{t}}}^{{\prime} },{y}_{{\rm{t}}}^{{\prime} }),\end{eqnarray} \tag{ C28 }$$

$$\begin{eqnarray}&&{{\mathbb{q}}}^{\top }={\int }_{{\boldsymbol{\lambda }}}{{\boldsymbol{G}}}^{\top }({x}_{{\rm{q}}}^{{\prime} },{y}_{{\rm{q}}}^{{\prime} })\cdot d{\boldsymbol{r}}({x}_{{\rm{q}}}^{{\prime} },{y}_{{\rm{q}}}^{{\prime} })\end{eqnarray} \tag{ C29 }$$

to be the line integrals of G along each of the curves ${ \mathcal P }$, ${ \mathcal T }$, and ${ \mathcal Q }$, respectively,^{, ,} where the coordinates along each curve are parametrized in terms of φふぁい as follows:

$$\begin{eqnarray}&&\begin{array}{l}{x}_{{\rm{p}}}^{{\prime} }={r}_{{\rm{o}}}\cos \varphi \\ {x}_{{\rm{t}}}^{{\prime} }=\cos \theta ^{\prime} \cos \varphi -b\sin \theta ^{\prime} \sin \varphi \\ {x}_{{\rm{q}}}^{{\prime} }=\cos \varphi \\ {y}_{{\rm{p}}}^{{\prime} }={b}_{{\rm{o}}}+{r}_{{\rm{o}}}\sin \varphi \\ {y}_{{\rm{t}}}^{{\prime} }=\sin \theta ^{\prime} \cos \varphi +b\cos \theta ^{\prime} \sin \varphi \\ {y}_{{\rm{q}}}^{{\prime} }=\sin \varphi ,\end{array}\end{eqnarray} \tag{ C30 }$$

and we define

$$\begin{eqnarray}\displaystyle \begin{array}{rcl}{\displaystyle \int }_{{\boldsymbol{\varphi }}} & \equiv & {\displaystyle \int }_{{\varphi }_{0}}^{{\varphi }_{1}}+{\displaystyle \int }_{{\varphi }_{2}}^{{\varphi }_{3}}+\cdots +{\displaystyle \int }_{{\varphi }_{N-2}}^{{\varphi }_{N-1}}\\ & \equiv & \sum _{i=0}^{\tfrac{N}{2}-1}{\displaystyle \int }_{{\varphi }_{2i}}^{{\varphi }_{2i+1}}\end{array}\end{eqnarray} \tag{ C31 }$$

to be the sum of definite integrals between pairs of limits φふぁい_i arranged in a vector φふぁい of length N. For future reference, it will also be useful to define the operator

$$\begin{eqnarray}&&{\rm{\Delta }}{\boldsymbol{x}}\equiv \sum _{i=0}^{\tfrac{N}{2}-1}\left({x}_{2i+1}-{x}_{2i}\right),\end{eqnarray} \tag{ C32 }$$

which sums the difference of successive pairs of values in a vector ${\boldsymbol{x}}={\left({x}_{0},{x}_{1},{x}_{2},{x}_{3},\cdot \cdot \cdot ,\,{x}_{N-1}\right)}^{\top }$. This will come in handy when computing definite integrals. Specifically, if g is the antiderivative of some function f, we may use the fundamental theorem of calculus to compute the integral of f over the interval(s) given by the vector of limit pairs φふぁい :

$$\begin{eqnarray}\begin{array}{rcl}{\int }_{{\boldsymbol{\varphi }}}f(\varphi )d\varphi & = & {\displaystyle \int }_{{\varphi }_{0}}^{{\varphi }_{1}}f(\varphi )d\varphi +{\displaystyle \int }_{{\varphi }_{2}}^{{\varphi }_{3}}f(\varphi )d\varphi +\cdot \cdot \cdot \\ & & +\,{\displaystyle \int }_{{\varphi }_{N-2}}^{{\varphi }_{N-1}}f(\varphi )d\varphi \\ & = & {\rm{\Delta }}{\boldsymbol{g}}({\boldsymbol{\varphi }}),\end{array}\end{eqnarray} \tag{ C33 }$$

where g is the vector given by

$$\begin{eqnarray}&&\begin{array}{l}{\boldsymbol{g}}({\boldsymbol{\varphi }})=\left(g({\varphi }_{0}),g({\varphi }_{1}),g({\varphi }_{2}),g({\varphi }_{3}),\cdot \cdot \cdot ,\right.\\ \ \ {\left.g({\varphi }_{N-2}),g({\varphi }_{N-1})\right)}^{\top }.\end{array}\end{eqnarray} \tag{ C34 }$$

Note that most of the cases (1–10) involve integration over a single closed region, so Equations (C31) and (C32) reduce to

$$\begin{eqnarray}&&{\int }_{{\boldsymbol{\varphi }}}\equiv {\int }_{{\varphi }_{0}}^{{\varphi }_{1}}\end{eqnarray} \tag{ C35 }$$

and

$$\begin{eqnarray}&&{\rm{\Delta }}{\boldsymbol{x}}\equiv {x}_{1}-{x}_{0}.\end{eqnarray} \tag{ C36 }$$

For cases 11–14, we must integrate over two disjoint regions, so we sum over two pairs of limits. In the next three sections, we derive the solutions to each of the primitive integrals ${{\mathbb{p}}}^{\top }$, ${{\mathbb{t}}}^{\top }$, and ${{\mathbb{q}}}^{\top }$. ^{54 , 55 , 56 , 57 , 58 , 59}

In this section we present a solution to Equation (C27). The first order of business is to derive expressions for the integration limits ϕ . Depending on the integration case, these limits will correspond to the point of intersection between the limb of the occultor and the limb of the occulted body and/or the point of intersection between the limb of the occultor and the terminator of the occulted body. The former is given by (see Equation (24) in Luger et al. 2019a)

$$\begin{eqnarray}&&{\phi }_{0}=\displaystyle \frac{\pi }{2}\pm \left(\arcsin \left(\displaystyle \frac{1-{r}_{{\rm{o}}}^{2}-{b}_{{\rm{o}}}^{2}}{2{b}_{{\rm{o}}}{r}_{{\rm{o}}}}\right)-\displaystyle \frac{\pi }{2}\right),\end{eqnarray} \tag{ C37 }$$

where the sign is chosen such that the point $({r}_{{\rm{o}}}\cos {\phi }_{0},{b}_{{\rm{o}}}+{r}_{{\rm{o}}}\sin {\phi }_{0})$ is on the dayside of the occulted body, and the latter (of which there may be multiple) is given by

$$\begin{eqnarray}&&{{\boldsymbol{\phi }}}_{{\bf{1}}}=\theta ^{\prime} +\arctan 2\left(b\sqrt{1-{\boldsymbol{x}}{{\prime\prime} }^{2}}-{y}_{{\rm{o}}}^{{\prime\prime} },{\boldsymbol{x}}^{\prime\prime} -{x}_{{\rm{o}}}^{{\prime\prime} }\right),\end{eqnarray} \tag{ C38 }$$

where x '' are the roots of the quartic (Equation (C3)). These angles are then wrapped to the range [0, 2πぱい) and sorted into the vector ϕ such that the integration is always performed in a counter-clockwise sense about the center of the integration region. The left panel in Figure 16 shows a configuration in which the lower integration limit ϕ₀ = 1582 corresponds to the point of intersection between the limbs of the two bodies and the upper integration limit ϕ₁ = 2233 corresponds to the limb-terminator intersection. Both angles are measured counter-clockwise from the line $x^{\prime} ={b}_{{\rm{o}}}$. ^{60 , 61}

In order to evaluate the integral in Equation (C27), we follow the reparametrization tricks of Appendix D.2.3 in Luger et al. (2019a). The algebra is long and tedious, so we merely present the result (alongside the usual validation links). The nth component of ${{\mathbb{p}}}^{\top }$ is

$$\begin{eqnarray}{{\mathbb{p}}}_{n}=\left\{\begin{array}{ll}2{\left(2{r}_{{\rm{o}}}\right)}^{l+2}\left\{\begin{array}{ll}V\left(\displaystyle \frac{\mu +4}{4},\displaystyle \frac{\nu }{2},0;{\mathbb{i}}\right) & \tfrac{\mu }{2}\,\mathrm{even}\\ V\left(\displaystyle \frac{\mu +2}{4},\displaystyle \frac{\nu }{2},0;{\mathbb{u}}\right) & \displaystyle \frac{\mu }{2}\,\mathrm{odd}\end{array}\right. & \mu ,\nu \,\mathrm{even}\\ {{\mathbb{p}}}_{2} & \mu =\nu =1\\ \beta {\left(2{r}_{{\rm{o}}}\right)}^{l-1}\left\{\begin{array}{ll}\left\{\begin{array}{ll}V\left(\displaystyle \frac{l-2}{2},0,0;{\mathbb{j}}\right)-2V\left(\displaystyle \frac{l-2}{2},0,1;{\mathbb{j}}\right) & l\,\mathrm{even}\\ V\left(\displaystyle \frac{l-3}{2},1,0;{\mathbb{j}}\right)-2V\left(\displaystyle \frac{l-3}{2},1,1;{\mathbb{j}}\right) & l\,\mathrm{odd}\ne 1\end{array}\right. & \mu =1\\ \left\{\begin{array}{ll}2V\left(\displaystyle \frac{\mu -1}{4},\tfrac{\nu -1}{2},0;{\mathbb{j}}\right) & \tfrac{\mu -1}{2}\,\mathrm{even}\\ 2V\left(\tfrac{\mu -1}{4},\tfrac{\nu -1}{2},0;{\mathbb{w}}\right) & \tfrac{\mu -1}{2}\,\mathrm{odd}\end{array}\right. & \mu \gt 1\end{array}\right. & \mu ,\nu \,\mathrm{odd}\end{array}\right.,\end{eqnarray} \tag{ C39 }$$

where $\beta ={\left(1-{\left({b}_{{\rm{o}}}-{r}_{{\rm{o}}}\right)}^{2}\right)}^{\tfrac{3}{2}}$ and we define the Vieta operator

$$\begin{eqnarray}&&V\left(u,v,w;{\boldsymbol{x}}\right)\equiv \sum _{i=0}^{u+v}{{ \mathcal A }}_{u,v,i}{x}_{u+w+i}\end{eqnarray} \tag{ C40 }$$

as the dot product of a vector x and the vector of Vieta's theorem coefficients, where (see Equation (D34) in Luger et al. 2019a)

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal A }}_{u,v,i} & = & \displaystyle \sum _{j=\max (0,u-i)}^{\min (u+v-i,u)}\left(\displaystyle \genfrac{}{}{0em}{}{u}{j}\right)\left(\displaystyle \genfrac{}{}{0em}{}{v}{u+v-i-j}\right)\\ & & \times {\left(-1\right)}^{u+j}{\left(\displaystyle \frac{{b}_{{\rm{o}}}-{r}_{{\rm{o}}}}{2{r}_{{\rm{o}}}}\right)}^{u+v-i-j}.\end{array}\end{eqnarray} \tag{ C41 }$$

The vectors ${\mathbb{i}}$, ${\mathbb{j}}$, ${\mathbb{u}}$, and ${\mathbb{w}}$ are solutions to specific integrals, which we compute recursively below. As in Luger et al. (2019a) the n = 2 term of ${{\mathbb{p}}}^{\top }$, ${{\mathbb{p}}}_{2}$, is handled separately; we also compute this below. ^{62 , 63 , 64}

Note that several of the cases in Equation (C39) are identical to those in Equation (D35) of Luger et al. (2019a), provided we replace their integrals ${ \mathcal I }$ and ${ \mathcal J }$ with our integrals ${\mathbb{i}}$ and ${\mathbb{j}}$, respectively. The integrals themselves are similar, except for a change in the limits of integration, which are no longer symmetric about zero. As we will see, this leads to the dependence of these expressions on incomplete elliptic integrals. Note also that the integrals ${\mathbb{u}}$ and ${\mathbb{w}}$ are new, as certain cancellations in Luger et al. (2019a) resulted in the corresponding cases contributing zero net flux (last case in Equation (D35) of Luger et al. 2019a).

C.6.1. The Vector ${\mathbb{i}}$

The components of the vector ${\mathbb{i}}$ are given by the integral

$$\begin{eqnarray}&&{{\mathbb{i}}}_{v}({\boldsymbol{\alpha }})={\int }_{{\boldsymbol{\alpha }}}{\sin }^{2v}\varphi d\varphi ,\end{eqnarray} \tag{ C42 }$$

for $v\in [0,{v}_{\max }]$, where we define the helper angle

$$\begin{eqnarray}&&{\boldsymbol{\alpha }}\equiv \displaystyle \frac{{\boldsymbol{\phi }}}{2}+\displaystyle \frac{\pi }{4}.\end{eqnarray} \tag{ C43 }$$

The integral in the expression above is the same as that in Equation (D38) of Luger et al. (2019a), except for a change in the limits of integration. As in Luger et al. (2019a), we can compute the vector ${\mathbb{i}}$ recursively given a trivial lower boundary condition:

$$\begin{eqnarray}\begin{array}{rcl}{{\mathbb{i}}}_{0}({\boldsymbol{\alpha }}) & = & {\rm{\Delta }}{\boldsymbol{\alpha }}\\ {{\mathbb{i}}}_{v}({\boldsymbol{\alpha }}) & = & \displaystyle \frac{1}{2v}\left((2v-1){{\mathbb{i}}}_{v-1}({\boldsymbol{\alpha }})-{\rm{\Delta }}\left({\sin }^{2v-1}{\boldsymbol{\alpha }}{\cos }^{2v-1}{\boldsymbol{\alpha }}\right)\right),\end{array}\end{eqnarray} \tag{ C44 }$$

where the last expression is valid for all v > 0. We find that this algorithm is generally stable, except when $\sin {\boldsymbol{\alpha }}$ is small. In that limit, we evaluate ${{\mathbb{i}}}_{{v}_{\max }}({\boldsymbol{\alpha }})$ by numerical integration of Equation (C42) using Gauss–Legendre quadrature with 100 points. We then recurse downward by substituting v → v + 1 in Equation (C44) and solving for ${{\mathbb{i}}}_{v}({\boldsymbol{\alpha }})$. ⁶⁵

C.6.2. The Vector ${\mathbb{j}}$

The components of the vector ${\mathbb{j}}$ are given by the integral

$$\begin{eqnarray}&&{{\mathbb{j}}}_{v}({k}^{2},{\boldsymbol{\alpha }})={\int }_{{\boldsymbol{\alpha }}}{\sin }^{2v}\varphi {\left(1-\displaystyle \frac{{\sin }^{2}\varphi }{{k}^{2}}\right)}^{\displaystyle \frac{3}{2}}d\varphi ,\end{eqnarray} \tag{ C45 }$$

where

$$\begin{eqnarray}&&{k}^{2}\equiv \displaystyle \frac{1-{r}_{{\rm{o}}}^{2}-{b}_{{\rm{o}}}^{2}+2{b}_{{\rm{o}}}{r}_{{\rm{o}}}}{4{b}_{{\rm{o}}}{r}_{{\rm{o}}}}.\end{eqnarray} \tag{ C46 }$$

The integral in this expression is again the same as that in Equation (D39) of Luger et al. (2019a), except for a change in the limits of integration. In that paper, we computed all terms $\{{{\mathbb{j}}}_{0},\cdot \cdot \cdot ,{{\mathbb{j}}}_{{v}_{\max }}\}$ from a three-term recurrence relation and two boundary conditions. In the case of upward recursion, the boundary conditions ${{\mathbb{j}}}_{0}$ and ${{\mathbb{j}}}_{1}$ were computed analytically from the complete elliptic integrals K(k²) and E(k²). In cases where upward recursion was not numerically stable, we evaluated ${{\mathbb{j}}}_{{v}_{\max }}$ and ${{\mathbb{j}}}_{{v}_{\max }-1}$ via a quickly convergent series expansion and recursed downward.

In order to solve Equation (C45), it is possible to replace the complete elliptic integrals K(k²) and E(k²) in the lower boundary conditions (Equation (D46) in Luger et al. 2019a) with the incomplete elliptic integrals

$$\begin{eqnarray}&&F(\psi \,| \,m)\equiv {\int }_{0}^{\psi }\displaystyle \frac{d\varphi }{\sqrt{1-m{\sin }^{2}\varphi }}\end{eqnarray} \tag{ C47 }$$

and

$$\begin{eqnarray}&&E(\psi | m)\equiv {\int }_{0}^{\psi }\sqrt{1-m{\sin }^{2}\varphi }\,d\varphi ,\end{eqnarray} \tag{ C48 }$$

which we compute from the el2 parameterization of Bulirsch (1965), then use the same upward recursion relation to obtain analytic solutions for all ${{\mathbb{j}}}_{v}$:

$$\begin{eqnarray}\begin{array}{rcl}{{\mathbb{j}}}_{0}({k}^{2},{\boldsymbol{\alpha }}) & = & \displaystyle \frac{1}{3}\left(2\left(2-\displaystyle \frac{1}{{k}^{2}}\right){\rm{\Delta }}{\boldsymbol{E}}\left({\boldsymbol{\alpha }}\left|\,\displaystyle \frac{1}{{k}^{2}}\right.\right)\right.\\ & & +\left.\left(\displaystyle \frac{1}{{k}^{2}}-1\right){\rm{\Delta }}{\boldsymbol{F}}\left({\boldsymbol{\alpha }}\left|\,\displaystyle \frac{1}{{k}^{2}}\right.\right)+{\rm{\Delta }}{{\boldsymbol{z}}}_{0}({k}^{2},{\boldsymbol{\alpha }})\right)\\ {{\mathbb{j}}}_{1}({k}^{2},{\boldsymbol{\alpha }}) & = & \displaystyle \frac{1}{15}\left(\left(-3{k}^{2}+13-\displaystyle \frac{8}{{k}^{2}}\right){\rm{\Delta }}{\boldsymbol{E}}\left({\boldsymbol{\alpha }}\left|\,\displaystyle \frac{1}{{k}^{2}}\right.\right)\right.\\ & & \left.+\left(3{k}^{2}-7+\displaystyle \frac{4}{{k}^{2}}\right){\rm{\Delta }}{\boldsymbol{F}}\left({\boldsymbol{\alpha }}\left|\,\displaystyle \frac{1}{{k}^{2}}\right.\right)+{\rm{\Delta }}{{\boldsymbol{z}}}_{1}({k}^{2},{\boldsymbol{\alpha }})\right)\\ {{\mathbb{j}}}_{v}({k}^{2},{\boldsymbol{\alpha }}) & = & \displaystyle \frac{1}{2v+3}\left(2\left(v+1+(v-1){k}^{2}\right){{\mathbb{j}}}_{v-1}({k}^{2},{\boldsymbol{\alpha }})\right.\\ & & \left.-(2v-3){k}^{2}{{\mathbb{j}}}_{v-2}({k}^{2},{\boldsymbol{\alpha }})+{\rm{\Delta }}{{\boldsymbol{z}}}_{v}({k}^{2},{\boldsymbol{\alpha }})\right),\end{array}\end{eqnarray} \tag{ C49 }$$

where the last expression is valid for all v > 1 and

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{z}}}_{0}({k}^{2},{\boldsymbol{\alpha }}) & = & \displaystyle \frac{\sin {\boldsymbol{\alpha }}\cos {\boldsymbol{\alpha }}\,{\boldsymbol{q}}({k}^{2},{\boldsymbol{\alpha }})}{{k}^{2}}\\ {{\boldsymbol{z}}}_{1}({k}^{2},{\boldsymbol{\alpha }}) & = & \left(3{\sin }^{2}{\boldsymbol{\alpha }}+4-6{k}^{2}\right){{\boldsymbol{z}}}_{0}({k}^{2},{\boldsymbol{\alpha }})\\ {{\boldsymbol{z}}}_{v}({k}^{2},{\boldsymbol{\alpha }}) & = & {k}^{2}{\sin }^{2v-3}{\boldsymbol{\alpha }}\cos {\boldsymbol{\alpha }}\,{\boldsymbol{q}}{\left({k}^{2},{\boldsymbol{\alpha }}\right)}^{5},\end{array}\end{eqnarray} \tag{ C50 }$$

with

$$\begin{eqnarray}&&{\boldsymbol{q}}({k}^{2},{\boldsymbol{\alpha }})=\sqrt{1-\displaystyle \frac{{\sin }^{2}{\boldsymbol{\alpha }}}{{k}^{2}}}.\end{eqnarray} \tag{ C51 }$$

Note that when k² < 1 we use the reciprocal-modulus transformation to evaluate the elliptic integrals:

$$\begin{eqnarray}\begin{array}{rcl}F\left(\psi \left|\,\displaystyle \frac{1}{{k}^{2}}\right.\right) & = & k\,F(\beta \,| \,{k}^{2})\\ E\left(\psi \left|\,\displaystyle \frac{1}{{k}^{2}}\right.\right) & = & \displaystyle \frac{E(\beta \,| \,{k}^{2})-(1-{k}^{2})F(\beta \,| \,{k}^{2})}{k}\end{array}\end{eqnarray} \tag{ C52 }$$

with ^{66 , 67}

$$\begin{eqnarray}&&\beta =\arcsin \left(\displaystyle \frac{\sin \psi }{k}\right).\end{eqnarray} \tag{ C53 }$$

In practice, however, we find that this procedure is even more numerically unstable than it was in Luger et al. (2019a). To address this, we express the recurrence structure of the problem as a tridiagonal system with one lower boundary condition ${{\mathbb{j}}}_{0}$ and one upper boundary condition ${{\mathbb{j}}}_{{v}_{\max }}$:

$$\begin{eqnarray}\left(\begin{array}{cccccc}{a}_{0} & 1 & & & & \\ {b}_{1} & {a}_{1} & 1 & & & \\ & {b}_{2} & {a}_{2} & 1 & & \\ & & {b}_{0} & {a}_{3} & 1 & \\ & & & \ddots & \ddots & \ddots \\ & & & & {b}_{{v}_{\max }} & {a}_{{v}_{\max }}\end{array}\right)\left(\begin{array}{c}{{\mathbb{j}}}_{1}\\ {{\mathbb{j}}}_{2}\\ {{\mathbb{j}}}_{3}\\ {{\mathbb{j}}}_{4}\\ \cdot \cdot \cdot \\ {{\mathbb{j}}}_{{v}_{\max }-1}\end{array}\right)=\left(\begin{array}{c}{c}_{0}-{b}_{0}{{\mathbb{j}}}_{0}\\ {c}_{1}\\ {c}_{2}\\ {c}_{3}\\ \cdot \cdot \cdot \\ {c}_{{v}_{\max }}-{{\mathbb{j}}}_{{v}_{\max }}\end{array}\right),\end{eqnarray} \tag{ C54 }$$

where the recursion coefficients are given by

$$\begin{eqnarray}\begin{array}{rcl}{a}_{v}(k) & = & -2\displaystyle \frac{(v+1)+(v-1){k}^{2}}{2v+3}\\ {b}_{v}(k) & = & \displaystyle \frac{(2v-3){k}^{2}}{2v+3}\\ {c}_{v}({k}^{2},{\boldsymbol{\alpha }}) & = & {\rm{\Delta }}\left(\displaystyle \frac{{{\boldsymbol{z}}}_{v}({k}^{2},{\boldsymbol{\alpha }})}{2v+3}\right).\end{array}\end{eqnarray} \tag{ C55 }$$

Solving this matrix system yields values for all intermediate $\{{{\mathbb{j}}}_{1},\cdot \cdot \cdot ,{{\mathbb{j}}}_{{v}_{\max }-1}\}$. While efficient algorithms exist for solving tridiagonal problems, we obtain far better numerical stability by instead performing traditional LU decomposition. We find that this algorithm is stable in all the regimes that we tested. ^{68 , 69}

We evaluate the upper boundary condition ${{\mathbb{j}}}_{{v}_{\max }}$ by numerical integration of Equation (C45) via Gauss–Legendre quadrature with 100 points. While the lower boundary condition may be computed analytically from Equation (C49), in practice we achieve better precision via numerical integration (as above), with negligible effects on computational performance.

C.6.3. The Vector ${\mathbb{u}}$

The components of the vector ${\mathbb{u}}$ are given by the integral

$$\begin{eqnarray}&&{{\mathbb{u}}}_{v}({\boldsymbol{\alpha }})={\int }_{{\boldsymbol{\alpha }}}\cos \varphi {\sin }^{2v+1}\varphi d\varphi .\end{eqnarray} \tag{ C56 }$$

This integral has an analytic solution for all v: ⁷⁰

$$\begin{eqnarray}&&{{\mathbb{u}}}_{v}({\boldsymbol{\alpha }})=\displaystyle \frac{{\rm{\Delta }}{\sin }^{2v+2}{\boldsymbol{\alpha }}}{2v+2}.\end{eqnarray} \tag{ C57 }$$

C.6.4. The Vector ${\mathbb{w}}$

The components of the vector ${\mathbb{w}}$ are given by the integral

$$\begin{eqnarray}&&{{\mathbb{w}}}_{v}({k}^{2},{\boldsymbol{\alpha }})={\int }_{{\boldsymbol{\alpha }}}\cos \varphi {\sin }^{2v+1}\varphi {\left(1-\displaystyle \frac{{\sin }^{2}\varphi }{{k}^{2}}\right)}^{\displaystyle \frac{3}{2}}d\varphi .\end{eqnarray} \tag{ C58 }$$

We may compute it by either upward or downward recursion. In both cases, we compute each of the ${{\mathbb{w}}}_{v}$ from

$$\begin{eqnarray}&&{{\mathbb{w}}}_{v}({k}^{2},{\boldsymbol{\alpha }})={\rm{\Delta }}{{\boldsymbol{b}}}_{v}({k}^{2},{\boldsymbol{\alpha }}).\end{eqnarray} \tag{ C59 }$$

In the upward case, we start with the lower boundary conditions

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{b}}}_{0}({k}^{2},{\boldsymbol{\alpha }}) & = & \displaystyle \frac{{\sin }^{2}{\boldsymbol{\alpha }}}{5}\left(\displaystyle \frac{1-{\boldsymbol{q}}{\left({k}^{2},{\boldsymbol{\alpha }}\right)}^{3}}{1-{\boldsymbol{q}}{\left({k}^{2},{\boldsymbol{\alpha }}\right)}^{2}}+{\boldsymbol{q}}{\left({k}^{2},{\boldsymbol{\alpha }}\right)}^{3}\right)\\ {{\boldsymbol{c}}}_{0}({k}^{2},{\boldsymbol{\alpha }}) & = & {\sin }^{4}{\boldsymbol{\alpha }}\,\displaystyle \frac{{\boldsymbol{q}}{\left({k}^{2},{\boldsymbol{\alpha }}\right)}^{5}}{1-{\boldsymbol{q}}{\left({k}^{2},{\boldsymbol{\alpha }}\right)}^{2}},\end{array}\end{eqnarray} \tag{ C60 }$$

and recurse upward in b and c simultaneously:

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{b}}}_{v}({k}^{2},{\boldsymbol{\alpha }}) & = & \displaystyle \frac{1}{2v+5}\left(\displaystyle \frac{2v{\sin }^{2}{\boldsymbol{\alpha }}}{1-{\boldsymbol{q}}{\left({k}^{2},{\boldsymbol{\alpha }}\right)}^{2}}\right.\\ & & \left.\times \,{{\boldsymbol{b}}}_{v-1}({k}^{2},{\boldsymbol{\alpha }})-{{\boldsymbol{c}}}_{v-1}({k}^{2},{\boldsymbol{\alpha }})\right)\\ {{\boldsymbol{c}}}_{v}({k}^{2},{\boldsymbol{\alpha }}) & = & {\sin }^{2}{\boldsymbol{\alpha }}\,{{\boldsymbol{c}}}_{v-1}({k}^{2},{\boldsymbol{\alpha }})\end{array}\end{eqnarray} \tag{ C61 }$$

for v > 0. In the case of downward recursion, we start with the upper boundary conditions

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{b}}}_{{v}_{\max }} & = & \displaystyle \frac{{\sin }^{2{v}_{\max }+2}{\boldsymbol{\alpha }}}{4{v}_{\max }+10}\left({{\boldsymbol{f}}}_{{v}_{\max }}({k}^{2},{\boldsymbol{\alpha }})+2{\boldsymbol{q}}{\left({k}^{2},{\boldsymbol{\alpha }}\right)}^{3}\right)\\ {{\boldsymbol{c}}}_{{v}_{\max }} & = & \displaystyle \frac{1}{2}{\boldsymbol{q}}{\left({k}^{2},{\boldsymbol{\alpha }}\right)}^{5}{\sin }^{2{v}_{\max }}{\boldsymbol{\alpha }},\end{array}\end{eqnarray} \tag{ C62 }$$

where

$$\begin{eqnarray}&&\begin{array}{l}{{\boldsymbol{f}}}_{v}({k}^{2},{\boldsymbol{\alpha }})\\ \equiv \displaystyle \frac{3}{v+1}{\,}_{2}{{\boldsymbol{F}}}_{1}\left(-\displaystyle \frac{1}{2},v+1;v+2;1-{\boldsymbol{q}}{\left({k}^{2},{\boldsymbol{\alpha }}\right)}^{2}\right),\end{array}\end{eqnarray} \tag{ C63 }$$

and ₂ F ₁(a, b; c; z ) is the Gauss hypergeometric function, ⁷¹ which we compute via its series definition. We recurse downward in b and c simultaneously:

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{b}}}_{v}({k}^{2},{\boldsymbol{\alpha }}) & = & \displaystyle \frac{1-{\boldsymbol{q}}{\left({k}^{2},{\boldsymbol{\alpha }}\right)}^{2}}{{\sin }^{2}{\boldsymbol{\alpha }}}\left(1+\displaystyle \frac{5}{2v+2}\right)\\ & & \times \,{{\boldsymbol{b}}}_{v+1}({k}^{2},{\boldsymbol{\alpha }})+\displaystyle \frac{{{\boldsymbol{c}}}_{v+1}({k}^{2},{\boldsymbol{\alpha }})}{v+1}\\ {{\boldsymbol{c}}}_{v}({k}^{2},{\boldsymbol{\alpha }}) & = & \displaystyle \frac{{{\boldsymbol{c}}}_{v+1}({k}^{2},{\boldsymbol{\alpha }})}{{\sin }^{2}{\boldsymbol{\alpha }}}.\end{array}\end{eqnarray} \tag{ C64 }$$

The final integral we must solve is that corresponding to ${{\mathbb{p}}}_{2}$ (μみゅー = νにゅー = 1). As in Luger et al. (2019a), this is the integral of the linear limb darkening term, whose solution must be handled separately due to the fact that the corresponding antiderivative in Equation (C25) is not a polynomial in x, y, and z(x, y); also see Agol et al. (2020). The integral we must solve is

$$\begin{eqnarray}\begin{array}{rcl}{{\mathbb{p}}}_{2} & = & {\int }_{{\boldsymbol{\phi }}}\displaystyle \frac{1}{3}\left(\displaystyle \frac{1-z{\left({r}_{{\rm{o}}}\cos \varphi ,{b}_{{\rm{o}}}+{r}_{{\rm{o}}}\sin \varphi \right)}^{3}}{1-z{\left({r}_{{\rm{o}}}\cos \varphi ,{b}_{{\rm{o}}}+{r}_{{\rm{o}}}\sin \varphi \right)}^{2}}\right)\\ & & \times \,\left({r}_{{\rm{o}}}^{2}+{b}_{{\rm{o}}}{r}_{{\rm{o}}}\sin \varphi \right)d\varphi ,\end{array}\end{eqnarray} \tag{ C65 }$$

where z is the usual Cartesian coordinate (Equation (A6)). The solution is tricky, but fortunately a similar integral was solved in Equation (34) of Pál (2012). Adapting their solution to our formalism, we obtain

$$\begin{eqnarray}\begin{array}{rcl}{{\mathbb{p}}}_{2} & = & \displaystyle \frac{1}{3}\left({c}_{0}+{c}_{1}{\rm{\Delta }}{\boldsymbol{F}}\left({\boldsymbol{\alpha }}\left|\,\displaystyle \frac{1}{{k}^{2}}\right.\right)\right.\\ & & +\,{c}_{2}{\rm{\Delta }}{\boldsymbol{E}}\left({\boldsymbol{\alpha }}\left|\,\displaystyle \frac{1}{{k}^{2}}\right.\right)\\ & & +\left.\,{c}_{3}{\rm{\Delta }}{\boldsymbol{\Pi }}\left(n;{\boldsymbol{\alpha }}\left|\,\displaystyle \frac{1}{{k}^{2}}\right.\right)\right),\end{array}\end{eqnarray} \tag{ C66 }$$

where ⁷²

$$\begin{eqnarray}\begin{array}{rcl}{c}_{0} & = & {\rm{\Delta }}\left\{-\arctan 2\left(-({b}_{{\rm{o}}}-{r}_{{\rm{o}}})\cos {\boldsymbol{\alpha }},({b}_{{\rm{o}}}+{r}_{{\rm{o}}})\sin {\boldsymbol{\alpha }}\right)+{\boldsymbol{\alpha }}\right.\\ & & -\,\displaystyle \frac{4}{3}{b}_{{\rm{o}}}{r}_{{\rm{o}}}\sin {\boldsymbol{\alpha }}\cos {\boldsymbol{\alpha }}\sqrt{1-{\left({b}_{{\rm{o}}}-{r}_{{\rm{o}}}\right)}^{2}-4{b}_{{\rm{o}}}{r}_{{\rm{o}}}{\sin }^{2}\alpha }\\ & & +\left.{\boldsymbol{\delta }}({b}_{{\rm{o}}},{r}_{{\rm{o}}},{\boldsymbol{\alpha }})\right\}\\ {c}_{1} & = & \displaystyle \frac{1+{b}_{{\rm{o}}}^{4}-{b}_{{\rm{o}}}^{2}(5+2{r}_{{\rm{o}}}^{2})+{r}_{{\rm{o}}}^{4}+{r}_{{\rm{o}}}^{2}}{3\sqrt{1-{\left({b}_{{\rm{o}}}-{r}_{{\rm{o}}}\right)}^{2}}}\\ {c}_{2} & = & \displaystyle \frac{({b}_{{\rm{o}}}^{2}+7{r}_{{\rm{o}}}^{2}-4)\sqrt{1-{\left({b}_{{\rm{o}}}-{r}_{{\rm{o}}}\right)}^{2}}}{3}\\ {c}_{3} & = & \displaystyle \frac{{b}_{{\rm{o}}}+{r}_{{\rm{o}}}}{({b}_{{\rm{o}}}-{r}_{{\rm{o}}})\sqrt{1-{\left({b}_{{\rm{o}}}-{r}_{{\rm{o}}}\right)}^{2}}},\end{array}\end{eqnarray} \tag{ C67 }$$

and ⁷³

$$\begin{eqnarray}\delta ({b}_{{\rm{o}}},{r}_{{\rm{o}}},\alpha )=\left\{\begin{array}{ll}-2\pi & \alpha \gt \displaystyle \frac{3\pi }{2}\,\mathrm{and}\,{b}_{{\rm{o}}}\gt {r}_{{\rm{o}}}\\ +2\pi & \alpha \gt \displaystyle \frac{3\pi }{2}\,\mathrm{and}\,{b}_{{\rm{o}}}\lt {r}_{{\rm{o}}}\\ 0 & \mathrm{otherwise}.\end{array}\right..\end{eqnarray} \tag{ C68 }$$

The quantities F ( αあるふぁ ∣ 1/k²) and E ( αあるふぁ ∣ 1/k²) are the same incomplete elliptic integrals as those in Appendix C.6.2, ⁷⁴ while

$$\begin{eqnarray}&&{\rm{\Pi }}(n;\psi \,| \,m)\equiv {\int }_{0}^{\psi }\displaystyle \frac{d\varphi }{(1-n{\sin }^{2}\varphi )\sqrt{1-m{\sin }^{2}\varphi }}\end{eqnarray} \tag{ C69 }$$

is the incomplete elliptic integral of the third kind, with

$$\begin{eqnarray}&&n=-\displaystyle \frac{4{b}_{{\rm{o}}}{r}_{{\rm{o}}}}{{\left({r}_{{\rm{o}}}-{b}_{{\rm{o}}}\right)}^{2}}.\end{eqnarray} \tag{ C70 }$$

While stable algorithms exist to evaluate Πぱい(n; ψぷさい ∣ m) (e.g., Bulirsch 1969), ⁷⁵ we find that the parameterization above has poor numerical stability, particularly in the vicinity of the singular points b_o = r_o and b_o = 1 + r_o. In practice, we find that numerical evaluation of Equation (C65) via Gaussian quadrature is more numerically stable and just as computationally efficient as the procedure outlined above.

In this section we present a solution to Equation (C28), the line integral along the day/night terminator of the occulted body. As before, the first thing we must do is derive expressions for the integration limits ξくしー . Depending on the integration case, these limits will corresponds to the point of intersection between the terminator and the limb of the occultor and/or the point of intersection between the terminator and the limb of the occulted body. The former is given by

$$\begin{eqnarray}&&{{\boldsymbol{\xi }}}_{{\bf{0}}}=\arctan 2\left(\sqrt{1-{\boldsymbol{x}}{{\prime\prime} }^{2}},{\boldsymbol{x}}^{\prime\prime} \right),\end{eqnarray} \tag{ C71 }$$

where x '' are the roots of the quartic (Equation (C3)). ⁷⁶ The latter is given by

$$\begin{eqnarray}{\xi }_{1}=\left\{\begin{array}{ll}0 & {\left(1-{x}_{{\rm{o}}}^{{\prime\prime} }\right)}^{2}+{{y}_{{\rm{o}}}^{{\prime\prime} }}^{2}\lt {r}_{{\rm{o}}}^{2}\\ \pi & \mathrm{otherwise}.\end{array}\right..\end{eqnarray} \tag{ C72 }$$

As before, these angles are then wrapped to the range [0, 2πぱい) and sorted into the vector ξくしー such that the integration is performed counter-clockwise about the center of the integration region. ⁷⁷ The middle panel of Figure 16 shows a case where ξくしー₀ = 965 corresponds to the point of intersection between the occultor limb and the terminator and ξくしー₁ = 0° corresponds to the point where the terminator extends onto the backside of the body. Note, importantly, that unlike ϕ , the angle ξくしー is not measured between the horizontal and a point on the curve of ${ \mathcal T }$. Recall that ξくしー is an angular parameter of the ellipse, so it is measured in the same way as the eccentric anomaly in a Keplerian orbit: it is the angle between the semimajor axis of the ellipse and the perpendicular projection of a point on the ellipse onto the unit circle (see Figure 16).

The solution to Equation (C28) involves repeated application of the binomial theorem. If we define the quantities

$$\begin{eqnarray*}&&{Z}_{j,k}^{u,v}(b,\theta ^{\prime} )=\left(\displaystyle \genfrac{}{}{0em}{}{u}{j}\right)\left(\displaystyle \genfrac{}{}{0em}{}{v}{k}\right){\left(-1\right)}^{v-k}{b}^{u+v-j-k}{\sin }^{v+j-k}\theta ^{\prime} {\cos }^{u-j+k}\theta ^{\prime} \end{eqnarray*}$$

and

$$\begin{eqnarray}\delta (b,\xi )=\left\{\begin{array}{ll}0 & 0\leqslant \xi \lt \displaystyle \frac{\pi }{2}\\ \pi & \displaystyle \frac{\pi }{2}\leqslant \xi \lt \pi \\ 2| b| {b}_{{\rm{c}}} & \pi \leqslant \xi \lt \displaystyle \frac{3\pi }{2}\\ \pi +2| b| {b}_{{\rm{c}}} & \displaystyle \frac{3\pi }{2}\leqslant \xi \lt 2\pi ,\end{array}\right.,\end{eqnarray} \tag{ C73 }$$

we may express the solution ⁷⁸ to the ${{\mathbb{t}}}_{n}$ integral as ⁷⁹

$$\begin{eqnarray}&&{{\mathbb{t}}}_{n}=\left\{\begin{array}{ll}\begin{array}{l}b\cos \theta ^{\prime} \displaystyle \sum _{j=0}^{\tfrac{\nu }{2}}\displaystyle \sum _{k=0}^{\tfrac{\mu +2}{2}}{Z}_{j,k}^{\tfrac{\nu }{2},\tfrac{\mu +2}{2}}(b,\theta ^{\prime} )\,{{\mathbb{H}}}_{j+k+1,l+1-j-k}({\boldsymbol{\xi }})\\ -\,\sin \theta ^{\prime} \displaystyle \sum _{j=0}^{\tfrac{\nu }{2}}\displaystyle \sum _{k=0}^{\tfrac{\mu +2}{2}}{Z}_{j,k}^{\tfrac{\nu }{2},\tfrac{\mu +2}{2}}(b,\theta ^{\prime} )\,{{\mathbb{H}}}_{j+k,l+2-j-k}({\boldsymbol{\xi }})\end{array} & \mu ,\nu ,\mathrm{even}\\ \begin{array}{l}\displaystyle \frac{1}{3}{\rm{\Delta }}\left\{\arctan \left(\displaystyle \frac{| b| \sin {\boldsymbol{\xi }}}{\cos {\boldsymbol{\xi }}}\right)\right.\\ \,\,\,\,-\,\mathrm{sgn}\left(\sin {\boldsymbol{\xi }}\right)\left(\arctan \left(\displaystyle \frac{{\left(\tfrac{\sin {\boldsymbol{\xi }}}{1+\cos {\boldsymbol{\xi }}}\right)}^{2}+2{b}^{2}-1}{2| b| {b}_{{\rm{c}}}}\right)+| b| {b}_{{\rm{c}}}\cos {\boldsymbol{\xi }}\right)\\ \left.+\,{\boldsymbol{\delta }}(b,{\boldsymbol{\xi }})\right\}\end{array} & \mu =\nu =1\\ \begin{array}{l}-{{bb}}_{c}^{3}\sin \theta ^{\prime} \displaystyle \sum _{j=0}^{l-2}{Z}_{0,j}^{0,l-2}(b,\theta ^{\prime} )\,{{\mathbb{H}}}_{j+1,l+1-j}({\boldsymbol{\xi }})\\ \ -\,{b}_{{\rm{c}}}^{3}\cos \theta ^{\prime} \displaystyle \sum _{j=0}^{l-2}{Z}_{0,j}^{0,l-2}(b,\theta ^{\prime} )\,{{\mathbb{H}}}_{j,l+2-j}({\boldsymbol{\xi }})\end{array} & \nu \,\mathrm{odd},\mu =1,l\,\mathrm{even}\\ \begin{array}{l}-{{bb}}_{c}^{3}{\sin }^{2}\theta ^{\prime} \displaystyle \sum _{j=0}^{l-3}{Z}_{0,j}^{0,l-3}(b,\theta ^{\prime} )\,{{\mathbb{H}}}_{j+2,l-j}({\boldsymbol{\xi }})\\ \ -\,{{bb}}_{{\rm{c}}}^{3}{\cos }^{2}\theta ^{\prime} \displaystyle \sum _{j=0}^{l-3}{Z}_{0,j}^{0,l-3}(b,\theta ^{\prime} )\,{{\mathbb{H}}}_{j,l+2-j}({\boldsymbol{\xi }})\\ -\,{b}_{{\rm{c}}}^{5}\sin \theta ^{\prime} \cos \theta ^{\prime} \displaystyle \sum _{j=0}^{l-3}{Z}_{0,k}^{0,l-3}(b,\theta ^{\prime} )\,{{\mathbb{H}}}_{j+1,l+1-j}({\boldsymbol{\xi }})\end{array} & \nu \,\mathrm{odd},\mu =1,l\,\mathrm{odd}\\ \begin{array}{l}{{bb}}_{c}^{3}\cos \theta ^{\prime} \displaystyle \sum _{j=0}^{\tfrac{\nu }{2}}\displaystyle \sum _{k=0}^{\tfrac{\mu +2}{2}}{Z}_{j,k}^{\tfrac{\nu -1}{2},\displaystyle \frac{\mu -3}{2}}(b,\theta ^{\prime} )\,{{\mathbb{H}}}_{j+k+1,l+1-j-k}({\boldsymbol{\xi }})\\ -\,{b}_{{\rm{c}}}^{3}\sin \theta ^{\prime} \displaystyle \sum _{j=0}^{\tfrac{\nu }{2}}\sum _{k=0}^{\tfrac{\mu +2}{2}}{Z}_{j,k}^{\tfrac{\nu -1}{2},\displaystyle \frac{\mu -3}{2}}(b,\theta ^{\prime} )\,{{\mathbb{H}}}_{j+k,l+2-j-k}({\boldsymbol{\xi }})\end{array} & \mathrm{otherwise}\end{array}\right..\end{eqnarray} \tag{ C74 }$$

The solution to Equation (C74) depends on the matrix ${\mathbb{H}}$, whose components are given by the integral

$$\begin{eqnarray}&&{{\mathbb{H}}}_{u,v}({\boldsymbol{\xi }})={\int }_{{\boldsymbol{\xi }}}{\cos }^{u}\varphi {\sin }^{v}\varphi d\varphi .\end{eqnarray} \tag{ C75 }$$

The ${\mathbb{H}}$ integral is the same as that in Equation (D27) of Luger et al. (2019a), except for a change in the limits of integration. We can compute this integral recursively given four lower boundary conditions:

$$\begin{eqnarray}\begin{array}{rcl}{{\mathbb{H}}}_{\mathrm{0,0}}({\boldsymbol{\xi }}) & = & {\rm{\Delta }}{\boldsymbol{\xi }}\\ {{\mathbb{H}}}_{\mathrm{1,0}}({\boldsymbol{\xi }}) & = & {\rm{\Delta }}\sin {\boldsymbol{\xi }}\\ {{\mathbb{H}}}_{\mathrm{0,1}}({\boldsymbol{\xi }}) & = & -{\rm{\Delta }}\cos {\boldsymbol{\xi }}\\ {{\mathbb{H}}}_{\mathrm{1,1}}({\boldsymbol{\xi }}) & = & -\displaystyle \frac{{\rm{\Delta }}{\cos }^{2}{\boldsymbol{\xi }}}{2}.\end{array}\end{eqnarray} \tag{ C76 }$$

The remaining terms may be computed ⁸⁰ by upward recursion using the relations

$$\begin{eqnarray}&&{{\mathbb{H}}}_{u,v}({\boldsymbol{\xi }})=\displaystyle \frac{-{\rm{\Delta }}\left({\cos }^{u+1}{\boldsymbol{\xi }}{\sin }^{v-1}{\boldsymbol{\xi }}\right)+(v-1){{\mathbb{H}}}_{u,v-2}({\boldsymbol{\xi }})}{u+v}\end{eqnarray} \tag{ C77 }$$

for u < 2, v ≥ 2 and

$$\begin{eqnarray}&&{{\mathbb{H}}}_{u,v}({\boldsymbol{\xi }})=\displaystyle \frac{{\rm{\Delta }}\left({\cos }^{u-1}{\boldsymbol{\xi }}{\sin }^{v+1}{\boldsymbol{\xi }}\right)+(u-1){{\mathbb{H}}}_{u-2,v}({\boldsymbol{\xi }})}{u+v}\end{eqnarray} \tag{ C78 }$$

for all remaining terms. ^{81 , 82}

The final line integral we must solve is the integral along the boundary of the occulted body, Equation (C29). Fortunately, this is also the easiest of the three. The limits of integration λらむだ correspond to the point at which the terminator crosses from the dayside to the nightside,

$$\begin{eqnarray}{\lambda }_{0}=\left\{\begin{array}{ll}\theta ^{\prime} & {\cos }^{2}\theta ^{\prime} +{\left(\sin \theta ^{\prime} -{b}_{{\rm{o}}}\right)}^{2}\lt {r}_{{\rm{o}}}^{2}\\ \theta ^{\prime} +\pi & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ C79 }$$

and the point of intersection between the limb of the occultor and the limb of the occulted body, ⁸³

$$\begin{eqnarray}&&{\lambda }_{1}=\displaystyle \frac{\pi }{2}\pm \left(\arcsin \left(\displaystyle \frac{1-{r}_{{\rm{o}}}^{2}+{b}_{{\rm{o}}}^{2}}{2{b}_{{\rm{o}}}}\right)-\displaystyle \frac{\pi }{2}\right),\end{eqnarray} \tag{ C80 }$$

where the sign is chosen such that the point $(\cos {\lambda }_{1},\sin {\lambda }_{1})$ is on the dayside of the occulted body. ⁸⁴ As before, these angles are wrapped to the range [0, 2πぱい) and placed in the vector λらむだ such that the line integral is taken in the counter-clockwise direction about the center of the integration region. The right panel of Figure 16 shows a case where λらむだ₀ = 75° and λらむだ₁ = 1305. Both angles are measured counter-clockwise from the x'-axis.

Given λらむだ , the solution to Equation (C29) is straightforward:

$$\begin{eqnarray}{{\mathbb{q}}}_{n}=\left\{\begin{array}{ll}{{\mathbb{H}}}_{\tfrac{\mu +4}{2},\displaystyle \frac{\nu }{2}}({\boldsymbol{\lambda }}) & \mu ,\nu ,\mathrm{even}\\ \displaystyle \frac{1}{3}{\rm{\Delta }}{\boldsymbol{\lambda }} & \qquad \mu =\nu =1\\ 0 & \mathrm{otherwise},\end{array}\right.\end{eqnarray} \tag{ C81 }$$

where the matrix ${\mathbb{H}}$ is given by Equation (C75). ⁸⁵