Traces of exomoons in computed flux and polarization phase curves of starlight reflected by exoplanets
Javier Berzosa Molina Loic Rossi Daphne M. Stam
 







Faculty of Aerospace Engineering, Delft University of Technology,
           Kluyverweg 1, 2629 HS Delft, The Netherlands

J. Berzosa Molina
j.berzosamolina@gmail.com






   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 will contribute to exoplanetary 
    signals but are virtually undetectable with current methods.
  
   We identify and analyze traces of exomoons in the temporal variation of
    total and polarized fluxes of starlight reflected by 
    an Earth–like exoplanet and its spatially unresolved moon across all 
    phase angles, with both orbits viewed in an edge–on geometry. 
  
   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 
    effects of mutual transits and mutual eclipses.
    Our computations pertain to a wavelength of 450 nm.
  
   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.
  
   Traces of an exomoon due to planetary and lunar transits and eclipses
    show up in F and P of sunlight reflected by planet–moon systems
    and could be searched for in exoplanet flux and/or polarization phase 
    functions.




Flux and polarization signals of exoplanets with moons 
Traces of exomoons in 

      computed flux and polarization phase curves 
 
                               of starlight reflected by exoplanets
    J. Berzosa MolinaNow at GMV AD, Calle de Isaac Newton, 11, 
        28760 Tres Cantos, Madrid, Spain 
        L. Rossi D. M. Stam
    Received 30 April 2018; accepted 23 July, 2018
================================================================================================================================================






§ INTRODUCTION



Since the detection of the first planets beyond 
our Solar System <cit.>, 
the number of discoveries has steadily increased, 
yiel­ding almost 4000 confirmed exoplanets and 2500 unconfirmed, 
candidate exoplanets to this day <cit.>. 
Exoplanet space telescopes, such as ESA's CHEOPS 
(CHaracterising ExOPlanet Satellite) and Plato (PLAnetary Transits and 
Oscillations of stars), and NASA's Transiting Exoplanet Survey Satellite (TESS), 
are dedicated to find exoplanets around bright, nearby stars. The relative 
small distances to these stars and their planets combined with the high sensitivity 
of these missions and the upcoming JWST/NASA (James Webb Space Telescope) 
and ARIEL/ESA <cit.> missions will allow 
the search for lunar companions and planetary rings.

The continuous increase in instrument precision, stability and spatial 
resolution has enabled a new generation of 
ground–based instruments, such as the Gemini Planet Imager (GPI) 
instrument <cit.> on the Gemini 
North telescope, the Spectro–Polarimetric High–contrast Exoplanet 
Research (SPHERE) instrument <cit.> on 
ESO's Very Large Telescope (VLT) and the proposed Exoplanet Ima­ging 
Camera and Spectrograph (EPICS) <cit.> 
on the European Extremely Large Telescope (E-ELT), which is under 
construction by ESO. These high–contrast instruments use 
direct imaging of planetary radiation to not only detect but also
characterize exoplanetary systems through a combination of
spectroscopy and polarimetry techniques. Using GPI and SPHERE, 
respectively, <cit.> and <cit.> announced 
the discoveries of young, and thus hot Jovian planets, whose atmospheric 
pro­perties and orbits were characterized using near–infrared spectroscopy. 

Direct detection of extrasolar bodies presents a major cha­llenge 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 of the planet properties that cannot be obtained through other 
methods, such as lower atmospheric composition and, for rocky planets, 
their surface coverage.

Despite not having been exploited yet, great hope is placed in 
the polarimetry capabilities of current and future telescopes as 
powerful tools for detecting and characterizing exoplanets 
<cit.>. 
Previous works in this 
field involve the modeling of stellar polarization during 
planetary transits <cit.>, 
the modeling of light curves and 
polarization of starlight reflected signals in the visible range 
of Earth-like planets <cit.>, and 
giant Jupiter-like planets 
<cit.>, as well as the mode­ling of 
exoplanetary atmospheres in the infrared <cit.>, 
which demonstrated the usefulness of direct observations on exoplanet 
characterization. More recently, <cit.> reported 
linear polarization observations of the hot Jupiter system HD 
189733, and <cit.> 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 exo­moon detection: <cit.> 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 symme­tric cloudy planets. 
Studying exomoons 
can improve our understanding of in particular: 





  * 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 
<cit.>. 
Triton might have been captured by Neptune
<cit.>, while collisions are thought to have 
altered the relative alignment between Uranus and its moons 
<cit.>. Indeed, studying Solar System 
moons gives essential insights on formation mechanisms and evolution 
<cit.>. Exomoon 
research would allow refining planet formation theories 
in a way not achievable by studying exoplanets alone.



  * Extra-solar system characterization: studying 
exomoons will not only provide information on lunar orbits and 
physical properties, but will also allow constraining planet 
characteristics such as i.e. mass, oblateness, 
and rotation axis <cit.>. A signal of a planet-moon system could 
be interpreted as that of a planet alone, resulting in e.g. an 
overestimation of the planet mass and effective temperature 
<cit.>, and/or an anomalous composition 
from spectroscopy <cit.>. 
Extrasolar system charac­terization would indeed require 
analysis of all of its elements, i.e. planets, 
moons, rings, and exozodiacal dust.  



  * Exoplanet and exomoon habitability:  
a moon may influence its planet's habitability  
<cit.>, and moons of giant exoplanets within the stellar 
habitable zone (HZ) might host habitable environments 
<cit.>. 
<cit.> and <cit.> mention 
the role of a moon's orbit on the presence of liquid, life–supporting 
water. Indeed, tidal heating could maintain surface temperatures 
compatible with life on large moons around cold giant planets 
<cit.>. 
<cit.> show that small moons 
could retain atmospheres over limited time periods, 
while Ganymede–sized moons in a stellar HZ could hold  
atmospheres and surface water indefinitely. Although radia­tion in 
a giant planet's magnetic field and eclipses could threaten 
local conditions for life 
<cit.>,
exomoons are interes­ting 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 Varia­tions 
(TTV) and Transit Duration Variations (TDV) data 
<cit.>,
and Hippke's search using the Orbital Sampling Effect (OSE) 
<cit.>, 
the search for exomoons is in its starting phase. Mars–sized and 
possibly even Ganymede-sized satellites could be 
traceable in archived Kepler data <cit.>. 
Unfortunately, as of yet no exomoons have been confirmed.
    
In this paper, we use numerical simulations to 
show how an exomoon could influence 
the flux and degree of polarization of the starlight that is reflected 
by an Earth-like exoplanet, using the following outline. 
In Sect. <ref>, 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. <ref>, we present computed flux and polarization 
phase functions at 450 nm, for an Earth-like planet 
(with a Lambertian reflecting surface and a gaseous atmosphere) with a 
Moon–like satellite (with a Lambertian reflecting surface)
in an edge–on geometry.
Finally, in Sect. <ref>, we summarize and discuss our findings
and their implications.




§ COMPUTING THE REFLECTED STARLIGHT






 §.§ Stokes vectors and polarization



We describe the flux and polarization of starlight that is reflected 
by a body, with a Stokes vector <cit.>:

    𝐅 = [ [ F; Q; U; V ]],

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 usua­lly much smaller than the linearly polarized fluxes 
<cit.>, and because
ignoring V does not 
lead to significant errors in the computation of F, Q, and U 
<cit.>. The light of the star is assumed to be 
unpola­rized <cit.>, and is 
given by 𝐅_ 0 = F_01, with 1 the 
unit column vector and π F_0 the flux measured perpendi­cular 
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 𝐅 is defined as

    P = √(Q^2+U^2)/F,

