© ESO 2018

1 Introduction

Since the detection of the first planets beyond our solar system (Campbell et al. 1988; Wolszczan & Frail 1992), the number of discoveries has steadily increased, yielding almost 4000 confirmed exoplanets and 2500 unconfirmed, candidate exoplanets to this day (Han et al. 2014). Exoplanet space telescopes, such as ESA’s CHEOPS (CHaracterising ExOPlanet Satellite) and PLATO (PLAnetary Transits and Oscillations of stars), and NASA’s TESS (Transiting Exoplanet Survey Satellite, are dedicated to finding exoplanets around bright, nearby stars. The relatively small distances to these stars and their planets combined with the high sensitivity of these missions and the upcoming NASA’s JWST (James Webb Space Telescope) and ESA’s ARIEL (Atmospheric Remote-Sensing Infrared Exoplanet Large-survey; Tinetti et al. 2016) missions will allow the search for lunar companions and planetary rings.

The continuous increase in instrument precision, stability, and spatial resolution has enabled a new generation of ground-based instruments, such as the Gemini Planet Imager (GPI) instrument (see Macintosh et al. 2014) on the Gemini North telescope, the Spectro-Polarimetric High-contrast Exoplanet Research (SPHERE) instrument (see Beuzit et al. 2006) on ESO’s Very Large Telescope (VLT), and the proposed Exoplanet Imaging Camera and Spectrograph (EPICS; see Keller et al. 2010; Gratton et al. 2011) on the European Extremely Large Telescope (E-ELT), which is under construction by ESO. These high-contrast instruments use direct imaging of planetary radiation to not only detect but also characterize exoplanetary systems through a combination of spectroscopy and polarimetry techniques. Using GPI and SPHERE, respectively, Macintosh et al. (2015) and Wagner et al. (2016) announced the discoveries of young, and thus hot, Jovian planets, whose atmospheric properties and orbits were characterized using near-infrared spectroscopy.

Direct detection of extrasolar bodies presents a major challenge as their observed radiation, both emitted and reflected, is very weak compared to that of the host star. In addition, the angular distance from the planet to the star is extremely small. Consequently, the vast majority of exoplanets have only been detected indirectly. In contrast, direct imaging of the planet can reveal a wealth of information on the planet properties that cannot be obtained through other methods, such as lower atmospheric composition and, for rocky planets, their surface coverage.

Despite not yet having been exploited, great hope is placed in the polarimetry capabilities of current and future telescopes as powerful tools for detecting and characterizing exoplanets (see Seager et al. 2000; Hough et al. 2003; Hough & Lucas 2003; Saar & Seager 2003; Stam 2003, 2008; Stam et al. 2004; Snik & Keller 2013). Previous works in this field involve the modeling of stellar polarization during planetary transits (Carciofi & Magalhães 2005; Kostogryz et al. 2011, 2015; Wiktorowicz & Laughlin 2014; Sengupta 2016), the modeling of light curves and polarization of starlight reflected signals in the visible range of Earth-like planets (Stam 2008; Karalidi et al. 2012; Rossi & Stam 2017) and giant Jupiter-like planets (Seager et al. 2000; Stam et al. 2004), as well as the modeling of exoplanetary atmospheres in the infrared (De Kok et al. 2011; Marley & Sengupta 2011), which demonstrated the usefulness of direct observations for exoplanet characterization. More recently, Bott et al. (2016) reported linear polarization observations of the hot Jupiter system HD 189733, and Ginski et al. (2018) announced the detection of planetary thermal radiation that is polarized upon reflection by circumstellar dust. Indeed, polarization has also been proposed as a means for exomoon detection: Sengupta & Marley (2016) studied the effects of a satellite transiting its hot host planet in the polarization signal of (infrared) thermally emitted radiation for the case of homogeneous, spherically symmetric cloudy planets. Studying exomoons can improve our understanding in particular:

• 1.

  Planet formation. Solar system moons appear to support diverse formation histories. For instance, Titan might have formed from circumplanetary debris, while the Moon, Phobos, and Deimos suggest a cumulative bombardment (Rosenblatt et al. 2016; Rufu et al. 2017). Triton might have been captured by Neptune (Agnor & Hamilton 2006), while collisions are thought to have altered the relative alignment between Uranus and its moons (Morbidelli et al. 2012). Indeed, studying solar system moons gives essential insights into formation mechanisms and evolution (see Heller 2017, and references therein). Exomoon research would enable the refinement of planet formation theories in a waynot achievable by studying exoplanets alone.

• 2.

  Extra-solar system characterization. Studying exomoons will not only provide information on lunar orbits and physical properties, but will also allow us to constrain planet characteristics such as mass, oblateness, and rotation axis (Barnes & Fortney 2003; Kipping et al. 2010; Schneider et al. 2015). A signal of a planet–moon system could be interpreted as that of a planet alone, resulting, for example, in an overestimation of the planet mass and effective temperature (Williams & Knacke 2004), and/or an anomalous composition from spectroscopy (Schneider et al. 2015). Extrasolar system characterization would indeed require analysis of all of its elements: planets, moons, rings, and exozodiacal dust.

• 3.

  Exoplanet and exomoon habitability. A moon may influence its planet’s habitability (Benn 2001), and moons of giant exoplanets within the stellar habitable zone (HZ) might host habitable environments (Canup & Ward 2006). Reynolds et al. (1987) and Heller & Barnes (2015) mention the role of a moon’s orbit on the presence of liquid, life-supporting water. Indeed, tidal heating could maintain surface temperatures compatible with life on large moons around cold giant planets (Scharf 2006). Lehmer et al. (2017) show that small moons could retain atmospheres over limited time periods, while Ganymede-sized moons in a stellar HZ could hold atmospheres and surface water indefinitely. Although radiation in a giant planet’s magnetic field and eclipses could threaten local conditions for life (Heller 2012; Heller & Barnes 2013; Forgan & Yotov 2014), exomoons are interesting targets in the search for extraterrestrial life.

Led by Kipping’s Hunt for Exomoons with Kepler, which uses a combination of photometric transits, transit timing variations (TTV) and transit duration variations (TDV) data (Sartoretti & Schneider 1999; Szabó et al. 2006; Simon et al. 2007, 2015; Kipping 2009; Kipping et al. 2015), and Hippke’s search using the orbital sampling effect (OSE; Heller 2014; Hippke 2015; Heller et al. 2016), the search for exomoons is in its initial phase. Mars-sized and possibly even Ganymede-sized satellites could be traceable in archived Kepler data (Heller et al. 2014). Unfortunately, as yet no exomoons have been confirmed.

In this paper, we use numerical simulations to show how an exomoon could influence the flux and degree of polarization of the starlight that is reflected by an Earth-like exoplanet, using the following outline. In Sect. 2, we describe the numerical code to compute the various geometries of the exoplanet–exomoon system that are required for our radiative transfer computations and the radiative transfer code to compute the reflected fluxes and polarization for a given exoplanet–exomoon system. In Sect. 3, we present computed flux and polarization phase functions at 450 nm for an Earth-like planet (with a Lambertian reflecting surface and a gaseous atmosphere) with a Moon-like satellite (with a Lambertian reflecting surface) in an edge-on geometry. Finally, in Sect. 4, we summarize and discuss our findings and their implications.

2 Computing the reflected starlight

2.1 Stokes vectors and polarization

We describe the flux and polarization of starlight that is reflected by a body, with a Stokes vector, (see, e.g., Hansen & Travis 1974), (1)

with F the total flux, Q and U the linearly polarized fluxes, and V the circularly polarized flux, all with dimensions W m⁻². In principle, these fluxes are wavelength dependent. However, we will not explicitly include the wavelength in the dimensions, because we focus on a single wavelength region. Fluxes Q and U are defined with respect to a reference plane, for which we use the planetary (or lunar) scattering plane, which contains the observer, and the centers of the planet (or moon) and the star. We do not compute the circularly polarized flux V, because it is usually much smaller than the linearly polarized fluxes (see Hansen & Travis 1974; Kawata 1978; Rossi & Stam 2018), and because ignoring V does not lead to significant errors in the computation of F, Q, and U (see Stam & Hovenier 2005). The light of the star is assumed to be unpolarized (see Kemp et al. 1987), and is given by F₀ = F₀1, with 1 the unit column vector and πF₀ the flux measured perpendicular to the light’s propagation direction. If the orbit of the barycenter of the planet–moon system around the star is eccentric, the incident flux varies along the orbit. Our standard stellar flux, πF₀, is defined with respect to the periapsis of the orbit of the system’s barycenter.

The degree of linear polarization, P, of vector F is defined as (2)

and the direction of polarization, χ, with respectto the reference plane can be computed from (3)

where χ is chosen such that 0 ≤ χ < π, while cos2χ and Q have the same sign (Hansen & Travis 1974; Hovenier et al. 2004).

2.2 Disk-integrated reflected Stokes vectors

We computethe reflected Stokes vector F of the spatially unresolved planet–moon system as a summation of the reflected Stokes vectors F^p and F^m of, respectively,the planet and the moon (the pair is spatially resolved from the star): (4)

Vectors F^p and F^m are disk-integrated vectors that include the effects of eclipses and transits. They are normalized such that the total fluxes reflected by the planet and moon at a phase angle α = 0° and withoutshadows and/or eclipses on their disks, equal the planet’s and moon’s geometric albedos, respectively (see Stam et al. 2006). Furthermore, R_p and R_m are the radiiof the (spherical) planet and moon, respectively.

Vectors F and F^p in Eq. (4) are defined with respect to the planetary scattering plane, while F^m is defined with respect to the lunar scattering plane. Depending on the orientation of the lunar orbit, the lunar scattering plane can have a different orientation than the planetary scattering plane. Matrix L in Eq. (4) rotates F^m from the lunar to the planetary scattering plane. It is given by (see Hovenier & van der Mee 1983) (5)

with ψ the rotation angle measured in the clockwise direction from the lunar to the planetary scattering plane when looking toward the moon (0 ≤ ψ < π).¹

To compute the disk-integrated vectors F^p and F^m, we divide the disks of the planet and the moon as seen by the observer into a grid of equally sized, square pixels (see Fig. 1). The number of pixels on the planetary disk is N_p and that on the lunar disk is N_m. A given pixel will contribute to a disk signal when its center is within the disk radius. Obviously, the larger the number of pixels (and the smaller each pixel), the better the approximation of the curved limb of the disk, the terminator, and the shadows, such as those due to eclipses (see Appendix C for insight into the effect of the number of pixels on the computed signals). The disk-integrated vectors are obtained by summing up the contributions of the individual pixels across the disk, fully taking into account shadows and eclipses, that is, (6)

where “x” is either “p” or “m”. Factor π∕N_x is the surface area per pixel. Furthermore, is the reflected Stokes vector for the ith pixel on the planet (x = p) or moon (x = m), the computation of which is described in Sect. 2.3. Matrix L is a rotation matrix (see Eq. (5)) that is used for rotating the local Stokes vector that is defined with respect to the local reference plane, or to the planetary or lunar scattering plane. Factor b_i accounts for the visibility of pixel i: if b_i = 1, the pixel is visible to the observer, and if b_i = 0, it is invisible due to a transiting body. Factor c_i accounts for the dimming of the local incident stellar flux due to a (partial) eclipse: c_i = 0.0 indicates that pixel i is eclipsed and receives no flux, and c_i = 1.0 indicates that the pixel is not eclipsed. For partial (penumbral) eclipses, 0.0 < c_i < 1.0. The computation of factors b_i and c_i is described in Sect. 2.4. Factor d_i, finally, indicates the decrease of the standard incident stellar flux πF₀ due to an increase of the distance to the star, according to (7)

where r_ref is the reference distance at which the standard stellar flux is defined and r_is is the actual distance between pixel i and the star.

Fig. 1 Three-dimensional view (left panel) and projection onto the x0 –y0 plane (right panel) of the discretized planet or moon. The z0-axis points toward the observer. The orientation of the x0 and y0 axes with respect to the disk of the planet or moon can be chosen arbitrarily.

2.3 Locally reflected starlight

The Stokes vector of the starlight that is reflected by pixel i on the planet or moon is computed using (see Hansen & Travis 1974) (8)

with θ_i the angle between the local zenith direction and the local direction to the observer, θ_0i the angle between the local zenith direction and the local direction to the star, and ϕ_i − ϕ_0i the local azimuthal difference angle, that is, the angle between the plane containing the local zenith direction and the local direction to the observer and the plane containing the local zenith direction and the local direction to the star (see de Haan et al. 1987; Rossi et al. 2018). Furthermore, is the first column of the local reflection matrix of the planet or moon. Only the first column is needed because the incident starlight is assumed to be unpolarized (cf. Sect. 2.1). For a given pixel, the illumination and viewing angles, and thus , depend on the position of the planet or moon with respect to the star and to each other. Local reflection matrix also depends on the local composition and structure of the atmosphere and/or surface of the reflecting body. We compute reflected starlight for an Earth–Moon-like planetary system, keeping the reflection models for the Earth and the moon simple to avoid introducing too many details that increase computational times while not adding insight into the observable signals.

Our model planet has a flat, Lambertian (i.e., isotropically and depolarizing) reflecting surface with an albedo, a_surf, of 0.3. The surface is overlaid by an atmosphere that is assumed to consist of only gas. We compute the atmospheric optical thickness at a given wavelength λ, using a model atmosphere consisting of 32 layers, with the ambient pressure and temperature according to a mid-latitude summer profile (McClatchey et al. 1972). The surface pressure is 1.0 bars. The molecular scattering optical thickness of an atmospheric layer at wavelength λ is calculated according to (9)

with the molecular scattering cross section (in m²) and N the molecular column number density (in m⁻²) of the atmospheric layer. The molecular scattering cross section is calculated according to (10)

with N_L being Loschmidt’s number at standard pressure and temperature, n the wavelength-dependent refractive index of dry air under standard pressure and temperature, and δ the depolarization factor of the atmospheric gas (see Stam 2008, and references therein for the values that have been chosen for the various parameters). To calculate the molecular column number density N, we assume hydrostatic equilibrium in each atmospheric layer; thus, (11)

with Δp the difference between the pressure at the bottom and at the top of the atmospheric layer (in Pa), m the average molecular mass in the layer (in kg), and g the acceleration of gravity (in m s⁻²). The atmospheric optical thickness at a given wavelength λ is calculated by adding the values of for all atmospheric layers at that wavelength (note that for a model atmosphere containing only gas, the radiative transfer of incident sunlight only depends on the total optical thickness, not on the vertical distribution of the optical thickness). The total atmospheric optical thickness at 450 nm, the wavelength of our interest, is 0.14. At this wavelength, there is no significant absorption by atmospheric gases in the Earth’s atmosphere (see Stam 2008, for sample spectra). The single scattering albedo of the gaseous molecules can thus be assumed to equal 1.0. At this short wavelength, the horizontal inhomogeneities of the Earth’s surface and the contributions of clouds and aerosol to the reflected signal are relatively small (see Stam 2008, for simulations of the Earth’s signal at 440 nm). Our model moon has no atmosphere above its flat, Lambertian (i.e., isotropic and depolarizing) reflecting surface, with a_surf = 0.1 (Williams 2017).

The computation of the local illumination and viewing geometries θ_i, θ_0i, and ϕ_i − ϕ_0i is described in Appendix A. Given these angles and the planet’s atmosphere–surface model, we use PyMieDAP² (Rossi et al. 2018), an efficient radiative transfer code based on the adding–doubling algorithm described by de Haan et al. (1987). PyMieDAP fully includes polarization for all orders of scattering, and assumes a locally plane–parallel atmosphere–surface model to compute for every pixel on the planet. The computed locally reflected Stokes vector, , is defined with respect to the local meridian plane, that is, the plane through the local zenith and the local direction toward the observer. For each illuminated pixel on the moon, . A detailed description of PyMieDAP including benchmark results can be found in Rossi et al. (2018).

Results of our radiative transfer code have been compared against results presented in, for example, Stam (2008), Stam et al. (2004; 2006; who all used the same adding–doubling code, but an entirely different disk integration algorithm), and Karalidi et al. (2012; who used their own version of an adding–doubling code and an independent disk–integration method). Buenzli & Schmid (2009) and Stolker et al. (2017) each compared their own independently implemented Monte Carlo radiative transfer codes successfully against results from the code used by Stam et al. (2004, 2006).

2.4 Computing transits and eclipses

As described in Eq. (6), the contribution of the light reflected by a pixel i on the planet or the moon to the disk-integrated Stokes vector F depends on factors b_i and c_i, which account for the pixel’s visibility and dimming, respectively. The values of these factors depend on so-called mutual events, specifically transits, in which one body is (partially) blocking the light that is reflected toward the observer by another body, and eclipses, in which one body is casting a (partial) shadow on the illuminated and visible disk of another body. Limiting ourselves to systems in which a single star is orbited by a planet with a single moon, we distinguish the following four mutual events (cf. Fig. 2):

• 1.

  A planetary eclipse. The moon is between the star and the planet, casting its shadow on the planet, the extent of which depends on the planet–star and moon–star distances, on the stellar, planetary, and lunar radii, and on their orbital positions.

• 2.

  A lunar eclipse. The planet is between the star and the moon, casting its shadow on the moon, the extent of which depends on the planet–star and moon–star distances, on the stellar, planetary, and lunar radii, and on their orbital positions.

• 3.

  A planetary transit. The planet is between the moon and the observer, blocking the view of the moon, the extent of which depends on the planetary and lunar radii, and their orbital positions.

• 4.

  A lunar transit. The moon is between the planet and the observer, occulting a region of the planetary disk, the extent of which depends on the planetary and lunar radii, and their orbital positions.

We exclude planetary and lunar transits of the star, that is, the epochs in which these bodies move in front of or behind the star. Numerical simulations of transiting planets with moons have been published by Kipping (2011). Modeling the transmission and scattering of starlight in the planetary atmosphere during those epochs (which is not included in the work by Kipping 2011) requires a fully spherical atmosphere model instead of a locally plane-parallel one (de Kok & Stam 2012) and falls outside the scope of this paper.

For our computation of the effects of transits of the planet in front of the moon and vice versa on the flux and polarization of the reflectedstarlight, we assume that the bodies are at an “infinite” distance from the observer. For our computation of the effects of the eclipses (i.e., the shadow of one body darkening regions on the other body) on the reflected flux and polarization, we follow the mathematical description of eclipses in the Moon–Earth system as developed by Link (1969), taking into account the sizes of the planet and the moon, their distances and positions with respect to the star, and the size of the stellar disk. The latter is crucial for the modeling of the umbra, antumbra, and penumbra shadow regions (for an example of the umbra and penumbra, see Fig. 2). The contribution of the starlight reflected by pixels in the antumbral or penumbral region of the planet or moon to the total signal is weighted by the depth of the shadow (i.e., factor c_i in Eq. (6)). We ignore stellar limb darkening and the transmission of starlight through the planetary atmosphere during a lunar eclipse.

A detailed description of our numerical computation of eclipses and the factor c_i in Eq. (6) is provided in Appendix B. This computation requires the positions of the planet and the moon with respect to the star across time, and thus the dynamics of the three-body system. The basics of this dynamics is outlined in the next section.

Fig. 2 The mutual events between a planet and its moon. Panel a: a lunar transit. Panel b: a planetary transit. Panel c: a planetary eclipse. Panel d: a lunar eclipse. The positive z0 -axis (in panels aand b) points to the observer. The white–black scale indicates full–null illumination of a body. An arrow indicates the lunar motion around the planet. The distances and radii are not to scale.

2.5 Computing the orbits of the planet and moon

We compute the position vectors of the planet and moon as functions of time for determining the factors b_i, c_i, and d_i of each pixel i and for evaluating the disk integration according to Eq. (6). Both the motions of the planet and its moon around the star depend on their mutual gravitational interactions. Assuming each body attracts as a point mass and neglecting the gravity of other planets and/or moons in the system, our star–planet–moon system is a classical, generic, three-body problem.

A precise computation of the orbital positions in the generic three-body problem requires the numerical propagation of a given set of initial conditions. Instead, we use the “nested two-body” approximation described by Kipping (2011, 2010), which assumes that the orbits of the planet and moon around the planet–moon system barycenter, and the orbit of this barycenter around the star, can all be described by Keplerian orbits. The advantages of the nested two-body approximation are as follows: (1) The solution can be described analytically. (2) The computational time is significantly shorter than with numerical integrations. (3) It provides better insight into the computed orbits as the elements of all orbits can be specified. (4) Unlike the circular restricted three-body problem simplification (see, e.g., Wakker 2015), it can handle eccentric orbits.

As demonstrated by Kipping (2010), the nested two-body approximation is excellent for the generic three-body problem provided 𝔇 ≤ 0.531, where 𝔇 is the moon–planet separation in units of the planet’s Hill’s sphere radius (see, e.g., De Pater & Lissauer 2015). As follows from Domingos et al. (2006) and Kipping (2011), stable, prograde orbiting moons should fulfill 𝔇 ≤ 0.4895, while retrograde orbiting moons can be stable up to 𝔇 ≈ 0.9309. The nested two-body approximation can thus be applied to all prograde orbiting moons, while retrograde orbiting moons are only partially covered, depending on 𝔇. We will limit ourselves to prograde orbiting moons, as we do not expect any influence of the moon’s orbital direction on the magnitude of reflected flux and polarization features, except on their timing.

Figure 3 shows the geometry of the planet–moon system with the reference frames describing the orbit of the planet–moon barycenter around the star, and the orbit of the moon around the barycenter. The nested two-body approximation assumes that the motions of the planet–moon barycenter around the star and that of the moon around the barycenter are independent. An orthonormal, right-handed coordinate system S₁ = {x₁, y₁, z₁} is the reference frame for the observation of the planet–moon–star system, with the star at the origin, and the z₁-axis pointing toward the observer. Plane p₁, through x₁ and y₁, is the plane on the observer’s sky, onto which the pixels (Fig. 1) are projected. Axes x₁ and y₁ have an arbitrary (but fixed) orientation. Coordinate system S₂ = {x₂, y₂, z₂} is the reference frame for the orbit of the barycenter, which lies in plane p₂, through x₂ and y₂. The lunar orbit itself lies in the x₃ –y₃ plane of coordinate system S₃ = {x₃, y₃, z₃}, which is centered at the barycenter’s position. The various orbital elements in these coordinate systems are indicated as follows:

With the orbital parameters of the barycenter and the moon, we compute the true anomalies ν_bs of the barycenter and ν_mb of the moon around the barycenter, respectively, at any time t using Kepler’s equation for elliptical orbits: (12)

where E is the eccentric anomaly and M the mean anomaly of the orbit. The true anomalies of the barycenter and the moon are related to the eccentric anomaly E through (13)

We solve for ν(t) for each orbit in the appropriate reference system by applying the Newton–Raphson method (Wakker 2015) to Eqs. (12) and (13).

We compute the position of the barycenter, r_bs, in coordinate system S₂, and the position of the moon around the barycenter, r_mb, in coordinate system S₃. The absolute position of the moon in coordinate system S₂ is then obtained through (14)

As formulated by Murray & Correia (2010), the position of the barycenter and the moon in S₂ at time t can be put through a series of transformation matrices to yield the position of the barycenter and the moon in the observer’s coordinate system S₁. Further details of these transformation matrices can be found in Kipping (2011, 2010).

After having computed the positions of the planet and the moon in S₁ at time t, we compute the positions of the pixels across the planetary and lunar disks (see Fig. 1), and the angles β_i that are used to rotate locally computed Stokes vectors to the planetary and lunar scattering planes (Eq. (6)), respectively, Then we calculate parameter d_i, which accounts for the change of the standard incident flux due to the changing distance to the star (see Eq. (6)), and the local illumination and viewing angles required for the computation of the Stokes vector of reflected starlight for each pixel seen by the observer. Details on these computations can be found in Appendix A. For each t, we also compute angle ψ to rotate F^m, the disk-integrated Stokes vector for the moon, to the planetary scattering plane (Eq. (4)).

Fig. 3 (a) The planet–moon barycenter orbit around the star (b) The moon orbit around the barycenter. The reference frames and angles used to describe the Keplerian orbits of the planet–moon system barycenter around the parent star (panel a) and the Keplerian orbit of the moon around the planet-moon system barycenter (panel b). Plane p1 (panel a) is the plane of the sky as seen by the observer on the positive z1 -axis. Plane p2 (panels a and b) is the barycenter’s orbital plane, and plane p3 (panel b) is the orbital plane of the moon around the barycenter. Angle i is the orbitalinclination angle, ω the argument of periastron, Ω the right ascension of the ascending node, and ν the true anomaly. Subscript bs refers to the barycenter of the planet–moon system around the star, and mb to the moon around the barycenter. Vectors rbs and rmb are the position vectors of the barycenter and the moon, respectively. The barycenter and the moon are indicated by b and m, respectively.

2.6 Our baseline planet–moon system

In this paper, we focus on planet–moon systems in edge-on geometries, in which the inclination angle of the barycenter’s orbit is 90°, because exoplanets in (near) edge-on orbits are prime targets for space telescopes such as TESS, JWST, PLATO, and CHEOPS, which will all employ the transit method to detect and/or characterize exoplanets, as well as for follow-up missions including telescopes aimed at directly detecting planet signals.

Table 1 lists the orbital elements of our baseline planet–moon system. Both the barycenter’s and the lunar orbit are assumed to be circular (e = 0.0), and their semi-major axes match those of the Earth–Moon system (Williams 2017). We neglect the Earth–barycenter distance, so that a_bs = 1 AU. The semi-major axis of the lunar orbit, a_mb, is computed from the Moon–Earth semi-major axis, a_mp = 2.5696 × 10⁻³ AU (Williams 2017), as (15)

with m_p and m_m the masses of the Earth and Moon, respectively. Because of the edge-on geometry, the right ascensions of the ascending nodes of the orbits of the barycenter and the moon are set to zero. Because both orbits are assumed to be circular, their perihelions are undefined. The barycenter’s argument of perihelion is chosen precisely behind the star at time t = t₀ = 0, that is, ω_b = 270°. For the moon, ω_m is set to zero. The observational and orbital geometry at t = 0 is sketched in Fig. 4. We use R_p = 6371.0 km and R_m = 1737.4 km for the baseline radii of the planet and the moon, respectively.

3 Numerical results

Here, we present the computed total flux F, the linearly polarized fluxes Q and U, and the degree of polarization P of starlight that is reflected by our model planet–moon system across time. As a trade-off between spatial resolution, radiometric and polarimetric accuracy, and computational time, we use 50 and 14 pixels along the equators of the planet and moon, respectively (see Appendix C), resulting in N_p = 1956 and N_m= 156 (Eq. (6)). In Sect. 3.1, we analyze the individual contributions of the planet and the moon, and in Sect. 3.2, the results for the spatially unresolved planet–moon system. In Sect. 3.3, we take a closer look at particular transit and eclipse events.

Fig. 4 The orbital geometry of our edge-on planet–moon system at time t = 0 as seen from the positive y-axis (left panel) and from the observer’s position at the positive z-axis (right panel). Indexes s, p, m, and b refer to the positions of the star, the planet, the moon, and the planet–moon system barycenter, respectively. Distances and radii are not to scale.

3.1 Reflection by the spatially resolved planet and moon

In order to understand the traces of eclipses and transits in the flux and polarization of starlight reflected by spatially unresolved planet–moon systems, we first discuss the disk–resolved signals of the planet and the moon separately. Figure 5 shows the elements of the locally reflected Stokes vectors F^p and F^m and the direction of polarization χ, with respect to the planetary and lunar scattering planes, respectively, at phase angles α of 0° and 50°.

At α =0° (Fig. 5a), both the planet and the moon would be behind the star and thus invisible, but their disk-resolved signals give insight into the reflection processes. For both bodies, total flux F is maximum at the sub-stellar/sub-observer region and decreases toward the terminator (which coincides with the limb at this phase angle). Because of the Lambertian reflection of the lunar surface and the lack of atmosphere around the moon, the reflected flux is unpolarized and χ is undefined for the moon. The linearly polarized fluxes Q and U of the planet are due to Rayleigh scattering in the planet’s atmosphere. At the sub-stellar region, both Q and U are zero because of symmetry. The general increase of Q and U toward the limb is due to polarized second order scattered light, which is also apparent from the direction of polarization χ. Because of its definition, Q(U) equals zero along the lines at angles of 45° (0°) and − 45° (90°) with the horizontal. Integrated across the planetary disk, P would equal zero. We note that because of the Lambertian reflection of the surface of the planet, Q and U are independent of planetary surface albedo a_surf, while P will generally decrease with increasing a_surf because of the increasing flux F (see, e.g., Stam 2008, for sample computations).

At α = 50° (Fig. 5b), the total flux F of the moon is maximum at the sub-stellar region and decreases toward the terminator, due to the isotropic surface reflection and the absence of an atmosphere. The planet also shows a decrease of F toward the terminator, but the location of the flux maximum is more diffuse and more toward the limb than on the moon, because light that is incident on the planet is scattered in the atmosphere in addition to being reflected by the surface; the reflected flux thus also depends on the optical path lengths through the atmosphere,which in turn depend on the local illumination and viewing angles. The planet’s polarized fluxes Q and U and angle χ are mostly determined by starlight that has been singly scattered by the atmospheric gas molecules. For our choice of reference plane, U is antisymmetric with respect to this plane (and U would thus equal zero when integrated across the disk), and Q is symmetric. The negative values for Q in Fig. 5b indicate that the reflected light is polarized perpendicular to the reference plane, which is indeed also clear from the polarization angle χ, and what is expected for a Rayleigh scattering atmosphere.

3.2 Reflection by the spatially unresolved planet and moon

Figure 7 shows the disk-integrated F, Q, U, and P as functions of the planetary phase angle α_p. For our edge-on system, the curves cover only half of the barycenter’s orbital period. For comparison, we have also included curves for the planet and the moon as isolated bodies, thus without any mutual events. The total flux of the planet–moon system is lower than the sum of the total fluxes of the isolated planet and moon because the latter have not been scaled to the actual radii of the planet and the moon. Indeed, the moon’s flux at α_p = 0° equals the moon’s geometric albedo, 0.067, which matches the theoretical geometric albedo of a Lambertian reflecting body with surfacealbedo of 0.1 (see Stam et al. 2006). The geometric albedo of the unresolved planet–moon system (with both bodies at α_p = 0° and next to each other) is about 0.33. In Fig. 7, the observable planet–moon flux at α_p = 0° is zero, because both bodies are then located behind the star as seen from the observer (in addition, the moon is located behind the planet, as can be seen from Fig. 4, and from situation 1 in Fig. 6).

The curves for F decrease smoothly with increasing α_p, apart from the occasional sharp dips due to eclipses and transits (to be discussed below) and reach zero close to α_p = 180°, where the planet and moon would both be in front of the star. The slightly different slope of the lunar flux phase function as compared tothat of the planet is due to the scattering of light in the latter’s atmosphere. As can be seen in Fig. 7, without the sharp dips, the smooth flux phase function of the planet–moon system does not reveal the presence of a moon, especially not without accurate information on the planetary radius, orbital distance, and atmospheric and surface properties.

Because the lunar surface is completely depolarizing, the moon’s polarized fluxes Q and U are zero at each α_p. The disk-integrated U of the light reflected by the planet is zero due to symmetry (see Fig. 5). Polarized flux Q and degree of polarization P of this light both show a smooth dependence on α_p, apart from the occasional sharp peaks that will be discussed below. The degree of polarization is maximum at phase angles between 90° and 100° due to the atmospheric Rayleigh scattering. At about 165°, the direction of polarization changes from perpendicular (χ = 90°) to parallel (χ = 0°) to the reference plane, and P equals zero. The degree of polarization of the unresolved planet–moon system is somewhat lower than that of the isolated planet, because of the added unpolarized lunar flux. If the planetary atmosphere contained clouds, the shape of this continuum curve would depend on the optical thickness and altitude of the clouds, the microphysical properties of the particles, and the cloud coverage across the planetary disk (for sample curves, see Karalidi et al. 2012; Rossi & Stam 2017, and references therein).

While the smooth curves for the spatially unresolved planet–moon system shown in Fig. 7 do not provide direct evidence of the presence of a moon, the mutual events result in a series of dips and peaks in the reflected flux and polarization, respectively. Figure 6 illustrates the various events. Both the planet and the moon are initially (α_p = 0°) behind the star (position 1 in Fig. 6). Given the prograde lunar motion, the next event, when planet and moon are in view of the observer, is a lunar transit (position 2) and an eclipse of the star on the planet (3). After the first lunar period, the moon again disappears behind the planet (4), followed by an eclipse of the star on the moon (5). This sequence repeats along the barycenter’s orbit.

Both the planetary and lunar eclipses and transits temporarily reduce the flux F that the observer receives. Indeed, when the planet transits the moon, the system’s flux phase function equals that of the isolated planet. The dip in the system’s flux due to a lunar transit (moon in front of the planet) will depend on the radius of the moon as compared to that of the planet and on the lunar surface albedo: the lower the lunar surface albedo and/or the larger the lunar radius, the deeper the dip compared to the continuum. The depth of the dip in the system’s flux F due to an eclipse depends on the relative sizes of the moon and the planet, the reflection properties of the eclipsed body, and on the precise orbital geometry, especially because an eclipse shadow on the moon will not be completely black (and the total flux F thus slightly larger; cf. Fig. 2) due to starlight that is refracted through the limb of the planetary atmosphere and reaches the moon. This refraction is not included in our code (due to the wavelength dependence of Rayleigh scattering, the contribution of refracted light would be larger in the (near) infrared region of the spectrum than at 450 nm).

Because the moon reflects unpolarized light, neither a planetary transit (planet in front of the moon) nor an eclipse on the moon leads to a reduction of the polarized fluxes, as can be seen in Fig. 7. Because the planet reflects polarized light, a transit of the moon and an eclipse on the planet will both decrease Q (given the geometry of our system). Because P depends on F and Q, the dips in F due to less (unpolarized) lunar light being observed yield peaks in P. The peak value of P that is due to the planet transiting the moon equals P of the isolated planet at that value of α_p. In our computational results, peaks in P that are due to an eclipse on the moon would equal P of the isolated planet when the whole lunar disk would be in the planet’s umbra, because we neglect refracted starlight through the limb of the planetary atmosphere. Changes in P that are dueto the moon transiting the planet or due to the moon casting a shadow on the planet will depend on the total and polarized fluxes of the region of the planetary disk that is covered or darkened, and thus, for a given model planet and its atmosphere,on the relative sizes of the moon and the planet and the precise orbital geometry. This will be discussed further in Sect. 3.3.

The absolute depth of the dips in F and Q decreases with increasing α_p because the fraction of a body’s disk that is illuminated and visible decreases with increasing α_p. The amplitudes of features in P for our planet-moon model system are maximum when α_p ≈ 90°. This is particularly convenient for exomoon detection with direct imaging techniques, because that is the phase angle range where the angular distance between the planet–moon system and the parent star will be largest.

As can be seen in Fig. 7, the phase angle gap between a lunar transit and the subsequent planetary eclipse (or a planetary transit and a subsequent lunar eclipse) increases with increasing α_p. As the orbital speed of both bodies is constant in our baseline system with circular orbits, this also applies in the time domain. Indeed, the lunar and planetary transits have a characteristic period because an observer–planet–moon alignment occurs twice per lunar orbit (see Fig. 6). The time gap between two consecutive transits and eclipses, however, increases with increasing α_p because of the movement of the barycenter along its orbit.

Fig. 5 Reflected fluxes F, Q, and U, and the direction of polarization χ across the planet and the moon at α=0° (panel a) and 50° (panel b). Fluxes Q and U, and χ are defined with respect to the scattering planes of the planet and the moon, respectively. In order to facilitate the interpretation of the direction of polarization, we plot 180° − χ. Fluxes Q and U are zero across the moon’s disk because of the Lambertian reflection. All fluxes have been normalized such that the disk-integrated flux F at α = 0° equals the body’s geometric albedo. Absolute planetary fluxes per pixel are not comparable to their absolute lunar fluxes per pixel because of the different number of pixels across each disk.

Fig. 6 Sketch illustrating the sequence of planetary (1, 4, 8, …) and lunar (2, 6, …) transits, as well as planetary (3, 7, …) and lunar (1, 5, 9, …) eclipses for part of the barycenter’s orbit for an edge-on system. The positive z0 -axis points toward the observer. Position 1 corresponds to phase angle α = 0° and time t =0 s in our simulation (cf. Fig. 7).

Fig. 7 Total flux F, the linearly polarized fluxes Q and U, and the degree of polarization P of the spatially unresolved, edge-on, baseline planet–moon system, as functions of the planet’s phase angle αp. Also included are curves for the isolated planet (coinciding with Q and U of the planet–moon system) and the isolated moon (equal to zero in Q, U, and P). Fluxes have been normalized such that at αp= 0°, F equals the geometric albedo of the planet–moon system or each of the isolated bodies. The labels in the plot for F refer to the illustrations in Fig. 6.

3.3 Analysis of the mutual events

In this section, we analyze individual mutual events, their shape, symmetry, periodicity, magnitude, and duration. For this analysis, the change in flux F and degree of polarization P during an event are defined as

First, we will discuss the lunar and planetary transits, then the lunar and planetary eclipses.

Fig. 8 Sketch of the ingress and egress of the moon during a lunar transit for αp = 70° and αp =110°. The arrow indicates the direction of motion of the moon across the planetary disk.

3.3.1 Lunar transits

Figures 9a and b provide detailed views of Δ F and Δ P during the six lunar transits (moon in front of the planet) shown in Fig. 7, together with sketches of the geometry of the planet and the moon at the beginning and the end of the transit for α_p ≈ 80°. As expectedwith constant orbital speeds, the duration of a lunar transit event decreases with increasing α_p because of the decrease of the illuminated area on the planetary disk, and thus the shift of the time of ingress (see Fig. 8). Because egress takes place over the planetary limb, all curves in Figs. 9a and b have the same egress time. Also, the planet is relatively dark near the terminator, and thus yields a smooth flux decrease upon the lunar ingress, while it is bright near the limb (see Fig. 5), yielding a rapid increase of F upon the lunar egress.

The depth ΔF depends strongly on α_p, because with increasing α_p, the illuminated area, and hence also the covered area on the planet, decreases. The shape of Δ F also depends on α_p. At α_p = 0°, the curve would be symmetric. At larger values of α_p, the trace of the lunar nightside starts to appear in the curve. Because of the moon’s prograde orbit, the lunar dayside ingresses before the nightside. In the curves for α_p = 13.4° and 40.3°, the steeper decrease of F due to the ingress of the lunar nightside can be seen. The value of ΔF that is reached within a transit at a given α_p depends on the lunar albedo and on the area of the planetary disk that is covered, thus on the lunar radius. The lower the lunar surface albedo and/or the larger the lunar radius as compared to the planetary radius, the larger Δ F will be. As an example, Fig. 11a shows F at α_p = 67.2° for various values of the lunar radius expressed as fraction of the planetary radius (the value for the baseline model is approximately 0.3). As can be seen, the continuum flux increases with increasing lunar radius due to the increased amount of flux reflected by the moon, and indeed the lowest flux during the transit decreases and Δ F increases with increasing lunar radius.

We note that a change in the lunar radius implies a change in the lunar mass (assuming a similar composition) and, thus, a change in the lunar period around the planet. While the frequency of the events decreases nonlinearly with increasing lunar radii, we have aligned the mutual events in Fig. 11 in time to facilitate a comparison. Mutual transits show up every half lunar sidereal period. Because , relative timing variations of 1–13% are obtained for lunar–to–planet radius ratios from 0.1 to 0.7 (with our baseline value of approximately 0.3). In the case of eclipses, and assuming coplanar circular lunar and planetary orbits, the repetition period equals the lunar synodic period, for which timing variations of 1–14% are obtained for the same range of radii ratios.

Figure 9b shows ΔP during lunar transits. It can be seen that P can also decrease during a transit, which is not apparent from the curves in Fig. 7. The curves exhibit a strong variation in shapes, and with increasing α_p become increasingly asymmetric. The largest ΔP is found around α_p = 90°, where P of our model planet is highest (see Fig. 7). The precise shapes of the curves depend on the properties of the planet and its moon and the path of the transit across the planetary disk.

In our planet–moon system, the lunar transit occurs along the planet’s equator, where the antisymmetry of U yields a null net contribution, and the shape of ΔP thus depends on the variation of Q and F along the path (cf. Fig. 5). At α_p = 13.4°, Q is maximum near the planet’s limb and close to zero at the center of the body. The disk-integrated P is close to zero due to symmetry. During ingress and egress, the moon breaks the symmetry, and (slightly) increases the disk-integrated value of P (see Fig. 9b). With increasing α_p, the maximum ΔP increases to reach a maximum (in the figure) at α_p = 90.1°. The α_p where the maximum P is found corresponds roughly with the α_p of the minimum F. This increase in P appears to be driven by the decrease of F.

The negative values of ΔP in the curves for α_p = 67.2°, 90.1°, and 120.9° indicate that during that part of the transit, the decrease in Q is larger than that in F. This happens in particular when the illuminated part of the moon is transiting the illuminated part of the planetary disk, while the dark part of the moon is still transiting the dark part of the planet, and when the moon transits the limb of the planet, where Q is relatively large and where the transit thus strongly decreases Q. As mentioned earlier, the maximum ΔP (on the order of 2% for our planet–moon system) should be observable when α_p is between about 70° and 150°, where the angular distance between the planet–moon and the star is relatively large, which would facilitate the detection of lunar transits with direct detection methods.

Figure 11b shows the change in P at α_p = 67.2°, for various values of r, the lunar radius expressed as a fraction of the planetary radius. With increasing r, the continuum P decreases because more unpolarized flux reflected by the moon is added to the total flux. During a transit, P can be seento be very sensitive to the lunar radius. Indeed, with increasing r, Δ P changes from positive (when P during the event is higher than in the continuum) to negative (when P during the event is lower than in the continuum) because apparently the polarized flux reflected by the planet decreases more during the event than the total flux reflected by the planet with the moon in front of it.

Fig. 9 Changes in the total reflected flux ΔF (panels a and c), and the degree of polarization ΔP (panels b and d), as functions of the lunar true anomaly, νmb, and relative time for the lunar transits (top panels) and planetary transits (bottom panels) shown in Fig. 7. The time step of these simulations is 3 min.

3.3.2 Planetary transits

Figures 9c and d show the planetary transits (planet in front of the moon) from Fig. 7 in detail, for both F and P. At α_p = 0°, the planet and the moon are behind the star and would thus not be visible to the observer, but, similar to the lunar transits, planetary transits would yield symmetric events in ΔF and ΔP. With increasing α_p, the events become more asymmetric.

With a planetary transit, the shapes of the ΔF curves are very different from those of the lunar transits, because, for our orbital geometry, the moon will be completely covered during the planetary transit and while it is covered, the transit curve is flat. The depth of the transit depends on the total amount of reflected flux received from the (isolated) moon. For a given value of α_p, Δ F will increase linearly with the surface area of the lunar disk (thus, with the lunar radius squared), and/or with the lunar surface albedo. The start time of the transit depends on α_p, as α_p determines the extent of the illuminated area on the moon, and thus when it will be covered. Like with the lunar transits, the end of the transit, over the bright limb of the moon, is independent of α_p; it only depends on the lunar true anomaly.

For the planetary transits, ΔP is entirely due to the decrease in F, as the moon itself reflects only unpolarized light, and the planet thus blocks no polarized flux. As a result, the transits in Δ P are flat as long as the illuminated part of the lunar disk is covered. The maximum Δ P depends on both Δ F and the planet’s polarized flux Q, and thus on α_p for a given planet–moon model. In Fig. 9d, a maximum ΔP of about 2.5% occurs at α_p =80.6°. Through ΔF, Δ P will increase with the lunar surface albedo and/or the surface area of the lunar disk at a given α_p and for a given planet–moon model; a darker and/or smaller moon would yield a smaller Δ F and hence a smaller ΔP.

Fig. 10 Similar to Fig. 9, except for the planetary eclipses (top panels) and the lunar eclipses (bottom panels), both as functions of angle φms (see Fig. 12).

3.3.3 Planetary eclipses

Figures 10a and b show the curves of ΔF and ΔP during the planetary eclipse events shown before in Fig. 7. During these eclipses, the moon casts its shadow on the planet, and because in our geometry the lunar orbital plane coincides with the barycenter’s orbital plane, the shadow of the moon travels along the horizontal line crossing the center of the planetary disk. Figure 12 illustrates the geometries for the planetary eclipse, with φ_ms the angle between the star and the moon measured positive in the counterclockwise direction from the center of the planet. Angle φ_ms is used as a relative measure of eclipse events. Its relation with time is linear given the circular orbital motion of the bodies.

For planetary eclipses, the explanation regarding the asymmetry of Δ F and Δ P is the same as for lunar transits, with one important difference: while the transit events start increasingly later with increasing α_p and end at the same (relative) time, the planetary eclipses start at the same (relative) time and end increasingly earlier with increasing α_p. The reason for this difference is that eclipses depend on the position of the star with respect to the planet–moon system, and not on the position of the observer. The observer’s position does influence the fraction of the eclipse that is captured, as it determines the phase angle and hence the fraction of the illuminated part of the planetary disk across which the eclipse travels. Thus, with increasing α_p, the duration of an eclipse decreases, as is visible in Figs. 10a and b. The depth Δ F decreases with increasing α_p because less of the illuminated part of the planetary disk is visible.

The shape of the ΔF curves for the planetary eclipses (where the moon’s shadow moves across the planetary disk) appears to be more gradual than that of the lunar transit curves (where the moon itself moves across the planetary disk; cf. Fig. 9). This is because the planet first travels through the lunar penumbral shadow before entering the deep, umbral shadow cone. Because we discuss only half of the barycenter’s orbit, the ingress of the lunar eclipse shadow on the illuminated part of the planetary disk is through the terminator and the egress through the bright limb, as illustrated in Figs. 8 and 12. However, because the spatial extent of the penumbral and umbral shadows across the disk are smallest halfway through the total duration of the eclipse (as seen from a vantage point on the moon facing the planet), the egress of the lunar shadow yields a much smoother Δ F curve than observed during the egress of a lunar transit. The difference in the maximum Δ F during a lunar transit and a planetary eclipse is most apparent at the smaller phase angles, because with increasing α_p, the contribution of light reflected by the moon decreases. We note that, differently than for planetary transits, the value of Δ F during a planetary eclipse is independent of the lunar surface albedo. It will obviously increase with the radius of the moon relative to that of the planet. This can also be seen in Fig. 11c which shows F at α_p = 72.5° for various values of r, the lunar radius expressed as fraction of the planetary radius (for the baseline model, r is 0.3). Indeed, with increasing r, the continuum flux increases because of the added flux reflected by the moon, and the minimum flux during the event decreases because of the increasing extent of the lunar shadow.

The change in P during the planetary eclipses is shown in Fig. 10b. The increase in P during the ingress of the lunar shadow is due to the decrease in F and a decrease in Q. As the shadow progresses across the disk, its spatial extent decreases, and its influence on P decreases. The maximum of ΔP appears to be 1−2% for phase angles α_p ≈ 70°−160°. For the largest phase angles, the corresponding value of ΔF is relatively small, because only a narrow crescent of the planet is illuminated, so there Δ P is mostly due to a change in Q.

Figure 11d shows P at α_p = 72.5° for various values of r (the lunar radius expressed as a fraction of the planetary radius). With increasing r, the continuum P decreases because of the added unpolarized flux reflected by the moon, and P during the eclipse also decreases, apparently because the polarized flux Q decreases more than the total flux F.

Fig. 11 Changes in the total reflected flux ΔF (panels a and c) and ΔP (panels b and d) during lunar transits at α = 67.2° (top panels) and planetary eclipses at α = 72.5° (bottom panels) for various lunar-to-planetary radius ratios r. The horizontal axis shows the elapsed time since the concentric alignment of the planet and moon as seen from the star in the case of an eclipse and as seen from the observer in the case of a transit. The time step of these simulations is 3 min. The baseline lunar-to-planetary radius ratio r is about 0.3.

3.3.4 Lunar eclipses

As can be seen in Figs. 10c and d, lunar eclipses, when the planet casts its shadow on the moon, show a similar symmetry as the planetary transits, where the planet moves in front of the moon (Figs. 9c and d). Because our model moon is small compared to the planetary shadow, both the Δ F and Δ P curves are flat except during ingress and egress. The flux changes during ingress and egress of the lunar eclipse, respectively,are smoother than with the planetary transits, due to the extended penumbral region of the planet’s shadow.

As with the planetary eclipses discussed above, the moon’s ingress into the planetary shadow occurs at the same value of angle φ_ms (cf. Fig. 12) independent of phase angle α_p. The duration of the eclipse decreases with increasing α_p because of the decrease of the illuminated area on the moon with increasing α_p. The change in P during lunar eclipses, shown in Fig. 10, is similar to that during planetary transits (Fig. 9), as in both cases P changes during the event because F decreases and because there is no actual change in the amount of polarized flux from the system, because our model moon reflects only unpolarized light. The maximum ΔP value is about 2.7%, attained at α_p ≈ 87.1° in Fig. 10.

With increasing lunar radius as compared to the planetary radius, and depending on the distance between the moon and the planet and the distance to the star, the planetary shadow might cover only part of the lunar disk. In that case, F and P will no longer be flat during the eclipse, because there will be a contribution of unpolarized lunar flux that will vary in time. Indeed, the curves will become asymmetric (more similar to those for the planetary eclipses). Because the moon only reflects unpolarized flux, P will always increase during the eclipses.

4 Conclusions

We present numerically simulated flux and polarization phase functions of starlight that is reflected by an orbiting planet–moon system, including mutual events such as transits and eclipses. Most results presented in this paper apply to a Moon-sized, Lambertian (i.e., isotropically and depolarizing) reflecting moon orbiting an Earth-sized exoplanet with an Earth-like, gaseous atmosphere on top of a Lambertian reflecting surface (the surface pressure is 1 bar), in an edge-on configuration. Our results show that the flux and polarization phase functions of starlight reflected by such a planet–moon system contain traces of the moon in the form of periodic changes in the total flux F and degree of polarization P as the bodies shadow each other (eclipses) and/or hide one another from the observer’s view (transits) along their orbit around the star. These changes in F and P are only one order of magnitude smaller than the system’s continuum phase functions. The magnitude, shape, and duration of the obtained total flux signatures are comparable with the results by Cabrera & Schneider (2007), except that they do not include the influence of the penumbra.

During events that darken the planet, that is, lunar transits and planetary eclipses, the shape of the dip in F depends on the reflection properties of the regions on the planet along the path of the shadow. The change in P during such events strongly depends on the ratio of polarized–to–total reflected flux across the disk and along the path of the shadow. Indeed, Δ P ≈ 1−2.5% for 67° < α_p < 121°. For the planet–moon system used in this paper, we found the strongest changes in P during either the first (planetary eclipse) or the last (lunar transit) half of the event, compared to the duration of the event in F, in particular at intermediate phase angles. The asymmetry of the planet darkening events as imprinted on the change in polarization Δ P is due to the variation of the polarization across the planetary disk with our model atmosphere–surface: the polarized flux at the limb is higher than at the terminator.

During lunar darkening events, that is, planetary transits and lunar eclipses, the size difference between the planet and the moon yields a relatively symmetric change in F and, due to the nonpolarizing lunar reflection, a similarly relatively symmetric change in P, as the latter is only due to the decrease in total flux, upon a lunar darkening event, not to a change in polarized flux. The curves for planetary transits have steeper slopes during the ingress and egress phases than the curves for the lunar eclipses, because with the latter the moon travels through the penumbral shadow of the planet. The change in P depends on the size and albedo of the moon, and on the polarization signal of the planet, which itself depends on the atmosphere–surface model and the phase angle. Our simulations have been performed at 450 nm (in the blue) where the scattering by the gas in the Earth-like planetary atmosphere strongly contributes to the planet’s polarization signals. In particular, at intermediate phase angles, the polarization signal of a gaseous atmosphere is strong, and the change in polarization during the lunar darkening events can reach a few percent. Indeed, Δ P ≈ 1.25−2.66% for 54° < α_p < 108° during planetary transits. At these phase angles, the angular distance between the planet–moon system and the parent star is relatively large, so these angles are favorable for the detection of reflected starlight. For a planet with clouds in its atmosphere, the continuum flux phase function will have a similar shape as that for our cloud-free planet, except the total amount of reflected flux will be larger (of course not at wavelengths where atmospheric gases absorb the light). The polarization curve of a cloudy planet will show angular features due to the scattering of light by the cloud particles, such as the rainbow for liquid water droplets (see, e.g., Bailey 2007; Karalidi et al. 2012, and references therein).

The duration of a transit event depends on the orbital parameters, on the sizes of the planet and moon, and on the phase angle (the latter mostly for lunar transits). In our simulations, a typical planetary transit takes approximately 4 h both in flux and polarization. A lunar transit at an intermediate phase angle of 90°, takes about 2 h in flux. In polarization, the change in P is apparent during a shorter period than the change in flux F, due to the distribution of polarized flux across the planetary disk. The duration of eclipse events is somewhat longer than that of transit events due to the diverging shape of the shadow cone. In our simulations, eclipse events can take up to 6 h, where the polarization change in the planetary eclipses is only apparent during part of the time of the flux change.

The results presented in this paper correspond to half of the planetary orbit around the star. The results for the other half of the orbit will be similar, except that the curves will be mirrored with respect to the central event time, because transit and eclipse ingresses and egresses will happen over the other side of the darkened body. Our results show that measuring the temporal variation in F and/or P during transits and eclipses could provide extra information on the properties of a planet and/or moon and their orbits. Extracting such information, however, requires not only detecting such events but also measuring the shape of the variations in F and/or P. For the interpretation of such measurements, numerical simulations to map in more detail the influence of the physical characteristics of the moon and the planet (radius, albedo, atmosphere-surface properties) and their orbital characteristics (inclination angles, ellipticity) on the temporal variation in F and P are required. Such simulations will be targeted in future research.

Fig. 12 Geometrical definition of angle φms as moon m passes between planet p and star s. Seen from the top, the moon moves counterclockwise around the planet. A similar definition holds for angle φmp as the planet passes between the moon and the star. Distances between bodies and radii are not to scale.

Appendix A Local illumination and viewing angles

In this appendix, we describe the computation of the angles required to compute the starlight that is reflected by each pixel on the planet (cf. Eq. (8)): the phase angle α, the local viewing zenith angle θ_i, the local illumination zenith angle θ_0i, the local azimuthal difference angle ϕ_i − ϕ_0i, and the local rotation angle β_i. To add the computed Stokes vectors of the moon to those of the planet, we usually also need rotation angle ψ that redefines the lunar Stokes vector from the lunar scattering plane to the planetary scattering plane (cf. Eq. (4)).

A.1 Phase angle α_x

Phase angle α is the angle between the direction to the star and the observer as measured from the center of a body (see Fig. A.1). In principle, its value ranges from 0° to 180°, although the phase angle range accessible to an observer depends on the inclination angle of the orbit of a body. A body in an edge-on orbital geometry (inclination angle i = 90°) can attain phase angles between 0° (when it is located behind the star) and 180°, while a body in a face-on orbital geometry (i = 0°) can only be observed at α = 90°. Generally, given an orbital inclination angle i, the phase angle range is given by (A.1)

The phase angle of the planet or the moon at time t is computed as (A.2)

where subscript “x” refers to either “p” (planet) or “m” (moon), is the unit vector along the z-axis, pointing toward the observer, and vector r_xs connects thecenter of the planet or moon with the center of the star. Given the small separation between the planet and moon compared to their distances to the star, α_p is virtually thesame as α_m, yet our numerical model uses both values.

A.2 Local viewing zenith angle θ_i

The local viewing zenith angle θ_i is the angle between the zenith direction of pixel i and the direction toward the observer (see Fig. A.1). Angle θ_i takes values between 0° (at the sub-observer location) and 90° (at the limb). It depends on the location of the pixel on the disk of the planet or moon and is thus time-independent. It is computed according to (A.3)

where subscript “x” refers to either “p” or “m”, u_z is the unit vector along the z-axis that points toward the observer, and r_ix is the vector pointing to the center of the pixel from the center of either the planet or the moon.

A.3 Local illumination zenith angle θ_0i

The local illumination zenith angle θ_0i is defined as the angle between the local zenith direction of pixel i and the direction toward the star (see Fig. A.1). Angle θ_0i takes values between 0° (at the sub-stellar location) and 90° (at the terminator). The position of the star changes in time, so that the time-dependent local illumination zenith angle can be computed as (A.4)

where subscript “x” refers to either “p” or “m”, r_ix is the vector from the center of the planet or moon to the center of the pixel, and r_is is the vector from the center of the star to the center of the pixel on the planet or moon.

Fig. A.1 Angular geometry of the spherical triangle centered at pixel i and defined by the zenith direction unit vector, ui, the observer’s direction unit vector, uob, and the star direction unit vector, us. The sides of the spherical triangle are: the observer–zenith angle θi, the star–zenith angle θ0i, and the pixel-based phase angle αi, all centered at pixel i. The angle between sides θ0i and θi is the azimuthal difference angle ϕi − ϕ0i.

A.4 Local azimuthal difference angle ϕ_i − ϕ_0i

The azimuthal difference angle ϕ_i − ϕ_0i for pixel i on the planet or moon is the angle between the plane described by the local zenith direction and the direction toward the observer and the plane described by the local zenith direction and the direction toward the star⁶. As Fig. A.1 shows, ϕ_i − ϕ_0i follows from

where α_i is the angle between the direction to the observer and the direction to the star measured from the center of pixel i. Given that with “x” referring to either “p” or “m”, α_i can be approximated by the body’s phase angle α_x, and thus (A.5)

A.5 Local rotation angle β_i

Angle β_i is used to rotate a locally computed vector (see Eq. (6)) for pixel i on the planet or the moon from the local meridian plane to the scattering plane of the body, which is used as the reference plane for the disk-integrated signal of the body. The pixel grid across the planet is defined with respect to the planetary scattering plane, and β_i is thus time-independent for the planetary pixels. For a pixel i, β_i is computed according to (A.6)

where x_ip and y_ip are the coordinates of the center of the pixel (recall that the z-axis points toward the observer).

For the lunar pixels, the alignment between the lunar scattering plane and the lunar grid, and hence angle β_i, is time-dependent and requires us to redefine the pixel coordinates with respect to the lunar scattering plane. Indicating the redefined coordinates of lunar pixel i with subscript j, angle β_j is then computed as (A.7)

A.6 Scattering plane rotation angle ψ

Scattering plane rotation angle ψ is used to rotate a Stokes vector that is defined with respect to the lunar scattering plane to the planetary scattering plane, which we use as the reference plane for the planet–moon system. Angle ψ is measured in a clockwise direction from the lunar scattering plane to the planetary scattering plane. For the results presented in this paper, the moon and the planet orbit in the same, edge-on plane, and angle ψ equals zero. In the general case, however, it is computed using (A.8)

where r_ms is the vector from the star to the centerof the moon, and u_x and u_y are the unit vectors along the x-axis and y-axis in coordinate system S₁.

Appendix B Computing eclipses

An eclipse occurs when body A is between the star S and body B such that the shadow of A falls onto B. The effect of a planetary or lunar eclipse depends on the positions of the star, moon, and planet, and, due to the extended size of the star, thesize, shape, and depth of the shadow depend not only on the radii of the star and the eclipsing body, but also on the distances and angles involved. Computing eclipses has been discussed in great detail by Link (1969) for the Moon–Earth system, which we apply to our exoplanetary system. We model the umbral, antumbral, and penumbral shadow regions. Figure B.1 shows the geometries involved in the various types of eclipses. The equations used for computing the influence of eclipses are described here.

The fluxarriving at a pixel i of eclipsed body B at time t depends on the fraction of the stellardisk and the local stellar surface brightness, as seen from the center of the pixel. In Eq. (6), this is accounted for by factor c_i, the ratio between the actual flux on pixel i and e_i, the flux on the non-eclipsed pixel, (B.1)

with and S_Si being the stellar disk area as observed from pixel i when it is eclipsed and when it is non-eclipsed, respectively. We ignore stellar limb darkening and stellar light that travels through the atmosphere of the eclipsing body (if present).

To determine c_i, we first have to identify whether or not pixel i is eclipsed. Obviously, c_i = 1 for a non-eclipsed pixel. If the pixel is (partly) eclipsed, we have to determine the type of eclipse: umbral (i.e., total), penumbral (i.e. partial), or annular³. For an umbral eclipse, c_i = 0.0; for a penumbral and annular eclipse, we have to compute in order to determine c_i.

As can be seen in Fig. B.1, a pixel on body B is eclipsed when it falls within the penumbral cone of body A. Opening angle Ω of the penumbral cone is given by (B.2)

Indeed, pixels in eclipse on the disk of body B can be found at times when (B.3)

with angle ρ () given by (see Fig. B.1): (B.4)

Vector AO₄ is a function of the radii of the shadowed and eclipsing bodies and of Ω, as (B.5)

Except when body B is entirely found inside the penumbral cone, there will also be noneclipsed pixels on the disk. When Eq. (B.3) holds, a pixel-by-pixel search is performed, in which the center of each pixel is checked for total or umbral (Sect. B.1), annular (Sect. B.2), or penumbral (Sect. B.3) eclipse conditions.

B.1 Total or umbral eclipses

In the umbral zone (see Fig. B.1), pixels experience a total stellar eclipse. If the umbral zone is wide and the shadowed body B relatively small, such as in the Earth–Moon system, all pixels on the disk of B can be simultaneously in the umbra, and factor c_i = 0 for all pixels. This is the case when (B.6)

with (B.7)

and (B.8)

The diskof body B will only be partially inside the umbral shadow cone of body A when (B.9)

and (B.10)

Here (B.11)

with (B.12)

If the disk of body B falls partially inside the umbral cone, the pixels where β_i < Ψ are inside the umbra, and c_i = 0. Because c_i only applies to pixels on the illuminated part of the disk of B, this condition can be reformulated as (B.13)

Fig. B.1 Geometry of the umbral, antumbral, and penumbral shadow cones, when star S is eclipsed by body A, casting a shadow on body B. The shadow cast by A into space is rotationally symmetric around the axis through the center of the star and body A. The radii of the star, body A, and body B are denoted by RS, RA, and RB, respectively. Points O1, O2, O3, O4, and O5 denote auxiliary points: the umbral and antumbral cones have apex O1 and aperture2Ψ and the penumbral cone has apex O2 and aperture2Ω. The lower figure also shows angles ζ, ρ, ωi, βi, and θ, that are used in the computation of the eclipse shadow depth. Distances between bodies and radii are not to scale in order to emphasize the geometry of the system.

We note that angle β_i of a pixel can be derived from (B.14)

The pixels on the disk of B that are not in the umbral cone can be in the penumbral or annular eclipse zone, as described below.

B.2 Annular or antumbral eclipses

Instead of crossing the umbral cone, body B can cross the antumbral cone, where eclipsing body A does not completely cover the stellar disk as seen from body B, thus yielding a so-called annular eclipse. Body B is in the antumbral shadow cone when (B.15)

where ξ and Ψ follow from Eqs. (B.8) and (B.7). When Eq. (B.15) is satisfied, pixels on the disk of B are checked for their eclipsed status. A pixel is (partially) eclipsed if (B.16)

For each pixel in the antumbral cone, factor c_i is given by the fraction of the stellar disk that is visible (see Fig. B.2), that is, (B.17)

where (B.18)

Fig. B.2 Disks of star S and eclipsing body A seen from a pixel on body B in the antumbral zone. Angles αS and αA indicate the angular radii of the bodies.

Here, r_Si is the position vector of pixel i on body B with respect to the center of star S.

B.3 Penumbral eclipses

When Eq. (B.3) holds, all pixels that are not in the umbral or antumbral eclipse regions are examined for being in the penumbral shadow. Indeed, pixels with ω_i < Ω are within the penumbral cone, as can be seen in Fig. B.1. This inequality can be rewritten as (B.19)

where cosΨ follows from Eq. (B.7) and cosω_i is given by (B.20)

The magnitude of the eclipse at pixel i, c_i, can be calculated through the stellar and eclipsing body viewing angles, that is, the angular diameter of the bodies, 2α_S and 2α_A, and the eclipsing body-to-star angular distance as seen from the shadowed body, δ. Then, as follows from Fig. B.3, c_i can be computed using (B.21)

with (B.22)

and (B.23)

Here, l_S and l_A are the distances from the centers of the stellar disk and eclipsing body A, respectively, to the line between the two points where the stellar and eclipsing body intersect. The angles θ_S and θ_A (see Fig. B.3) follow from, for example, Heron’s formula: (B.24)

and (B.25)

with (B.26)

Pixels on body B that satisfy Eq. (B.3) but do not meet the conditions established by Eqs. (B.6), (B.10), (B.16), and (B.19) are not eclipsed and have c_i = 1.0.

Fig. B.3 Stellar (S) and eclipsing body (A) disks as seen from a pixel on the shadowed body during a penumbral eclipse. The stellar shadowed area is decomposed in two components A1 and A2. The angular radius of the bodies is represented by αS and αA, and δ is the angular separation between the bodies’ center. The central angles of the circular segments defined by the common cord of the intersecting stellar and eclipsing body disks are represented by θS and θA. The minimum distance from the star and eclipsing body centers to the common chord is defined as lS and lA, respectively.

Appendix C Number of pixels across the disk

The results of our numerical simulations depend on the number of pixels that is used to compute the flux and polarization signals of the disk of the planet and moon. The number of pixels across the disk of the planet or the moon, N^p and N^m, respectively,determines the spatial resolution of the locally reflected Stokes vectors and hence the numerical error in the integration across the disk, in particular at large phase angles. However, the larger the number of pixels, the smaller the error but the longer the computational time. We have investigated the optimal number of pixels, expressed in the number of pixels across the equator of the planet and the moon, and , respectively,using a trade-off between the errors and the computational time, with time steps of 24 h.

For the trade-off, we compute the flux F (Eq. (C.2)) and degree of polarization P (Eq. (C.2)) for consecutive values of N_eq as functions of phase angle α, both for the planet and the moon. In Fig. C.1, we show the maximum and mean differences encountered across the whole phase angle range of the planetand the moon, together with the average disk integration time and the total phase curve computation time. The differences are defined as

with n − 1 and n two consecutive values of N_eq.

Figure C.1 shows that the flux and polarization differences decrease with increasing N_eq, and that for the values of N_eq considered, the computed flux and polarization curves have not yet completely converged. However, further increasing N_eq increases the integration time across the planetary disk, as can be seen in Fig. C.1c. In the simulations presented in this paper, we decided to use and , which yields an average disk-integration time of ~0.8 s (thus an overall phase curve computation time of ~2.4 min with a 24 h temporal resolution). These values for and produce smooth curves for individual transit and eclipse events for temporal resolutions as small as 1 min.

Fig. C.1 Analysis for the number of pixels along the equator Neq of the planet (top panels) and moon (bottom panels). Shown are the maximum (solid line) and mean (dashed line) differences between results computed across the whole phase angle range and for consecutive values of Neq values, for the reflected flux F(α) relative to F(α = 0°) (panels a and d), and degree of polarization P (panels b and e, note that for the moon, P = 0). Also shown is the computational time (in minutes) for the computation of a full phase curve (with 24 h temporal resolution) and the average disk integration (panels c and f). For (top panels), we used values of 10, 20, 30, 40, 50, 60, and 70, and for (bottom panels), 3, 5, 8, 11, 14, 16, and 19.

References

All Tables

All Figures

	Fig. 1 Three-dimensional view (left panel) and projection onto the x0 –y0 plane (right panel) of the discretized planet or moon. The z0-axis points toward the observer. The orientation of the x0 and y0 axes with respect to the disk of the planet or moon can be chosen arbitrarily.
In the text

	Fig. 2 The mutual events between a planet and its moon. Panel a: a lunar transit. Panel b: a planetary transit. Panel c: a planetary eclipse. Panel d: a lunar eclipse. The positive z0 -axis (in panels aand b) points to the observer. The white–black scale indicates full–null illumination of a body. An arrow indicates the lunar motion around the planet. The distances and radii are not to scale.
In the text

Fig. 3 (a) The planet–moon barycenter orbit around the star (b) The moon orbit around the barycenter. The reference frames and angles used to describe the Keplerian orbits of the planet–moon system barycenter around the parent star (panel a) and the Keplerian orbit of the moon around the planet-moon system barycenter (panel b). Plane p1 (panel a) is the plane of the sky as seen by the observer on the positive z1 -axis. Plane p2 (panels a and b) is the barycenter’s orbital plane, and plane p3 (panel b) is the orbital plane of the moon around the barycenter. Angle i is the orbitalinclination angle, ω the argument of periastron, Ω the right ascension of the ascending node, and ν the true anomaly. Subscript bs refers to the barycenter of the planet–moon system around the star, and mb to the moon around the barycenter. Vectors rbs and rmb are the position vectors of the barycenter and the moon, respectively. The barycenter and the moon are indicated by b and m, respectively.

In the text

	Fig. 4 The orbital geometry of our edge-on planet–moon system at time t = 0 as seen from the positive y-axis (left panel) and from the observer’s position at the positive z-axis (right panel). Indexes s, p, m, and b refer to the positions of the star, the planet, the moon, and the planet–moon system barycenter, respectively. Distances and radii are not to scale.
In the text

Fig. 5 Reflected fluxes F, Q, and U, and the direction of polarization χ across the planet and the moon at α=0° (panel a) and 50° (panel b). Fluxes Q and U, and χ are defined with respect to the scattering planes of the planet and the moon, respectively. In order to facilitate the interpretation of the direction of polarization, we plot 180° − χ. Fluxes Q and U are zero across the moon’s disk because of the Lambertian reflection. All fluxes have been normalized such that the disk-integrated flux F at α = 0° equals the body’s geometric albedo. Absolute planetary fluxes per pixel are not comparable to their absolute lunar fluxes per pixel because of the different number of pixels across each disk.

In the text

	Fig. 6 Sketch illustrating the sequence of planetary (1, 4, 8, …) and lunar (2, 6, …) transits, as well as planetary (3, 7, …) and lunar (1, 5, 9, …) eclipses for part of the barycenter’s orbit for an edge-on system. The positive z0 -axis points toward the observer. Position 1 corresponds to phase angle α = 0° and time t =0 s in our simulation (cf. Fig. 7).
In the text

Fig. 7 Total flux F, the linearly polarized fluxes Q and U, and the degree of polarization P of the spatially unresolved, edge-on, baseline planet–moon system, as functions of the planet’s phase angle αp. Also included are curves for the isolated planet (coinciding with Q and U of the planet–moon system) and the isolated moon (equal to zero in Q, U, and P). Fluxes have been normalized such that at αp= 0°, F equals the geometric albedo of the planet–moon system or each of the isolated bodies. The labels in the plot for F refer to the illustrations in Fig. 6.

In the text

	Fig. 8 Sketch of the ingress and egress of the moon during a lunar transit for αp = 70° and αp =110°. The arrow indicates the direction of motion of the moon across the planetary disk.
In the text

	Fig. 9 Changes in the total reflected flux ΔF (panels a and c), and the degree of polarization ΔP (panels b and d), as functions of the lunar true anomaly, νmb, and relative time for the lunar transits (top panels) and planetary transits (bottom panels) shown in Fig. 7. The time step of these simulations is 3 min.
In the text

	Fig. 10 Similar to Fig. 9, except for the planetary eclipses (top panels) and the lunar eclipses (bottom panels), both as functions of angle φms (see Fig. 12).
In the text

Fig. 11 Changes in the total reflected flux ΔF (panels a and c) and ΔP (panels b and d) during lunar transits at α = 67.2° (top panels) and planetary eclipses at α = 72.5° (bottom panels) for various lunar-to-planetary radius ratios r. The horizontal axis shows the elapsed time since the concentric alignment of the planet and moon as seen from the star in the case of an eclipse and as seen from the observer in the case of a transit. The time step of these simulations is 3 min. The baseline lunar-to-planetary radius ratio r is about 0.3.

In the text

	Fig. 12 Geometrical definition of angle φms as moon m passes between planet p and star s. Seen from the top, the moon moves counterclockwise around the planet. A similar definition holds for angle φmp as the planet passes between the moon and the star. Distances between bodies and radii are not to scale.
In the text

Fig. A.1 Angular geometry of the spherical triangle centered at pixel i and defined by the zenith direction unit vector, ui, the observer’s direction unit vector, uob, and the star direction unit vector, us. The sides of the spherical triangle are: the observer–zenith angle θi, the star–zenith angle θ0i, and the pixel-based phase angle αi, all centered at pixel i. The angle between sides θ0i and θi is the azimuthal difference angle ϕi − ϕ0i.

In the text

Fig. B.1 Geometry of the umbral, antumbral, and penumbral shadow cones, when star S is eclipsed by body A, casting a shadow on body B. The shadow cast by A into space is rotationally symmetric around the axis through the center of the star and body A. The radii of the star, body A, and body B are denoted by RS, RA, and RB, respectively. Points O1, O2, O3, O4, and O5 denote auxiliary points: the umbral and antumbral cones have apex O1 and aperture2Ψ and the penumbral cone has apex O2 and aperture2Ω. The lower figure also shows angles ζ, ρ, ωi, βi, and θ, that are used in the computation of the eclipse shadow depth. Distances between bodies and radii are not to scale in order to emphasize the geometry of the system.

In the text

	Fig. B.2 Disks of star S and eclipsing body A seen from a pixel on body B in the antumbral zone. Angles αS and αA indicate the angular radii of the bodies.
In the text

Fig. B.3 Stellar (S) and eclipsing body (A) disks as seen from a pixel on the shadowed body during a penumbral eclipse. The stellar shadowed area is decomposed in two components A1 and A2. The angular radius of the bodies is represented by αS and αA, and δ is the angular separation between the bodies’ center. The central angles of the circular segments defined by the common cord of the intersecting stellar and eclipsing body disks are represented by θS and θA. The minimum distance from the star and eclipsing body centers to the common chord is defined as lS and lA, respectively.

In the text

Fig. C.1 Analysis for the number of pixels along the equator Neq of the planet (top panels) and moon (bottom panels). Shown are the maximum (solid line) and mean (dashed line) differences between results computed across the whole phase angle range and for consecutive values of Neq values, for the reflected flux F(α) relative to F(α = 0°) (panels a and d), and degree of polarization P (panels b and e, note that for the moon, P = 0). Also shown is the computational time (in minutes) for the computation of a full phase curve (with 24 h temporal resolution) and the average disk integration (panels c and f). For (top panels), we used values of 10, 20, 30, 40, 50, 60, and 70, and for (bottom panels), 3, 5, 8, 11, 14, 16, and 19.

In the text