Issue 
A&A
Volume 618, October 2018



Article Number  A162  
Number of page(s)  20  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201833320  
Published online  23 October 2018 
Traces of exomoons in computed flux and polarization phase curves of starlight reflected by exoplanets
^{1}
Faculty of Aerospace Engineering, Delft University of Technology,
Kluyverweg 1,
2629
HS Delft,
The Netherlands
email: D.M.Stam@tudelft.nl
^{2}
GMV AD, Calle de Isaac Newton 11, 28760 Tres Cantos,
Madrid,
Spain
email: j.berzosamolina@gmail.com
Received:
28
April
2018
Accepted:
23
July
2018
Context. Detecting moons around exoplanets is a major goal of current and future observatories. Moons are suspected to influence rocky exoplanet habitability, and gaseous exoplanets in stellar habitable zones could harbor abundant and diverse moons to target in the search for extraterrestrial habitats. Exomoons contribute to exoplanetary signals but are virtually undetectable with current methods.
Aims. We identify and analyze traces of exomoons in the temporal variation of total and polarized fluxes of starlight reflected by an Earthlike exoplanet and its spatially unresolved moon across all phase angles, with both orbits viewed in an edgeon geometry.
Methods. We compute the total and linearly polarized fluxes, and the degree of linear polarization P of starlight that is reflected by the exoplanet with its moon along their orbits, accounting for the temporal variation of the visibility of the planetary and lunar disks, and including the effects of mutual transits and mutual eclipses. Our computations pertain to a wavelength of 450 nm.
Results. Total flux F shows regular dips due to planetary and lunar transits and eclipses. Polarization P shows regular peaks due to planetary transits and lunar eclipses, and P can increase and/or slightly decrease during lunar transits and planetary eclipses. Changes in F and P will depend on the radii of the planet and moon, on their reflective properties, and their orbits, and are about one magnitude smaller than the smooth background signals. The typical duration of a transit or an eclipse is a few hours.
Conclusions. Traces of an exomoon due to planetary and lunar transits and eclipses show up in the F and P of sunlight reflected by planet–moon systems and could be searched for in exoplanet flux and/or polarization phase functions.
Key words: methods: numerical / polarization / radiative transfer / planetary systems / planets and satellites: detection
© 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 RemoteSensing Infrared Exoplanet Largesurvey; 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 groundbased instruments, such as the Gemini Planet Imager (GPI) instrument (see Macintosh et al. 2014) on the Gemini North telescope, the SpectroPolarimetric Highcontrast 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 (EELT), which is under construction by ESO. These highcontrast 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 nearinfrared 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 Earthlike planets (Stam 2008; Karalidi et al. 2012; Rossi & Stam 2017) and giant Jupiterlike 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.
Extrasolar 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, lifesupporting 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 Ganymedesized 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. Marssized and possibly even Ganymedesized 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 Earthlike 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 Earthlike planet (with a Lambertian reflecting surface and a gaseous atmosphere) with a Moonlike satellite (with a Lambertian reflecting surface) in an edgeon 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^{−2}. 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_{0} = F_{0}1, with 1 the unit column vector and πF_{0} 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_{0}, 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 Diskintegrated 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 diskintegrated 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 ≤ ψ < π).^{1}
To compute the diskintegrated 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 diskintegrated 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_{0} 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 Threedimensional view (left panel) and projection onto the x_{0} –y_{0} plane (right panel) of the discretized planet or moon. The z_{0}axis points toward the observer. The orientation of the x_{0} and y_{0} 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–Moonlike 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 midlatitude 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^{2}) and N the molecular column number density (in m^{−2}) 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 wavelengthdependent 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^{−2}). 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^{2} (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 diskintegrated 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 socalled 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 planeparallel 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 threebody 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 z_{0} 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, threebody problem.
A precise computation of the orbital positions in the generic threebody problem requires the numerical propagation of a given set of initial conditions. Instead, we use the “nested twobody” 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 twobody 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 threebody problem simplification (see, e.g., Wakker 2015), it can handle eccentric orbits.
As demonstrated by Kipping (2010), the nested twobody approximation is excellent for the generic threebody 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 twobody 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 twobody 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, righthanded coordinate system S_{1} = {x_{1}, y_{1}, z_{1}} is the reference frame for the observation of the planet–moon–star system, with the star at the origin, and the z_{1}axis pointing toward the observer. Plane p_{1}, through x_{1} and y_{1}, is the plane on the observer’s sky, onto which the pixels (Fig. 1) are projected. Axes x_{1} and y_{1} have an arbitrary (but fixed) orientation. Coordinate system S_{2} = {x_{2}, y_{2}, z_{2}} is the reference frame for the orbit of the barycenter, which lies in plane p_{2}, through x_{2} and y_{2}. The lunar orbit itself lies in the x_{3} –y_{3} plane of coordinate system S_{3} = {x_{3}, y_{3}, z_{3}}, 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_{2}, and the position of the moon around the barycenter, r_{mb}, in coordinate system S_{3}. The absolute position of the moon in coordinate system S_{2} is then obtained through (14)
As formulated by Murray & Correia (2010), the position of the barycenter and the moon in S_{2} 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_{1}. 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_{1} 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 diskintegrated 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 planetmoon system barycenter (panel b). Plane p_{1} (panel a) is the plane of the sky as seen by the observer on the positive z_{1} axis. Plane p_{2} (panels a and b) is the barycenter’s orbital plane, and plane p_{3} (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 r_{bs} and r_{mb} 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 edgeon geometries, in which the inclination angle of the barycenter’s orbit is 90°, because exoplanets in (near) edgeon 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 followup 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 semimajor axes match those of the Earth–Moon system (Williams 2017). We neglect the Earth–barycenter distance, so that a_{bs} = 1 AU. The semimajor axis of the lunar orbit, a_{mb}, is computed from the Moon–Earth semimajor axis, a_{mp} = 2.5696 × 10^{−3} AU (Williams 2017), as (15)
with m_{p} and m_{m} the masses of the Earth and Moon, respectively. Because of the edgeon 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} = 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.
Orbital elements of the barycenter of our planet–moon system and of the orbit of the moon.
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 tradeoff 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 edgeon planet–moon system at time t = 0 as seen from the positive yaxis (left panel) and from the observer’s position at the positive zaxis (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 diskresolved signals give insight into the reflection processes. For both bodies, total flux F is maximum at the substellar/subobserver 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 substellar 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 substellar 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 diskintegrated F, Q, U, and P as functions of the planetary phase angle α_{p}. For our edgeon 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 diskintegrated 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 planetmoon 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 diskintegrated 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 edgeon system. The positive z_{0} 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, edgeon, 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 diskintegrated P is close to zero due to symmetry. During ingress and egress, the moon breaks the symmetry, and (slightly) increases the diskintegrated 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 lunartoplanetary 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 lunartoplanetary 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 Moonsized, Lambertian (i.e., isotropically and depolarizing) reflecting moon orbiting an Earthsized exoplanet with an Earthlike, gaseous atmosphere on top of a Lambertian reflecting surface (the surface pressure is 1 bar), in an edgeon 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 Earthlike 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 cloudfree 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, atmospheresurface 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 edgeon 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 faceon 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 zaxis, 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 subobserver 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 timeindependent. It is computed according to (A.3)
where subscript “x” refers to either “p” or “m”, u_{z} is the unit vector along the zaxis 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 substellar location) and 90° (at the terminator). The position of the star changes in time, so that the timedependent 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, u_{i}, the observer’s direction unit vector, u_{ob}, and the star direction unit vector, u_{s}. The sides of the spherical triangle are: the observer–zenith angle θ_{i}, the star–zenith angle θ_{0i}, and the pixelbased 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^{6}. 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 diskintegrated signal of the body. The pixel grid across the planet is defined with respect to the planetary scattering plane, and β_{i} is thus timeindependent 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 zaxis points toward the observer).
For the lunar pixels, the alignment between the lunar scattering plane and the lunar grid, and hence angle β_{i}, is timedependent 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, edgeon 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 xaxis and yaxis in coordinate system S_{1}.
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 noneclipsed 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 noneclipsed, 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 noneclipsed 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^{3}. 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_{4} 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 pixelbypixel 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)
The diskof body B will only be partially inside the umbral shadow cone of body A when (B.9)
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 R_{S}, R_{A}, and R_{B}, respectively. Points O_{1}, O_{2}, O_{3}, O_{4}, and O_{5} denote auxiliary points: the umbral and antumbral cones have apex O_{1} and aperture2Ψ and the penumbral cone has apex O_{2} 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 socalled 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)
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 bodytostar angular distance as seen from the shadowed body, δ. Then, as follows from Fig. B.3, c_{i} can be computed using (B.21)
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)
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 A_{1} and A_{2}. 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 l_{S} and l_{A}, 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 tradeoff between the errors and the computational time, with time steps of 24 h.
For the tradeoff, 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 diskintegration 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 N_{eq} 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 N_{eq} 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
 Agnor, C. B., & Hamilton, D. P. 2006, Nature, 441, 192 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Bailey, J. 2007, Astrobiology, 7, 320 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Barnes, J. W., & Fortney, J. J. 2003, ApJ, 588, 545 [NASA ADS] [CrossRef] [Google Scholar]
 Benn, C. R. 2001, The Moon and the Origin of Life, eds. C. Barbieri & F. Rampazzi (Dordrecht, the Netherlands: Springer), 61 [Google Scholar]
 Beuzit, J.L., Feldt, M., Dohlen, K., et al. 2006, The Messenger, 125, 29 [Google Scholar]
 Bott, K., Bailey, J., KedzioraChudczer, L., et al. 2016, MNRAS, 459, L109 [NASA ADS] [CrossRef] [Google Scholar]
 Buenzli, E., & Schmid, H. M. 2009, A&A, 504, 259 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cabrera, J., & Schneider, J. 2007, A&A, 464, 1133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Campbell, B., Walker, G. A., & Yang, S. 1988, ApJ, 331, 902 [NASA ADS] [CrossRef] [Google Scholar]
 Canup, R. M., & Ward, W. R. 2006, Nature, 441, 834 [NASA ADS] [CrossRef] [Google Scholar]
 Carciofi, A. C., & Magalhães, A. M. 2005, ApJ, 635, 570 [NASA ADS] [CrossRef] [Google Scholar]
 de Haan, J. F., Bosma, P. B., & Hovenier, J. W. 1987, A&A, 183, 371 [Google Scholar]
 de Kok, R. J., & Stam, D. M. 2012, Icarus, 221, 517 [NASA ADS] [CrossRef] [Google Scholar]
 De Kok, R., Stam, D., & Karalidi, T. 2011, ApJ, 741, 59 [NASA ADS] [CrossRef] [Google Scholar]
 De Pater, I., & Lissauer, J. J. 2015, Planetary Sciences (Cambridge, UK: Cambridge University Press), 22 [CrossRef] [Google Scholar]
 Domingos, R. C., Winter, O. C., & Yokoyama, T. 2006, MNRAS, 373, 1227 [NASA ADS] [CrossRef] [Google Scholar]
 Forgan, D., & Yotov, V. 2014, MNRAS, 441, 3513 [NASA ADS] [CrossRef] [Google Scholar]
 Ginski, C., Benisty, M., van Holstein, R. G., et al. 2018, A&A, 616, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gratton, R., Kasper, M., Vnaud, C., Bonavita, M., & Schmid, H. M. 2011, IAU Symp., 276, 343 [NASA ADS] [Google Scholar]
 Han, E., Wang, S. X., Wright, J. T., et al. 2014, PASP, 126, 827 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, J. E., & Travis, L. D. 1974, Space Sci. Rev., 16, 527 [NASA ADS] [CrossRef] [Google Scholar]
 Heller, R. 2012, A&A, 545, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heller, R. 2014, ApJ, 787, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Heller, R. 2017, in Handbook of Exoplanets, eds. H. J. Deeg & J. A. Belmonte (New York, NY: Springer), 35 [Google Scholar]
 Heller, R., & Barnes, R. 2013, Astrobiology, 13, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Heller, R., & Barnes, R. 2015, Int. J. Astrobiol., 14, 335 [CrossRef] [Google Scholar]
 Heller, R., Williams, D., Kipping, D., et al. 2014, Astrobiology, 14, 798 [NASA ADS] [CrossRef] [Google Scholar]
 Heller, R., Hippke, M., & Jackson, B. 2016, ApJ, 820, 88 [NASA ADS] [CrossRef] [Google Scholar]
 Hippke, M. 2015, ApJ, 806, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Hough, J. H., & Lucas, P. W. 2003, in Earths: DARWIN/TPF and the Search for Extrasolar Terrestrial Planets, eds. M. Fridlund, T. Henning, & H. Lacoste, ESA SP, 539, 11 [Google Scholar]
 Hough, J. H., Lucas, P. W., Bailey, J. A., & Tamura, M. 2003, A HighSensitivity Polarimeter for the Direct Detection and Characterization of Extrasolar Planets, Proc. SPIE, 4843, 517 [Google Scholar]
 Hovenier, J. W., & van der Mee, C. V. M. 1983, A&A, 128, 1 [NASA ADS] [Google Scholar]
 Hovenier, J. W., van der Mee, C. V., & Domke, H. 2004, Transfer of Polarized Light in Planetary Atmospheres: Basic Concepts and Practical Methods, (Berlin: Springer Science + Business Media), 318 [CrossRef] [Google Scholar]
 Karalidi, Stam, D. M., & Hovenier, J. W. 2012, A&A, 548, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kawata, Y. 1978, Icarus, 33, 217 [NASA ADS] [CrossRef] [Google Scholar]
 Keller, C. U., Schmid, H. M., Venema, L. B., et al. 2010, Proc. SPIE, 7735, 77356G [Google Scholar]
 Kemp, J. C., Henson, G. D., Steiner, C. T., & Powell, E. R. 1987, Nature, 326, 270 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D. M. 2009, MNRAS, 392, 181 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D. M. 2010, MNRAS, 409, L119 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D. M. 2011, The Transits of Extrasolar Planets with Moons, 1st edn., Springer Theses (Berlin: Springer), 127 [Google Scholar]
 Kipping, D. M., Fossey, S. J., Campanella, G., Schneider, J., & Tinetti, G. 2010, in Pathways Towards Habitable Planets, eds. V. Coudé du Foresto, D. M. Gelino, & I. Ribas (San Francisco: ASP), 139 [Google Scholar]
 Kipping, D. M., Schmitt, A. R., Huang, X., et al. 2015, ApJ, 813, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Kostogryz, N., Yakobchuk, T., Morozhenko, O., & VidMachenko, A. 2011, MNRAS, 415, 695 [NASA ADS] [CrossRef] [Google Scholar]
 Kostogryz, N., Yakobchuk, T., & Berdyugina, S. 2015, ApJ, 806, 97 [NASA ADS] [CrossRef] [Google Scholar]
 Lehmer, O. R., Catling, D. C., & Zahnle, K. J. 2017, ApJ, 839, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Link, F. 1969, Eclipse Phenomena in Astronomy (Berlin: Springer) [CrossRef] [Google Scholar]
 Macintosh, B., Gemini Planet Imager Instrument Team, Planet Imager Exoplanet Survey, G., & Observatory, G. 2014, in AAS Meeting Abstracts, 223, 229.02 [Google Scholar]
 Macintosh, B., Graham, J. R., Barman, T., et al. 2015, Science, 350, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Marley, M. S., & Sengupta, S. 2011, MNRAS, 417, 2874 [NASA ADS] [CrossRef] [Google Scholar]
 McClatchey, R., Fenn, R., Selby, J., Volz, F., & Garing, J. 1972, Optical Properties of the Atmosphere (US Air Force Cambridge Research Labs), AFCRL72.0497 [Google Scholar]
 Morbidelli, A., Tsiganis, K., Batygin, K., Crida, A., & Gomes, R. 2012, Icarus, 219, 737 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, C. D., & Correia, A. C. M. 2010, Keplerian Orbits and Dynamics of Exoplanets, ed. S. Seager, 15 [Google Scholar]
 Reynolds, R. T., McKay, C. P., & Kasting, J. F. 1987, Adv. Space Res., 7, 125 [Google Scholar]
 Rosenblatt, P., Charnoz, S., Dunseath, K. M., et al. 2016, Nat. Geosci., 9, 581 [NASA ADS] [CrossRef] [Google Scholar]
 Rossi, L., & Stam, D. 2017, A&A, 607, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rossi, L., & Stam, D. M. 2018, A&A, 616, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rossi, L., BerzosaMolina, J., & Stam, D. M. 2018, A&A, 616, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rufu, R., Aharonson, O., & Perets, H. B. 2017, Nat. Geosci., 10, 89 [Google Scholar]
 Saar, S., & Seager, S. 2003, ASP Conf. Ser., 294, 529 [Google Scholar]
 Sartoretti, P., & Schneider, J. 1999, A&AS, 134, 553 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Scharf, C. A. 2006, ApJ, 648, 1196 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, J., Lainey, V., & Cabrera, J. 2015, Int. J. Astrobiol., 14, 191 [CrossRef] [Google Scholar]
 Seager, S., Whitney, B. A., & Sasselov, D. D. 2000, ApJ, 540, 504 [NASA ADS] [CrossRef] [Google Scholar]
 Sengupta, S. 2016, ApJ, 152, 98 [NASA ADS] [CrossRef] [Google Scholar]
 Sengupta, S., & Marley, M. S. 2016, ApJ, 824, 76 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, A., Szatmáry, K., & Szabó, G. M. 2007, A&A, 470, 727 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Simon, A., Szabó, G. M., Kiss, L., Fortier, A., & Benz, W. 2015, PASP, 127, 1084 [NASA ADS] [CrossRef] [Google Scholar]
 Snik, F., & Keller, C. U. 2013, Astronomical Polarimetry: Polarized Views of Stars and Planets, eds. T. D. Oswalt & H. E. Bond (Dordrecht: Springer), 175 [Google Scholar]
 Stam, D. M. 2003, in Earths: DARWIN/TPF and the Search for Extrasolar Terrestrial Planets, eds. M. Fridlund, T. Henning, & H. Lacoste, ESA SP, 539, 615 [Google Scholar]
 Stam, D. M. 2008, A&A, 482, 989 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stam, D. M., & Hovenier, J. W. 2005, A&A, 444, 275 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stam, D. M., Hovenier, J. W., & Waters, L. B. F. M. 2004, A&A, 428, 663 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stam, D. M., de Rooij, W. A., Cornet, G., & Hovenier, J. W. 2006, A&A, 452, 669 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Stolker, T., Min, M., Stam, D. M., et al. 2017, A&A, 607, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Szabó, G. M., Szatmáry, K., Divéki, Z., & Simon, A. 2006, A&A, 450, 395 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tinetti, G., Drossart, P., Eccleston, P., et al. 2016, in The science of ARIEL, Proc. SPIE, 9904, 99041X [CrossRef] [Google Scholar]
 Wagner, K., Apai, D., Kasper, M., et al. 2016, Science, 353, 673 [NASA ADS] [CrossRef] [Google Scholar]
 Wakker, K. F. 2015, Fundamentals of Astrodynamics (Delft, NL: TU Delft Library) [Google Scholar]
 Wiktorowicz, S. J., & Laughlin, G. P. 2014, ApJ, 795, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Williams, D. 2017, NASA Planetary Fact Sheet – Metric, https://nssdc.gsfc.nasa.gov/planetary/factsheet/ [Google Scholar]
 Williams, D., & Knacke, R. 2004, Astrobiology, 4, 400 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Wolszczan, A., & Frail, D. A. 1992, Nature, 355, 145 [NASA ADS] [CrossRef] [Google Scholar]