and the direction of polarization, χ, with respect to the
reference plane can be computed from

    tan2χ = U/Q,

where χ is chosen such that 0 ≤χ < π, while 
cos2χ and Q have the same sign 
<cit.>.




 §.§ Disk–integrated reflected Stokes vectors



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

    𝐅 = 
           𝐅^ p + 
           R^2_ m/R^2_ p𝐋 (ψ) 𝐅^ m.


Vectors 𝐅^ p and 𝐅^ m are disk–integrated
vectors that include the effects of eclipses and transits. 
They are norma­lized such that the total fluxes reflected by 
the planet and moon at a phase angle α= 0^∘ and without
shadows and/or eclipses on their disks, 
equal the planet's and moon's geometric albedo's, res­pectively 
<cit.>.
Furthermore, R_ p and R_ m are the radii of the (spherical) 
planet and moon, respectively.

Vectors 𝐅 and 𝐅^ p in Eq. <ref>
are defined with 
respect to the planetary scattering plane, while 𝐅^ 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 scatte­ring plane.
Matrix L in Eq. <ref> rotates  
𝐅^ m from the lunar to the planetary scattering plane. 
It is given by <cit.>

    𝐋(ψ) = [ [     1     0     0     0;     0 cos2ψ sin2ψ     0;     0 sin2ψ cos2ψ     0;     0     0     0     1 ]] ,

with ψ the rotation angle measured in the clockwise direction from
the lunar to the planetary scattering plane when looking towards the moon 
(0 ≤ψ < π).[<cit.> define ψ  
while rotating in the anti–clockwise direction when looking 
towards the observer, which yields the same angle.]









To compute the disk–integrated vectors 𝐅^ p and 𝐅^ 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. <ref>).
The number of pixels on the planetary disk is N_ p 
and that on the lunar disk 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
App. <ref> for insight into the effect of the number of pixels on 
the computed signals).
The disk–integrated vectors are obtained by summing up the contributions
of the individual pixels across the disk, fully taking into account shadows and 
eclipses, i.e.

    𝐅^ x = 
        π/N_ x∑_i=1^N_ x b_i  c_i  d_i 𝐋(β_i) 
        𝐅^ x_i,

where 'x' is either 'p' or 'm'. Factor π/N_ x is the 
surface area per pixel. 
Furthermore, 𝐅^ x_i is the reflected Stokes vector 
for the i-th pixel on the planet (x=p) or moon (x=m), the computation 
of which is described in Sect. <ref>.
Matrix 𝐋 is a rotation matrix (see Eq. <ref>) that is
used for rotating the local Stokes vector 𝐅^ x_i
that is defined with respect to the local reference plane, 
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. <ref>.
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 

    d_i = ( r_ ref / r_i s)^2,

where r_ ref is the reference distance at which the standard 
stellar flux is defined and r_i s is the actual distance 
between pixel i and the star.




 §.§ The locally reflected starlight



The Stokes vector of the starlight that is reflected by
pixel i on the planet or moon is computed using <cit.>:

    𝐅_i^ x(θ_i,θ_0i,ϕ_i-ϕ_0i) = 
              cosθ_0i 𝐑_1i^ x(θ_i,θ_0i,ϕ_i-ϕ_0i) 
              F_0,

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, i.e. 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 
<cit.>.  
Furthermore, 𝐑^ x_1i 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 unpola­rized (cf. Sect. <ref>). 
For a given pixel, the illumination and viewing angles, and thus 
𝐑^ x_1i, depend on the position of the planet or moon
with res­pect to the star and to each other. 
Local reflection matrix 𝐑^ x_i also depends 
on the local composition and structure of the atmosphere and/or 
surface of the reflecting body.
We compute reflected starlight for an Earth–Moon–like planetary 
system, keeping the reflection models for the Earth and the moon simple 
to avoid introducing too many details that increase computational times while 
not adding insight into the observable signals.

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

    b^ m_ sca(λ) = σ^ m_ sca(λ) 
                                       N,

with σ^ m_ sca 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

    σ^ m_ sca(λ) = 
            24 π^3/N_ L^2(n(λ)^2 - 1)^2/(n(λ)^2 + 2)^26 + 3 δ(λ)/6 - 7 δ(λ)1/λ^4,

with N_ L, Loschmidt's number at standard pressure and temperature,
n, the wavelength dependent refractive index of dry air under
standard pressure and temperature, and δ, the
depolarization factor of the atmospheric gas 
<cit.>.
To calculate the molecular column number density N, we assume 
hydrostatic equilibrium in each atmospheric layer, thus

    N = δ p/m g,

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 b^ m_ sca 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
<cit.>. The single scattering
albedo of the gaseous molecules can thus be assumed to equal 1.0.
And, 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
<cit.>.
Our model moon has no atmosphere above its flat, Lambertian  
(i.e. isotropic and depolarizing)
reflecting surface with a_ surf= 0.1 <cit.>.

The computation of the local illumination and viewing geometries 
θ_i, θ_0i, and ϕ_i-ϕ_0i is described in 
Appendix <ref>.
Given these angles and the planet's atmosphere–surface model, 
we use PyMieDAP[PyMieDAP is freely available under the GNU GPL 
license at http://gitlab.com/loic.cg.rossi/pymiedap] 
<cit.>, an efficient radiative transfer code based 
on the adding–doubling algorithm described by <cit.>. 
PyMieDAP fully includes polarization for all orders of scattering,
and assumes a locally plane–parallel atmosphere–surface model
to compute 𝐑_1i^ p for every pixel on the planet.
The computed locally reflected Stokes vector, 𝐅_i^ p, 
is defined with respect to the local meridian plane, 
i.e. the plane through the local zenith and the local direction towards the 
observer. For each illuminated pixel on the moon, 
𝐑_1i^ m = a_ surf1.
A detailed description of PyMieDAP including benchmark results
can be found in <cit.>. 

Results of our radiative transfer code have been compared against results 
presented in e.g. 
<cit.> (who all used the same adding–doubling code, 
but an entirely different disk–integration algorithm), and
<cit.> (who used their own version of an adding–doubling code
and an independent disk–integration method).
<cit.> and <cit.>, 
each compared their own,
independently implemented Monte Carlo radiative transfer codes
successfully against results from the code used by 
<cit.> and <cit.>.












 §.§ Computing transits and eclipses



As described in Eq. <ref>, the contribution of the light reflected 
by a pixel i on the planet or the moon to the disk–integrated 
Stokes vector 𝐅, 
depends on the factors b_i and c_i, that account for the pixel's 
visibility and dimming, respectively. The values of these factors depend on 
so–called mutual events, specifically, transits, in
which one body is (partially) blocking the light that is reflected
towards 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. <ref>):


  * 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.

  * 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.

  * 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. 

  * 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, i.e.the epochs in which these bodies move in front or behind the star. 
Numerical simulations of transiting planets with moons 
have been ­published by <cit.>. 
Modeling the transmission and scattering of starlight in the 
planetary atmosphere
during those epochs (which is not included in the work by <cit.>),
requires a fully spherical atmosphere model instead of a locally
plane-parallel one <cit.> 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 reflected
starlight, we assume that the bodies are at 'infinite' distance of 
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 fo­llow the mathematical description of 
eclipses in the Moon–Earth system­ as developed by <cit.>,
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 en penumbra shadow regions (for an example of the umbra and
penumbra, see Fig. <ref>). 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. <ref>). 
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. <ref> 
is provided in Appendix <ref>.
This computation requires the positions of the planet and the moon
with respect to the star across time, and thus the dynamics of the 
three–body system. The basics of this dynamics is outlined in the 
next section.




 §.§ Computing the orbits of the planet & 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 evalua­ting the disk–integration according to Eq. <ref>.
Both the motions of the planet and its moon around the star
depend on their mutual gravitational interactions. 
Assuming each body attracts as a point mass and neglecting the gravity of other
planets and/or moons in the system, our star–planet–moon 
system is a classical, generic three–body problem. 

A precise computation of the orbital positions in the generic 
three–body problem 
requires the numerical propagation of a given set of initial conditions. 
Instead, we use the 'nested two–body' appro­ximation
described by <cit.>, which assumes 
that the orbits of the planet and moon around the
planet–moon system barycenter, and the orbit of this 
barycenter around the star can all be described by Keplerian orbits.
The advantages of the nested two–body approximation are:
1. the solution can be described analytically;
2. the computational time is significantly shorter than with 
numerical integrations;
3. it provides better insight in the computed orbits as the
elements of all orbits can be specified;
4. unlike the circular restricted three–body problem 
simplification <cit.>, 
it can handle eccentric orbits.

As demonstrated by <cit.>, the nested two–body
approximation is excellent for the generic three–body problem 
provided ≤ 0.531, where  is the 
moon–planet separation in units of the planet's Hill's sphere radius
<cit.>. 
As follows from <cit.> and <cit.>, 
stable, prograde orbiting moons should ­fulfill 
≤ 0.4895, while retrograde orbiting moons 
can be stable up to ≈ 0.9309. 
The nested two–body approximation can thus be applied to all 
prograde orbiting moons, while retrograde orbiting moons are only 
partially ­covered, depending on . 
We will limit ourselves to prograde orbiting moons, as we do not 
expect any influence of the moon's orbital direction on the 
magnitude of reflected flux and polarization features, except 
on their timing.

Figure <ref> shows the geometry of the planet–moon 
system with the reference frames describing the orbit of the planet–moon 
barycenter around the star, and the orbit of the moon around the 
barycenter. The nested two–body approximation assumes that the motions of the
planet–moon barycenter around the star and that of the moon around the 
barycenter are independent.
Orthonormal, right–handed 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 towards the observer. 
Plane p_1, through x_1 and y_1, is the plane on the observer's
sky, onto which the pixels (Fig. <ref>) 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} that is centered at the 
barycenter's position. 
The various orbital elements in these coordinate systems
are indicated as follows:

    [                                              a                                semi–major axis;                                              e                            eccentricity [0, 1];                                              i                     inclination angle [0, 180];                                              ω                argument of periastron [0, 360];                                              Ω right ascension of the ascending node [0, 360];                                              ν                          true anomaly [0, 360] ]

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,
i.e.
    E(t) - e sinE(t) = M(t),
	
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:

    tanν(t)2 = √(1+e1-e)tanE(t)2.
	

We solve for ν(t) for each orbit in the appropriate reference system
by applying the Newton–Raphson method <cit.>
to Eqs. <ref> and <ref>.

We compute the position of
the barycenter, 𝐫_ bs, in coordinate system S_2,
and the position of the moon around the barycenter, 
𝐫_ mb, in coordinate system S_3.
The absolute position of the moon in coordinate system S_2
is then obtained through:

    𝐫_ ms(t) =  𝐫_ mb(t) + 
        𝐫_ bs(t) .


As formulated by <cit.>, 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 
coor­dinate system S_1. 
For further details on these transformation matrices, see  
<cit.>. 

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. <ref>), and the angles
β_i that are used to rotate locally computed Stokes vectors to 
the planetary and lunar scattering planes (Eq. <ref>),
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. <ref>),
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 
App. <ref>.
For each t, we also compute angle ψ to rotate F^ m,
the disk-integrated Stokes vector for the moon, to the planetary
scattering plane (Eq. <ref>).




 §.§ Our baseline planet–moon system



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








Table <ref> lists the orbital elements of our baseline 
planet–moon system.
Both the barycenter's and the lunar orbit are assumed­ to be circular (e=0.0), 
and their semi–major axes match those of the Earth–Moon system <cit.>. 
We neglect the Earth–barycenter distance, so that a_ bs= 1 AU. 
The semi–major axis of the lunar orbit, a_ mb, is computed from the 
Moon–Earth semi–major axis, a_ mp = 2.5696 · 10^-3 AU <cit.>, 
as follows: 

    a_ mb = a_ mpm_ p/m_ p + m_ m≈ 2.54 · 10^-3 AU ,

with m_ p and m_ m the masses of the Earth and Moon, 
respectively.
Because of the edge–on geometry, the right ascensions of the 
ascending nodes of the orbits of the barycenter and the moon are 
set to zero.
Because both orbits are assumed to be circular, their perihelions 
are undefined. The barycenter's argument of perihelion is chosen 
precisely behind the star at time t = t_0 = 0, 
i.e. ω_ b = 270. 
For the moon, ω_ m is set to zero.
The observational and orbital geometry at t=0 is sketched in 
Fig. <ref>.
We use R_ p= 6371.0 km and R_ m= 1737.4 km for
the baseline radii of the planet and the moon, respectively.














§ NUMERICAL RESULTS



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




 §.§ Reflection by the spatially resolved planet & 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 <ref> 
shows the elements of the locally reflected Stokes vectors 
𝐅^ p and 𝐅^ 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. <ref>a), 
both the planet and the moon would be behind the star
and thus invisible, but their disk–resolved signals give 
insight in the reflection processes. For both bodies, total flux F is 
maximum at the sub-stellar/sub-observer region
and decreases towards 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 χ undefined for the moon. The linearly polarized 
fluxes Q and U of the planet are due to Rayleigh scattering in the 
planet's atmosphere. At the sub-stellar region, both Q and U are zero
because of symmetry. The general increase of Q and U towards 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. 
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 <cit.>.

At α= 50^∘ (Fig. <ref>b), 
the total flux F of the moon is maximum at 
the sub-stellar region and decreases towards the terminator, 
due to the isotropic surface reflection and the absence of an atmosphere. 
The planet also shows a decrease of F towards the terminator,
but the location of the flux maximum is more diffuse and more towards
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 anti–symmetric
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. <ref>b 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.

















 §.§ Reflection by the spatially unresolved planet & moon



Figure <ref> shows the disk–integrated F, Q, U, and P
as functions of the planetary phase angle α_ p.  
For our edge–on system, the curves cover only half of the barycenter's 
orbital period.
For comparison, we have also included curves for the planet and the
moon as isolated bodies, thus without any mutual events.
The total flux of the planet–moon system is lower than the sum
of the total fluxes of the isolated planet and moon because the latter
have not been scaled to the actual radii of the planet and the moon.
Indeed the moon's flux at α_ p=0^∘, equals the moon's geometric
albedo, i.e. 0.067, which matches the theoretical geometric albedo 
of a Lambertian reflecting body with surface albedo of 0.1
<cit.>. The geometric albedo of the unresolved 
planet–moon system (with both bodies at α_ p=0^∘ and next to
each other) is about 0.33. 
Note that in Fig. <ref>, 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. <ref>, and from situation 1
in Fig. <ref>).

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 to that of the planet is due to 
the scattering of light in the latter's atmosphere. 
As can be seen in Fig. <ref>, 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, atmospheric and surface properties. 