Hovenier & van der Mee (1983) define ψ while rotating in the counterclockwise direction when looking toward the observer, which yields the same angle.
PyMieDAP is freely available under the GNU GPL license at http://gitlab.com/loic.cg.rossi/pymiedap.
All Tables
Orbital elements of the barycenter of our planet–moon system and of the orbit of the moon.
All Figures
Fig. 1 Threedimensional view (left panel) and projection onto the x_{0} –y_{0} plane (right panel) of the discretized planet or moon. The z_{0}axis points toward the observer. The orientation of the x_{0} and y_{0} 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 z_{0} 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 planetmoon system barycenter (panel b). Plane p_{1} (panel a) is the plane of the sky as seen by the observer on the positive z_{1} axis. Plane p_{2} (panels a and b) is the barycenter’s orbital plane, and plane p_{3} (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 r_{bs} and r_{mb} 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 edgeon planet–moon system at time t = 0 as seen from the positive yaxis (left panel) and from the observer’s position at the positive zaxis (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 diskintegrated 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 edgeon system. The positive z_{0} 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, edgeon, 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 lunartoplanetary 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 lunartoplanetary 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, u_{i}, the observer’s direction unit vector, u_{ob}, and the star direction unit vector, u_{s}. The sides of the spherical triangle are: the observer–zenith angle θ_{i}, the star–zenith angle θ_{0i}, and the pixelbased 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 R_{S}, R_{A}, and R_{B}, respectively. Points O_{1}, O_{2}, O_{3}, O_{4}, and O_{5} denote auxiliary points: the umbral and antumbral cones have apex O_{1} and aperture2Ψ and the penumbral cone has apex O_{2} 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 A_{1} and A_{2}. 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 l_{S} and l_{A}, respectively. 

In the text 
Fig. C.1 Analysis for the number of pixels along the equator N_{eq} 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 N_{eq} 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.