Because the lunar surface is completely depolarizing, 
the moon's polarized fluxes Q and U are zero at each 
α_ p. The disk–integrated U of the light 
reflected by the planet is zero due to symmetry (see Fig. <ref>).
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 would contain 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 
<cit.>.

While the smooth curves for the spatially unresolved planet–moon system
shown in Fig. <ref>, 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 <ref> illustrates the various events. 
Both the planet and the moon are initially (α_ p=0^∘) 
behind the star (position 1 in Fig. <ref>).
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 (cf. Fig. <ref>) (and the total flux F
thus slightly larger) 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. <ref>.
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 due to 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. <ref>.

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

As can be seen in Fig. <ref>, 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. <ref>). 
The time gap between two consecutive transits and eclipses, however, 
increases with increasing α_ p because of the movement 
of the barycenter along its orbit.




 §.§ Analysis of the mutual events



In this section, we analyze individual mutual events, i.e. 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 follows

    Δ F     =     F_ event - F_ continuum, 
    Δ P     =      P_ event - P_ continuum.

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












Lunar transits



Figures <ref>a sand b provide detailed views of Δ F and 
Δ P during the six lunar transits (moon in front of the planet) 
shown in Fig. <ref>,
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 expected with 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. <ref>). 
Because egress takes place over the planetary limb, all curves in 
Figs. <ref>a 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. <ref>), 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 night–side starts to appear in the curve. 
Because of the moon's prograde orbit, the lunar day–side ingresses before 
the night–side. In the curves for α_ p=13.4^∘,
and 40.3^∘, the steeper decrease of F due to the ingress of the
lunar night–side 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. <ref>a 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.

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 non–linearly 
with increasing lunar radii, we have aligned the mutual events in 
Fig. <ref> in time to facilitate a comparison. 
Mutual transits show up every half lunar sidereal period. 
Because T ∝√(1/(M_ m+M_ p)), relative timing variations 
of 1% to 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% to 14% are obtained for the same 
range of radii ratios.


















Figure <ref>b 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. <ref>.
The curves exhibit a strong variation in shapes, and with increasing 
α_ p, get increasingly asymmetric. 
The largest Δ P is found around α_ p=90^∘, 
where P of our model planet is highest
(see Fig. <ref>).
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. <ref>).
At α_ p=13.4^∘, Q is maximum near the planet's limb and 
close to zero at the center of the body. The disk–integrated P is close to 
zero due to symmetry. 
During ingress and egress, the moon breaks the symmetry, 
and (slightly) increases the disk–integrated value of P (see Fig. <ref>b). 
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 <ref>b 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 seen to be very sensitive to the lunar 
radius. Indeed, with increasing r, Δ P changes from positive 
(P during the event is higher than in the continuum) to negative 
(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.




Planetary transits 



Figures <ref>c and d show the planetary transits 
(planet in front of the moon) from Fig. <ref>
in detail, both for 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, similarly 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 both on the Δ F and on the planet's
polarized flux Q, and thus on α_ p for a given planet-moon
model.  
In Fig. <ref>d, 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. 




Planetary eclipses 



Figures <ref>a and b show the curves of Δ F 
and Δ P during the planetary eclipse events shown before in 
Fig. <ref>. 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 <ref> illustrates the geometries for the 
planetary eclipse, with φ_ ms the angle between the star 
and the moon measured positive in the counter clock–wise direction
from the center of the planet. Angle φ_ ms is used 
as 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. <ref>a 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. <ref>). 
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. <ref> and <ref>. However, because the 
spatial extent of the penumbral and umbral shadows across the disk are 
smallest halfway 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.
Note that differently than for lunar 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 Figs. <ref>c and d,
which show 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. <ref>b.
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 <ref>d shows P at 
α_ p=72.5^∘ for various values of r (the lunar radius
expressed as 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 decreases, too, apparently because the
polarized flux Q decreases more than the total flux F.













Lunar eclipses 



As can be seen in Figs. <ref>c 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
(Fig. <ref>c 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.

Like with the planetary eclipses discussed above, the moon's ingress 
into the planetary shadow occurs at the same value of 
angle φ_ ms (cf. Fig. <ref>) 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. <ref>, 
is similar to that during planetary transits (Fig. <ref>), 
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. <ref>.

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.




§ CONCLUSIONS



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








During events that darken the planet, i.e. the 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% - 1.8 % 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, i.e. 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 
non–polarizing 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, i.e. in the blue, where the 
scattering by the gas in the Earth–like planetary atmosphere 
strongly contributes to the planet's polarization signals. 
In particular at intermediate phase angles, the polarization
signal of a gaseous atmosphere is strong, and the change in polarization 
during the lunar darkening events can reach a few percent. Indeed, 
Δ P ≈ 1.25% - 2.66 % for 
54 < α_ p < 108 during planetary transits.
Note that at these phase angles, the angular distance between the planet–moon
system and the parent star is relatively large, so these angles are favorable for
the detection of reflected starlight.
For a planet with clouds in its atmosphere, the 
continuum flux phase function will have a similar shape as that for
our cloud–free planet, except the total amount of reflected flux 
will be larger (of course not at wavelengths where atmospheric gases absorb
the light). The polarization curve of a cloudy planet will show
angular features due to the scattering of light by the cloud particles,
such as the rainbow for liquid water droplets 
<cit.>.

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 ∼ 4 hours both in flux and polarization. 
A lunar transit at an intermediate phase angle of 90^∘, 
takes about 2 hours 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 hours, where the
polarization change in the planetary eclipses is only apparent
during part of the time of the flux change.

The results presented in this paper correspond to half of the planetary 
orbit around the star. The results for the other half of the orbit will be 
similar, except that the curves will be mirrored with respect to the central
event time, because transit and eclipse ingresses and egresses will happen
over the other side of the darkened body. 

Our results show that measuring the temporal variation in F and/or P
during transits and eclipses could provide extra information on the 
properties of a planet and/or moon and their orbits. Extracting such 
information, however, requires not only detecting such events but also
measuring the shape of the variations in F and/or P. 
For the interpretation of such measurements, numerical simulations 
to map in more detail the influence of the physical characteristics of 
the moon and the planet (radius, albedo, atmosphere-surface
properties) and their orbital characteristics (inclination angles, ellipticity)
on the temporal variation in F and P are required. 
Such simulations will be targeted in future research.




















1000


[Agnor & Hamilton(2006)]agnor2006neptune
Agnor, C. B. & Hamilton, D. P. 2006, Nature, 441, 192

[Bailey(2007)]2007AsBio...7..320B
Bailey, J. 2007, Astrobiology, 7, 320

[Barnes & Fortney(2003)]barnes2003measuring
Barnes, J. W. & Fortney, J. J. 2003, The Astrophysical Journal, 588, 545

[Benn(2001)]Benn2001
Benn, C. R. 2001, The Moon and the Origin of Life (Dordrecht: Springer
  Netherlands), 61–66

[Beuzit et al.(2006)Beuzit, Feldt, Dohlen, Mouillet,
  Puget, Antichi, Baruffolo, Baudoz, Berton, Boccaletti,
  Carbillet, Charton, Claudi, Downing, Feautrier, Fedrigo, Fusco,
  Gratton, Hubin, Kasper, Langlois, Moutou, Mugnier, Pragt,
  Rabou, Saisse, Schmid, Stadler, Turrato, Udry, Waters, &
  Wildi]beuzit06
Beuzit, J.-L., Feldt, M., Dohlen, K., et al. 2006, The Messenger, 125

[Bott et al.(2016)Bott, Bailey, Kedziora-Chudczer, Cotton, Lucas,
  Marshall, & Hough]bott2016polarization
Bott, K., Bailey, J., Kedziora-Chudczer, L., et al. 2016, Monthly Notices of
  the Royal Astronomical Society: Letters, 459, L109

[Buenzli & Schmid(2009)]buenzli2009grid
Buenzli, E. & Schmid, H. M. 2009, Astronomy & Astrophysics, 504, 259

[Cabrera & Schneider(2007)]cabrera2007detecting
Cabrera, J. & Schneider, J. 2007, Astronomy & Astrophysics, 464, 1133

[Campbell et al.(1988)Campbell, Walker, &
  Yang]campbell1988search
Campbell, B., Walker, G. A., & Yang, S. 1988, The Astrophysical Journal, 331,
  902

[Canup & Ward(2006)]canup2006common
Canup, R. M. & Ward, W. R. 2006, Nature, 441, 834

[Carciofi & Magalhães(2005)]2005ApJ...635..570C
Carciofi, A. C. & Magalhães, A. M. 2005, , 635, 570

[de Haan et al.(1987)de Haan, Bosma, & Hovenier]deHaan87
de Haan, J. F., Bosma, P. B., & Hovenier, J. W. 1987, AAP, 183, 371

[De Kok et al.(2011)De Kok, Stam, & Karalidi]de2011
De Kok, R., Stam, D., & Karalidi, T. 2011, The Astrophysical Journal, 741, 59

[de Kok & Stam(2012)]2012Icar..221..517D
de Kok, R. J. & Stam, D. M. 2012, , 221, 517

[De Pater & Lissauer(2015)]de2015planetary
De Pater, I. & Lissauer, J. J. 2015, Planetary sciences (Cambridge University
  Press), 22–55

[Domingos et al.(2006)Domingos, Winter, & Yokoyama]Domingos2006
Domingos, R. C., Winter, O. C., & Yokoyama, T. 2006, Monthly Notices of the
  Royal Astronomical Society, 373, 1227

[Forgan & Yotov(2014)]forgan2014effect
Forgan, D. & Yotov, V. 2014, Monthly Notices of the Royal Astronomical
  Society, 441, 3513

[Ginski et al.(2018)Ginski, Benisty, van Holstein,
  Juhász, Schmidt, Chauvin, de Boer, Wilby, Manara, Delorme,
  Ménard, Pinilla, Birnstiel, Flock, Keller, Kenworthy,
  Milli, Olofsson, Pérez, Snik, & Vogt]2018arXiv180502261G
Ginski, C., Benisty, M., van Holstein, R. G., et al. 2018, ArXiv
  e-prints 1805.02261

[Gratton et al.(2010)Gratton, Kasper, Vérinaud, Bonavita, &
  Schmid]gratton10
Gratton, R., Kasper, M., Vérinaud, C., Bonavita, M., & Schmid, H. M. 2010,
  Proceedings of the International Astronomical Union, 6, 343–348

[Han et al.(2014)Han, Wang, Wright, Feng, Zhao, Fakhouri, Brown, &
  Hancock]han2014exoplanet
Han, E., Wang, S. X., Wright, J. T., et al. 2014, Publications of the
  Astronomical Society of the Pacific, 126, 827

[Hansen & Travis(1974)]Hansen1974
Hansen, J. E. & Travis, L. D. 1974, Space Science Reviews, 16, 527

[Heller(2012)]heller2012exomoon
Heller, R. 2012, Astronomy & Astrophysics, 545, L8

[Heller(2014)]heller2014detecting
Heller, R. 2014, The Astrophysical Journal, 787, 14

[Heller(2017)]heller2017detection
Heller, R. 2017, arXiv preprint arXiv:1701.04706

[Heller & Barnes(2013)]heller2013exomoon
Heller, R. & Barnes, R. 2013, Astrobiology, 13, 18

[Heller & Barnes(2015)]heller2015runaway
Heller, R. & Barnes, R. 2015, International Journal of Astrobiology, 14, 335

[Heller et al.(2016)Heller, Hippke, & Jackson]heller2016modeling
Heller, R., Hippke, M., & Jackson, B. 2016, The Astrophysical Journal, 820, 88

[Heller et al.(2014)Heller, Williams, Kipping, Limbach, Turner,
  Greenberg, Sasaki, Bolmont, Grasset, Lewis, et al.]heller2014formation
Heller, R., Williams, D., Kipping, D., et al. 2014, Astrobiology, 14, 798

[Hippke(2015)]hippke2015detection
Hippke, M. 2015, The Astrophysical Journal, 806, 51

[Hough & Lucas(2003)]houghlucas03
Hough, J. H. & Lucas, P. W. 2003, in ESA Special Publication, Vol. 539,
  Earths: DARWIN/TPF and the Search for Extrasolar Terrestrial Planets, ed.
  M. Fridlund, T. Henning, & H. Lacoste, 11–17

[Hough et al.(2003)Hough, Lucas, Bailey, & Tamura]houghetal03
Hough, J. H., Lucas, P. W., Bailey, J. A., & Tamura, M. 2003, in , 517–523

[Hovenier et al.(2004)Hovenier, van der Mee, &
  Domke]hovenier2004transfer
Hovenier, J. W., van der Mee, C. V., & Domke, H. 2004, Transfer of polarized
  light in planetary atmospheres: basic concepts and practical methods, Vol.
  318 (Springer Science & Business Media)

[Hovenier & van der Mee(1983)]HoveniervanderMee83
Hovenier, J. W. & van der Mee, C. V. M. 1983, AAP, 128, 1

[Karalidi et al.(2012)Karalidi, Stam, D. M., & Hovenier, J.
  W.]Karalidi12
Karalidi, Stam, D. M., & Hovenier, J. W. 2012, A&A, 548, A90

[Kawata(1978)]KAWATA1978217
Kawata, Y. 1978, Icarus, 33, 217

[Keller et al.(2010)Keller, Schmid, Venema, Hanenburg, Jager,
  Kasper, Martinez, Rigal, Rodenhuis, Roelfsema, Snik, Verinaud, &
  Yaitskova]keller10
Keller, C. U., Schmid, H. M., Venema, L. B., et al. 2010, in ,
  77356G–77356G–13

[Kemp et al.(1987)Kemp, Henson, Steiner, &
  Powell]KempHensonSteinerEtAl1987
Kemp, J. C., Henson, G. D., Steiner, C. T., & Powell, E. R. 1987, Nature, 326,
  270

[Kipping(2009)]kipping2009transit
Kipping, D. M. 2009, Monthly Notices of the Royal Astronomical Society, 392,
  181

[Kipping(2010)]kipping2010weigh
Kipping, D. M. 2010, Monthly Notices of the Royal Astronomical Society:
  Letters, 409, L119

[Kipping(2011)]Kipping2011
Kipping, D. M. 2011, The Transits of Extrasolar Planets with Moons, 1st edn.,
  Springer Theses (Springer Berlin Heidelberg), 127–164

[Kipping et al.(2009)Kipping, Fossey, Campanella,
  Schneider, & Tinetti]Kipping09a
Kipping, D. M., Fossey, S. J., Campanella, G., Schneider, J., &
  Tinetti, G. 2009, ArXiv e-prints 0911.5170

[Kipping et al.(2015)Kipping, Schmitt, Huang, Torres, Nesvornỳ,
  Buchhave, Hartman, & Bakos]kipping2015hunt
Kipping, D. M., Schmitt, A. R., Huang, X., et al. 2015, The Astrophysical
  Journal, 813, 14

[Kostogryz et al.(2015)Kostogryz, Yakobchuk, &
  Berdyugina]kostogryz2015polarization
Kostogryz, N., Yakobchuk, T., & Berdyugina, S. 2015, The Astrophysical
  Journal, 806, 97

[Kostogryz et al.(2011)Kostogryz, Yakobchuk, Morozhenko, &
  Vid’Machenko]kostogryz2011polarimetric
Kostogryz, N., Yakobchuk, T., Morozhenko, O., & Vid’Machenko, A. 2011,
  Monthly Notices of the Royal Astronomical Society, 415, 695

[Lehmer et al.(2017)Lehmer, Catling, &
  Zahnle]lehmer2017longevity
Lehmer, O. R., Catling, D. C., & Zahnle, K. J. 2017, The Astrophysical
  Journal, 839, 32

[Link(1969)]link69
Link, F. 1969, Eclipse Phenomena in Astronomy (Springer Berlin Heidelberg)

[Macintosh et al.(2014)Macintosh, Gemini Planet Imager
  instrument Team, Planet Imager Exoplanet Survey, &
  Observatory]macintosh2014gemini
Macintosh, B., Gemini Planet Imager instrument Team, Planet Imager
  Exoplanet Survey, G., & Observatory, G. 2014, in American Astronomical
  Society Meeting Abstracts, Vol. 223, American Astronomical Society Meeting
  Abstracts #223, 229.02

[Macintosh et al.(2015)Macintosh, Graham, Barman, De Rosa,
  Konopacky, Marley, Marois, Nielsen, Pueyo, Rajan, Rameau, Saumon, Wang,
  Ammons, Arriaga, Artigau, Beckwith, Brewster, Bruzzone, Bulger, Burningham,
  Burrows, Chen, Chiang, Chilcote, Dawson, Dong, Doyon, Draper, Duchêne,
  Esposito, Fabrycky, Fitzgerald, Follette, Fortney, Gerard, Goodsell,
  Greenbaum, Hibon, Hinkley, Cotton, Hung, Ingraham, Johnson-Groh, Kalas,
  Lafreniere, Larkin, Lee, Line, Long, Maire, Marchis, Matthews, Max, Metchev,
  Millar-Blanchaer, Mittal, Morley, Morzinskim, Murray-Clay, Oppenheimer,
  Palmer, Patel, Patience, Perrin, Poyneer, Rafikov, Rantakyrö, Rice, Rojo,
  Rudy, Ruffio, Ruiz, Sadakuni, Saddlemyer, Salama, Savransky, Schneider,
  Sivaramakrishnan, Song, Soummer, Thomas, Vasisht, Wallace, Ward-Duong,
  Wiktorowicz, Wolff, & Zuckerman]Macintosh15
Macintosh, B., Graham, J. R., Barman, T., et al. 2015, Science

[Marley & Sengupta(2011)]MarleySengupta2011
Marley, M. S. & Sengupta, S. 2011, , 417, 2874

[McClatchey et al.(1972)McClatchey, Fenn, Selby, Volz, &
  Garing]1972McClatchey
McClatchey, R., Fenn, R., Selby, J., Volz, F., & Garing, J. 1972,
  Optical Properties of the Atmosphere, AFCRL-72.0497 (U.S. Air Force
  Cambridge Research Labs)

[Morbidelli et al.(2012)Morbidelli, Tsiganis, Batygin, Crida, &
  Gomes]morbidelli2012explaining
Morbidelli, A., Tsiganis, K., Batygin, K., Crida, A., & Gomes, R. 2012,
  Icarus, 219, 737

[Murray & Correia(2010)]Murray10
Murray, C. D. & Correia, A. C. M. 2010, Keplerian Orbits and Dynamics of
  Exoplanets, ed. S. Seager, 15–23

[Reynolds et al.(1987)Reynolds, McKay, &
  Kasting]reynolds1987europa
Reynolds, R. T., McKay, C. P., & Kasting, J. F. 1987, Advances in Space
  Research, 7, 125

[Rosenblatt et al.(2016)Rosenblatt, Charnoz, Dunseath,
  Terao-Dunseath, Trinh, Hyodo, Genda, & Toupin]rosenblatt2016accretion
Rosenblatt, P., Charnoz, S., Dunseath, K. M., et al. 2016, Nature Geoscience

[Rossi et al.(2018)Rossi, Berzosa-Molina, &
  Stam]rossi2018pymiedap
Rossi, L., Berzosa-Molina, J., & Stam, D. M. 2018, ArXiv e-prints
  1804.08357

[Rossi & Stam(2017)]Rossi2017
Rossi, L. & Stam, D. 2017, Astronomy & Astrophysics

[Rossi & Stam(2018)]RossiStam2018
Rossi, L. & Stam, D. M. 2018, Accepted for publication in A&A
  arXiv 1805.08686

[Rufu et al.(2017)Rufu, Aharonson, & Perets]rufu2017multiple
Rufu, R., Aharonson, O., & Perets, H. B. 2017, Nature Geoscience

[Saar & Seager(2003)]Saar03
Saar, S. & Seager, S. 2003, ASP Conf. Ser., 294, 529

[Sartoretti & Schneider(1999)]sartoretti1999detection
Sartoretti, P. & Schneider, J. 1999, Astronomy and Astrophysics Supplement
  Series, 134, 553

[Scharf(2006)]scharf2006potential
Scharf, C. A. 2006, The Astrophysical Journal, 648, 1196

[Schneider et al.(2015)Schneider, Lainey, &
  Cabrera]schneider2015next
Schneider, J., Lainey, V., & Cabrera, J. 2015, International Journal of
  Astrobiology, 14, 191

[Seager et al.(2000)Seager, Whitney, & Sasselov]Seageretal00
Seager, S., Whitney, B. A., & Sasselov, D. D. 2000, The Astrophysical Journal,
  540, 504

[Sengupta(2016)]sengupta2016polarimetric
Sengupta, S. 2016, The Astronomical Journal, 152, 98

[Sengupta & Marley(2016)]sengupta2016detecting
Sengupta, S. & Marley, M. S. 2016, The Astrophysical Journal, 824, 76

[Simon et al.(2015)Simon, Szabó, Kiss, Fortier, &
  Benz]simon2015cheops
Simon, A., Szabó, G. M., Kiss, L., Fortier, A., & Benz, W. 2015,
  Publications of the Astronomical Society of the Pacific, 127, 1084

[Simon et al.(2007)Simon, Szatmáry, &
  Szabó]simon2007determination
Simon, A., Szatmáry, K., & Szabó, G. M. 2007, Astronomy &
  Astrophysics, 470, 727

[Snik & Keller(2013)]Snik2013
Snik, F. & Keller, C. U. 2013, Astronomical Polarimetry: Polarized Views of
  Stars and Planets, ed. T. D. Oswalt & H. E. Bond (Dordrecht: Springer
  Netherlands), 175–221

[Stam et al.(2006)Stam, De Rooij, Cornet, &
  Hovenier]stam2006integrating
Stam, D., De Rooij, W., Cornet, G., & Hovenier, J. 2006, Astronomy &
  Astrophysics, 452, 669

[Stam(2003)]stam03
Stam, D. M. 2003, in ESA Special Publication, Vol. 539, Earths: DARWIN/TPF and
  the Search for Extrasolar Terrestrial Planets, ed. M. Fridlund,
  T. Henning, & H. Lacoste, 615–619

[Stam(2008)]Stam08
Stam, D. M. 2008, A&A, 482, 989

[Stam et al.(2006)Stam, de Rooij, Cornet, &
  Hovenier]Stam2006
Stam, D. M., de Rooij, W. A., Cornet, G., & Hovenier, J. W. 2006,
  , 452, 669

[Stam & Hovenier(2005)]stametal05
Stam, D. M. & Hovenier, J. W. 2005, A&A, 444, 275

[Stam et al.(2004)Stam, Hovenier, & Waters]Stametal04
Stam, D. M., Hovenier, J. W., & Waters, L. B. F. M. 2004, A&A, 428, 663

[Stolker et al.(2017)Stolker, Min, Stam, Mollière,
  Dominik, & Waters]2017Stolker
Stolker, T., Min, M., Stam, D. M., et al. 2017, , 607, A42

[Szabó et al.(2006)Szabó, Szatmáry, Divéki, &
  Simon]szabo2006possibility
Szabó, G. M., Szatmáry, K., Divéki, Z., & Simon, A. 2006,
  Astronomy & Astrophysics, 450, 395

[Tinetti et al.(2016)Tinetti, Drossart, Eccleston,
  Hartogh, Heske, Leconte, Micela, Ollivier, Pilbratt, Puig,
  Turrini, Vandenbussche, Wolkenberg, Pascale, Beaulieu, Güdel,
  Min, Rataj, Ray, Ribas, Barstow, Bowles, Coustenis, Coudé
  du Foresto, Decin, Encrenaz, Forget, Friswell, Griffin, Lagage,
  Malaguti, Moneti, Morales, Pace, Rocchetto, Sarkar, Selsis,
  Taylor, Tennyson, Venot, Waldmann, Wright, Zingales, &
  Zapatero-Osorio]2016SPIE.9904E..1XT
Tinetti, G., Drossart, P., Eccleston, P., et al. 2016, in ,
  Vol. 9904, Space Telescopes and Instrumentation 2016: Optical, Infrared, and
  Millimeter Wave, 99041X

[Wagner et al.(2016)Wagner, Apai, Kasper, Kratter, McClure,
  Robberto, & Beuzit]Wagner673
Wagner, K., Apai, D., Kasper, M., et al. 2016, Science, 353, 673

[Wakker(2015)]wakker2015fundamentals
Wakker, K. F. 2015, Fundamentals of astrodynamics, Institutional Repository (TU
  Delft Library)

[Wiktorowicz & Laughlin(2014)]wiktorowicz2014toward
Wiktorowicz, S. J. & Laughlin, G. P. 2014, The Astrophysical Journal, 795, 12

[Williams(2017)]NASAfacts
Williams. 2017, NASA Planetary Fact Sheet - Metric, retrieved on 10/01/2017
  from <https://nssdc.gsfc.nasa.gov/planetary/factsheet/>

[Williams & Knacke(2004)]williams2004looking Williams, D. & Knacke, R. 2004, Astrobiology, 4, 400

[Wolszczan & Frail(1992)]WolszczanFrail1992 
         Wolszczan, A. & Frail, D. A. 1992, Nature, 355, 145




 






§ 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. <ref>): 
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. <ref>). 








 §.§ 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. <ref>). In principle its value ranges from 
0^∘ to 180, although the phase angle range accessible to an observer
depends on the inclination angle of the orbit of a body. A body in an edge–on
orbital geometry (inclination angle i= 90^∘) can attain phase angles between
0^∘ (when it is located behind the star) 
and 180^∘, while a body in a face–on orbital geometry 
(i=0^∘) can only be observed at α=90^∘. Generally, given
an orbital inclination angle i, the phase angle range is given by

    90^∘ - i ≤α≤ 90^∘ + i.
   
The phase angle of the planet or the moon at time t is computed as

    α_ x(t) = arccos[ 𝐮^T_ z·( - 𝐫_ xs(t)/𝐫_ xs(t)) ] ,

where subscript 'x' refers to either 'p' (planet) or 'm' (moon),
𝐮^T_ z= [ 0, 0, 1 ] 
is the unit vector along the z–axis, pointing towards the observer,
and vector r_ xs connects the center 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 the same as 
α_ m, yet our numerical model uses both values.




 §.§ Local viewing zenith angle θ_i


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

    θ_i = arccos[ 𝐮^T_ z·𝐫_i x/𝐫_i x] ,

where subscript 'x' refers to either 'p' or 'm', 𝐮_ z 
is the unit vector along the z–axis that points towards the observer, 
and 𝐫_i x is the vector pointing to the center of the 
pixel from the center of either the planet or the moon.
 



 §.§ 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 towards the star (see Fig. <ref>).
Angle θ_0i takes values between 0^∘ (at the sub-stellar 
location) and 90^∘ (at the terminator).
The position of the star changes in time, so that the time–dependent 
local illumination zenith angle can be computed as 

    θ_0i(t) = arccos[  𝐫^T_i x(t)/𝐫_i x(t)·( -𝐫_i s(t)/𝐫_i s(t))
                              ] ,

where subscript 'x' refers to either 'p' or 'm',  
𝐫_i x is the vector from the center of the 
planet or moon to the center of the pixel, and 
𝐫_i s the vector from the center of the star to 
the center of the pixel on the planet or moon. 





 §.§ 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 towards the observer and the plane 
described by the local zenith direction and the direction towards the 
star[Only the difference between ϕ_i 
and ϕ_0i is important, as our pixels are horizontally homogeneous].
As Fig. <ref> shows, ϕ_i - ϕ_0i follows 
from 

    ϕ_i-ϕ_0i(t) = arccos( 
                        cosα_i(t)-cosθ_icosθ_0i(t)/sinθ_icosθ_0i(t))  .

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 
𝐫_i x≪𝐫_ xs with 'x' 
referring to either 'p' or 'm', α_i can be approximated 
by the body's phase angle α_ x, and thus

    ϕ_i-ϕ_0i(t) = arccos(    
              cosα_ p(t)-cosθ_icosθ_0i(t)/sinθ_icosθ_0i(t))  .





 §.§ Local rotation angle β_i


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

    β_i = arcsiny_i px^2_i p+y^2_i p ,

where x_i p and y_i p are the coordinates of the center
of the pixel (recall that the z–axis points towards the observer).

For the lunar pixels, the alignment between the 
lunar scattering plane and the lunar grid, and hence angle β_i,
is time-dependent and requires 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:

    β_j(t) = arcsiny_j m(t)/x^2_j m(t) + y^2_j m(t).





 §.§ 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 the clock–wise direction from the lunar 
scattering plane to the planetary scattering plane.
For the results presented in this paper, the moon and the
planet orbit in the same, edge–on plane, and angle ψ equals zero. 
In the general case, however, it is computed using

    ψ(t) = arctan( -(𝐮^T_ y·𝐫_ ms(t)) / 
                                  (𝐮^T_ x·𝐫_ ms(t))  ) ,

where 𝐫_ ms is the vector from the star to the center of the 
moon, and 𝐮_ x and 𝐮_ y are the unit vectors 
along the x-axis and y-axis in coordinate system S_1.






§ 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, the size, 
shape, and depth of the shadow depend not only with 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 <cit.> 
for the Moon–Earth system, which we apply to our exoplanetary system. 
We model the umbral, antumbral, and penumbral shadow regions.
Figure <ref> shows the geometries involved in the 
various types of eclipses.
The equations used for computing the influence of eclipses are 
described here.

The flux arriving at a pixel i of eclipsed body B at time t depends 
on the fraction of the stellar disk and the local stellar surface brightness,
as seen from the center of the pixel. In Eq. <ref>, 
this is accounted for by factor c_i, the ratio between the actual 
flux e'_i on pixel i and e_i, the flux on the non–eclipsed pixel: 

    c_i(t) = e'_i(t) / e_i(t) = S'_Si(t) / S_Si(t),

with S'_Si and S_Si the stellar disk area as observed from 
pixel i when it is eclipsed and when it is non–eclipsed, respectively.
Note that we ignore stellar limb darkening and stellar light that travels 
through the atmosphere of the eclipsing body (if present).

To determine c_i, we first have to identify whether or not 
pixel i is eclipsed. Obviously, c_i=1 for a non–eclipsed pixel. 
If the pixel is (partly eclipsed), we have to determine the type of eclipse: 
umbral (i.e. total), antumbral, or annular. 
[A so–called hybrid eclipse
is an eclipse epoch where different types of eclipses occur along different
parts of the path of the shadow across the eclipsed body. Hybrid eclipses
are covered with our algorithms].
For an umbral eclipse c_i = 0.0, for an antumbral and
annular eclipse, we have to compute S'_Si in order to determine c_i. 


  


As can be seen in Fig. <ref>, 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:

    sinΩ = R_S+R_A𝐫_ AS .


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

    sinρ(t) < sinΩ(t),

with angle ρ ([0,90)) given by (see Fig. <ref>) 

    cosρ = -(𝐫_ AB-𝐀𝐎_4)
                  ·𝐮_ AS𝐫_ AB-𝐀𝐎_4 .


Vector 𝐀𝐎_4 is a function of the radii of the shadowed and 
eclipsing bodies and of Ω, as follows:

    𝐀𝐎_4 = R_B+R_AsinΩ 𝐮_ AS .


Except when body B falls completely in the penumbral cone, there will 
also be non–eclipsed pixels on the disk.
When Eq. <ref> holds, a
pixel-by-pixel search is performed, in which the center of each pixel is
checked for total or umbral (Sect. <ref>), 
annular (Sect. <ref>), or penumbral  
(Sect. <ref>) eclipse conditions.




 §.§ Total or umbral eclipses












In the umbral zone (see Fig. <ref>), pixels experience a 
total stellar eclipse. If the umbral zonal 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

    cosξ > cosΨ ,

with  

    sinΨ = R_S-R_A𝐫_ AS ,

and 

    cosξ = (𝐫_ AB-𝐀𝐎_3)·𝐮_ AS/𝐫_ AB-𝐀𝐎_3  with𝐀𝐎_3 = -R_A-R_BsinΨ 𝐮_ AS .


The disk of body B will only be partially inside the umbral shadow cone 
of body A when 

    cosΨ < (𝐫_ AB-𝐀𝐎_5)·𝐮_ AS𝐫_ AB-𝐀𝐎_5

and

    ( 𝐎_5𝐁 > R_BtanΨor𝐎_1𝐁 < R_B )

Here

    𝐎_5𝐁 = 𝐫_ AB-𝐀𝐎_5 and𝐎_1𝐁 = 𝐫_ AB-𝐀𝐎_1 ,
	
with

    𝐀𝐎_5 = -R_B+R_AsinΨ 𝐮_ AS and𝐀𝐎_1 = -R_AsinΨ 𝐮_ AS .
	

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

    sinβ_i < sinΨ .


Note that angle β_i of a pixel can be derived from 

    cosβ_i = (𝐫_ Ai-𝐀𝐎_1)·𝐮_ AS𝐫_ Ai-𝐀𝐎_1 .


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




 §.§ Annular or antumbral eclipses



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

    cosξ < -cosΨ ,

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

    cosβ_i < -cosΨ .


For each pixel in the antumbral cone, factor c_i is given by the fraction
of the stellar disk that is visible (see Fig. <ref>), i.e.
    c_i = πα^2_Si-πα^2_Aiπα^2_Si 
           = 1-(α_Aiα_Si)^2 ,

where

    α_Si = arcsinR_Si𝐫_Si  andα_Ai = arcsinR_Ai𝐫_Si .


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








 §.§ Penumbral eclipses



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

    cosω_i > cosΨ ,

where cosΨ follows from Eq. <ref> and 
cosω_i is given by:

    cosω_i = -(𝐫_ Ai-𝐀𝐎_2)·𝐮_ AS𝐫_ Ai-𝐀𝐎_2  with𝐀𝐎_2 = R_AsinΩ 𝐮_ AS .


The magnitude of the eclipse 
at pixel i, c_i, can be calculated through the stellar and 
eclipsing body viewing angles, i.e. the angular diameter of the 
bodies, 2α_S and 2α_A, and the eclipsing body-to-star 
angular distance as seen from the shadowed body, δ. Then, as 
follows from Fig. <ref>, c_i can be computed using

    c_i = πα_S^2-A_1-A_2πα_S^2 = 
            1- A_1+A_2πα_S^2,

with

    A_1 = {[     θ_A-sinθ_A2α_A^2 ifδ≥ l_S;                              ; π - θ_A-sinθ_A2α_A^2 ifδ< l_S ].

and

    A_2 = {[     θ_S-sinθ_S2α_S^2 ifδ≥ l_A;                              ; π - θ_S-sinθ_S2α_S^2 ifδ< l_A ].








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. <ref>) follow from e.g. Heron's formula:

    θ_S = 2arcsin( 2α_Sδ√(S(S-α_S)(S-α_A)
                (S-δ))) ,

and

    θ_A = 2arcsin( 2α_Aδ√(S(S-α_S)(S-α_A)
                (S-δ))) ,

with

    S = 12(α_A+α_S+δ) .


Pixels on body B that satisfy Eq. <ref> but do 
not meet the conditions established by Eqs. <ref>, 
<ref>, <ref> and 
<ref> are not eclipsed and have c_i=1.0.






§ 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, N_ eq^ p 
and N_ eq^ m, respectively, using a trade–off between the errors and 
the computational time, with time steps of 24 hours. 

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

    |Δ F_n(α)|     =    |F_n(α)-F_n-1(α)|F_n-1(α=0), 
           
        |Δ P_n(α)|     =     |P_n(α)-P_n-1(α)|,

with n-1 and n two consecutive values of N_ eq.

Figure <ref> 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. <ref>c. 
In the simulations presented in this paper, we decided to use 
N_ eq^ p= 50 and N_ eq^ m= 14, which yields an 
average disk–integration time of ∼ 0.8 seconds (thus an overall phase curve 
computation time of ∼ 2.4 minutes with a 24 h temporal resolution).
These values for N_ eq^ p and N_ eq^ m produce 
smooth curves for individual transit and eclipse events for temporal 
resolutions as small as 1 minute. 









