© The Authors 2025

1 Introduction

In the Solar System, polarimetric observations have been used for the characterization of planetary atmospheres (e.g., Dollfus & Coffeen 1970; Hansen & Hovenier 1974; West et al. 1983; Smith & Tomasko 1984; Tomasko & Doose 1984; Braak et al. 2002; Schmid et al. 2006, 2011; Rossi et al. 2015; McLean et al. 2017); circumplanetary rings (e.g., Dollfus 1979; Lane et al. 1982); and the surfaces of atmosphere-less bodies (e.g., Dollfus & Auriere 1974; Dollfus 1975; Kiselev & Lupishko 2004; Rosenbush & Kiselev 2005), for example for the surface of Earth’s moon (e.g., Gehrels et al. 1964; Coyne & Pellicori 1970; Dollfus & Bowell 1971; Shkuratov et al. 1992). Polarimetry is now becoming a valuable method for characterizing extrasolar planets as well. With modern polarimeters such as HIPPI (Bailey et al. 2015), HIPPI-2 (Bailey et al. 2020), or POLISH2 (Wiktorowicz et al. 2023), the polarization can be measured with a precision of a few parts per million (ppm). This places the detection of the polarized flux from hot Jupiters within the realm of possibility. For instance, Bott et al. (2018) set an upper limit of 40 ppm for the polarization amplitude of the reflected light signal of WASP-18b. Bailey et al. (2021) observed the linearly polarized flux of planet-hosting systems and detected the reflected light polarization signal of 51 Pegasi b with a 2.8σ significance.

The potential of polarimetric observations of the starlight scattered by Jupiter-like exoplanets (e.g., Seager et al. 2000; Stam et al. 2004) or Earth-like exoplanets (e.g., Stam 2008; Fauchez et al. 2017) has already been demonstrated in numerous theoretical studies using semi-analytical approaches or numerical simulations. This includes studies on the influence of circumplanetary rings (Lietzow & Wolf 2023; Veenstra et al. 2025), different cloud compositions (Karalidi et al. 2011, 2012; Rossi & Stam 2018; Lietzow & Wolf 2022), and the percentage of the planet covered in clouds (Rossi & Stam 2017; Winning et al. 2024). The dependence on the atmospheric composition enables polarimetry to determine whether a terrestrial exoplanet has an Earth-like or Venus-like climate (Mahapatra et al. 2023). The polarized glint of a surface ocean will also be detectable with spectropolarimetry using future instruments (Kopparla et al. 2018; Trees & Stam 2019, 2022; Vaughan et al. 2023).

Polarimeters with even higher accuracy will be required to use exoplanet polarimetry to its full capacity. With the Habitable Worlds Observatory (HWO), a coronagraphic space mission planned in the 2040s, it will be possible to directly image exoplanets in the habitable zone around their host stars and consequently to detect polarized scattered radiation from Earth-like exoplanets (National Academies of Sciences, Engineering, and Medicine 2023; Vaughan et al. 2023).

In the Solar System, six out of eight planets are orbited by moons. Thus, it is assumed that many exoplanets also host satellites. Despite the ongoing search for such exomoons (e.g., Kipping et al. 2012; Teachey et al. 2018; Fox & Wiegert 2021; Kipping et al. 2022; Kipping & Yahalomi 2023; Yahalomi et al. 2024), there has been no confirmed detection yet. Whether the transit signals associated with the Neptune-sized exomoon candidates around the planets Kepler-1625b (Teachey et al. 2018) and Kepler-1708b (Kipping et al. 2022) exist or are merely an artifact resulting from the data reduction, remains controversial (Kreidberg et al. 2019; Heller & Hippke 2024).

The inevitable influence of an exomoon on the flux and polarization measurements of its exoplanet has been addressed previously in a few studies. The flux change during mutual events, where one body casts a shadow on the other or transits in front of it, was calculated for the simple case of uniformly illuminated visible planetary and lunar disks by Limbach et al. (2024), for two Lambertian reflecting spheres by Cabrera & Schneider (2007), and for occultations in general by Luger et al. (2022). According to Limbach et al. (2024), the HWO will be able to detect mutual events in systems analogous to the Earth and its moon up to a distance of 10 pc. Sengupta & Marley (2016) modeled the influence of exomoons on the polarized thermal flux of self-luminous exoplanets during transits of the satellite in front of the planet. Berzosa Molina et al. (2018), on the other hand, used the code PyMieDAP (Rossi et al. 2018) to calculate the influence of an exomoon on the polarization of the radiation scattered by Earth-like exoplanets. The authors studied a system with a non-inclined exomoon orbit at a single wavelength of 450 nm. However, their set of mutual events is limited to five or six phase angles per type of mutual event.

For this article we acquired selected general results for an arbitrary exoplanet-exomoon system using analytic calculations. We also performed three-dimensional Monte Carlo radiative transfer simulations to model the polarimetric signature of mutual events in systems consisting of an Earth-like exoplanet and an exomoon similar to the Earth’s moon at wavelengths of 400 nm, 550 nm, and 700 nm as well as for the full range of possible phase angles in the case of a non-inclined exomoon orbit. An advantage of Monte Carlo simulations is that simplifications such as the locally plane-parallel atmosphere as used in Berzosa Molina et al. (2018) are not needed. We can therefore model mutual events at phase angles close to 180°, which require a fully spherical description of the planetary atmosphere. Additionally, we investigated the influence of noncentral transits and eclipses, that is when the moon is not in the orbital plane of the planet, but is slightly above or below it during the mutual event. This corresponds to possible geometries in systems with inclined exomoon orbits. We then evaluated the potential of polarimetric observations for the characterization of the properties of the exomoon including its orbital parameters.

In Sect. 2, we introduce the description of polarized radiation and our radiative transfer simulations. Thereafter, in Sect. 3 we describe the parameters of the model used for this study that resembles the Earth and its moon. In Sect. 4, we present the results. Section 4.1 focuses on the influence of an exomoon on the flux and polarization phase curve of an exoplanet neglecting the influence of mutual events, which is the topic of Sect. 4.2. The usability of polarimetric observations for a possible exomoon detection is discussed in Sect. 5. We conclude this study with a summary of our results in Sect. 6.

2 Methods

In this section we introduce the Stokes formalism and further concepts that are used to describe polarized radiation in this paper. We also present the framework of our simulations with the three-dimensional Monte Carlo radiative transfer code POLARIS¹ (Reissl et al. 2016). Recently, the code was extended to and optimized for the radiative transfer in exoplanetary atmospheres by Lietzow et al. (2021). For this study, we extended the code further to include an exomoon in the model space.

The state of polarization of electromagnetic radiation is described using the Stokes vector S = (I, Q, U, V)^T. Here, I is the total intensity, Q and U describe the linearly polarized intensity and polarization direction, and V is the circularly polarized intensity. Using the entries of the Stokes vector, the degree of linear polarization P and the angle of polarization χ are defined as (e.g., Newton 1966; Hansen & Travis 1974) $$P=\frac{\sqrt{Q^2+U^2}}{I},\tan 2\chi =\frac{U}{Q}.$$(1)

Throughout this study, the contribution of the stellar flux and thus its impact on the net polarization of a spatially unresolved exoplanetary system is not considered. However, the requirements for polarimetric observations with the star included are briefly discussed in Sect. 5.

POLARIS solves the radiative transfer equation by dividing the radiative field into photon packages which are characterized by their Stokes vector and move through the model space on a random path (Reissl et al. 2016). They are emitted at a random place on the surface of the star. The radiation emitted by the star is assumed to be unpolarized, as a disk-integrated polarization of less than 1 ppm was measured for the quiet sun (Kemp et al. 1987). Radiation emitted by other inactive FGK-dwarfs also has an intrinsic degree of linear polarization around only 10 ppm (Cotton et al. 2017). Because the exoplanet (in short “planet”) and its exomoon (in short “moon”) are seen under a very small solid angle from a point on the stellar surface (5.7 × 10⁻⁹ sr for the Earth seen from the sun), photons are only emitted toward the planet and its atmosphere as described in Lietzow et al. (2021) or toward the surface of the moon.

The modification of the polarization state and thus the values of the Stokes parameters of a photon package due to scattering by gas particles in the exoplanetary atmosphere is specified by the Müller matrix F(Θ). The Stokes vector S′ after scattering is (e.g., Newton 1966; Hansen & Travis 1974) $${\mathit{S}}^{\prime }\propto \mathbf{F}(\mathrm{\Theta })\cdot \mathit{S},$$(2) where Θ is the scattering angle and S is the Stokes vector of the incoming radiation. The influence of absorption by atmospheric particles is ignored in this study. For the atmosphere of Earth, absorption is of minor importance at wavelengths of 400 nm and 550 nm, while absorption by water vapor and O₂ has a noticeable influence on the flux at 700 nm (see, e.g., Schubert & Walterscheid 2000). Neglecting this influence has the advantage that the wavelength dependence of our results for an Earth-like exoplanet can also be interpreted as a dependence on the optical depth and hence on the density of the planetary atmosphere.

If a photon package is reflected by the surface of the planet or moon, its energy is weighted by the surface albedo ω which is the fraction of incoming radiation that is scattered by a surface rather than absorbed (e.g., Nicodemus et al. 1977). The corresponding impact on the Stokes vector is described by a reflection matrix (see, e.g., Zhai et al. 2010).

3 Model of planet and moon

For our simulations, we chose a model consisting of an Earth-like exoplanet and an exomoon comparable to the moon of Earth. The reasoning behind this choice is that such a system has a much more favorable moon-to-planet continuum flux ratio of F_m/F_p ≈ 2% than the moons around the gas planets in the outer Solar System, which reach F_m/F_p of at most 10⁻³, with the only exception of the Neptune-Triton pair with 0.5% (see Appendix A). The parameters of the model are listed in Table 1. The planetary parameters were chosen to be comparable to the model used by Stam (2008) and Berzosa Molina et al. (2018).

The planetary orbit around the star in the model is circular (e_b = 0) with a radius of d_sp = 1 AU. It is seen edge-on (i_b = 90°). We neglect the distance between the planet–moon barycenter and the geometric center of the planet. The orbit of the moon is also circular, with a radius of d_pm = 380 000 km. In the baseline model, it is not inclined with respect to the orbital plane of the planet. The influence of an inclined orbit of the moon is topic of Sect. 4.2.4. The position of the planet with respect to the star and the observer is defined using the planetary phase angle α_p. The angle φ_ms specifies the position of the moon on its orbit around the planet as defined in Fig. 1. The moon orbits the planet prograde with respect to the planetary orbit.

The model atmosphere of the planet in our simulations consists solely of gas particles and does not include a cloud layer. Because gas particles are small compared to the considered wavelengths (400 nm, 550 nm, 700 nm), we describe their scattering properties by applying the concept of Rayleigh scattering. The scattering matrix F(Θ) for Rayleigh scattering can be found in Hansen & Travis (1974). The Rayleigh scattering cross section for anisotropic, randomly oriented gas particles at wavelength λ is given by Sneep & Ubachs (2005).

We use typical values for air and assume that the depolarization factor δ = 0.028 (Bates 1984) is constant across all considered wavelengths. The refractive index n(λ) is calculated according to Keady & Kilcrease (2000) as $$n(\lambda )-1=A\left(1+\frac{B}{{(\lambda /\text{μm})}^2}\right),$$(3) in which for air A = 28.71 × 10⁻⁵ and B = 5.67 × 10⁻⁵. The model atmosphere consists of 32 vertical layers with a total height of 100 km. We use the pressure profile for the mid-latitude summer profile of the atmosphere of Earth as given in McClatchey et al. (1972) and calculate the number density in each layer assuming hydrostatic equilibrium (similar to Berzosa Molina et al. 2018; Lietzow et al. 2021) and a molar mass of air of 28.964 × 10⁻³ kg mol⁻¹ (Schubert & Walterscheid 2000). The optical depth of the model atmosphere is τ_λ = 0.098 at a wavelength of 550 nm, identical to the result by Stam (2008).

Both bodies have diffuse reflecting (Lambertian) surfaces with albedoes of ω_p = 0.3 and ω_m = 0.1 for the planet and the moon, respectively. After a reflection by such a surface, the radiation is completely depolarized (see also Lietzow et al. 2021). Diffuse reflecting surfaces are commonly used as a model for exoplanetary surfaces (e.g., Stam 2008; Rossi & Stam 2017; Berzosa Molina et al. 2018). However, they are only an approximation of real planetary surfaces. Even radiation reflected by atmosphere-less bodies such as Mercury or Earth’s moon reaches degrees of linear polarization of a few percent at wavelengths in the visual spectrum (e.g., Lyot 1929; Coyne & Pellicori 1970; Dollfus & Auriere 1974).

Because the moon in the model is a diffuse reflecting sphere, its flux phase curve is calculated analytically according to Russell (1916) in the following sections. Because it does not have an atmosphere, it does not contribute to the polarized flux. In our study, we evaluate in which cases this assumption is justified and present analytical estimates of the systematical errors.

Fig. 1 Geometry of the non-inclined baseline system. The centers of the star, the planet, and the moon as well as the observer are all coplanar. The angles αp and αm are the phase angles of the planet and the moon, respectively. φms is the angle between the direction to the star and the direction to the moon for an observer at the center of the planet.

4 Results

In the following presentations of our modeling results, we normalize the flux of the planet F_p and the moon F_m such that they are equal to the geometrical albedoes A_p and A_m at a phase angle of 0°. For the planet, this is achieved by dividing the planetary flux phase curve F_p(α_p) by F_0p, the flux reflected by a flat fully diffuse reflecting disk with the same radius as the planet. If B_ν is the Planck function, ν the frequency, T_s the stellar effective temperature, r_s the stellar radius, and d_po the distance between the planet and the observer, then $$F_{0\mathrm{p}}=\pi B_{\nu }(T_{\mathrm{s}})\frac{r_{\mathrm{s}}^2{r}_{\mathrm{p}}^2}{d_{\mathrm{sp}}^2{d}_{\mathrm{po}}^2}.$$(4)

This is a common approach to make the results independent of the basic stellar and planetary parameters (e.g., Stam 2008; Berzosa Molina et al. 2018; Lietzow et al. 2021). For the exomoon F_0m is defined accordingly. The flux of the unresolved planet-moon system is normalized by dividing by F_0p which facilitates the comparability with the results for the flux of the planet alone.

Fig. 2 Phase curves of the planet, the moon, and the planet and moon together but without the influence of mutual events. Top: Flux F/F0 as a function of the phase angle α. For the planet and for planet and moon combined, F0 = F0p. In contrast, for the moon, F0 = F0m. In addition, we use the analytic result from Russell (1916) for a Lambertian sphere for the phase curve of the moon. Middle: Phase curves for the linear polarization P(α) for the planet and planet and moon combined, assuming Pm = 0 (i.e., of the diffuse reflecting moon surface). The stellar flux is not included in the calculation of P. Bottom: Difference P − Pp between the polarization phase curve of the unresolved system consisting of planet and moon and the planetary polarization phase curve.

4.1 Phase curves

First, we calculated phase curves of the flux and degree of linear polarization of the moon, the planet, and the two bodies combined neglecting the influence of possible eclipses, shadows, transits, or occultations which are discussed in detail in Sect. 4.2. Solely for the curves presented in this subsection, the Stokes vector of the planet S_p and the moon S_m were calculated separately and added to get the Stokes vector of the unresolved system. The systematic error caused by neglecting the flux that an observer measures due to the mutual illumination of the planet and the moon is less than 2 × 10⁻⁴ of the maximum lunar flux for the model system (see Appendix B). Also, we neglected the small difference in the planetary and the lunar phase angle and assumed α_p ≈ α_m ≈ α. For the planetary Stokes vector, we used results of POLARIS simulations. The disk-integrated results of our simulations are shown in Fig. 2.

The linear polarization phase curve of the planet resembles a typical curve for planets with Rayleigh scattering atmospheres, which is discussed in detail by Stam (2008). The maximum degree of polarization P of 42.4% at a wavelength of 400 nm is reached at a phase angle of 96°. For the longer wavelengths, the maximum shifts to larger phase angles of 109° for 550 nm, where P reaches 22.9%, and 120° for 700 nm, where P is 11.9%.

The direction of polarization is mostly perpendicular to the orbital and scattering plane of the planet. Only at very large phase angles, it is parallel to the scattering plane. This results in a zero-crossing of P(α) when the direction of polarization changes. For 400 nm this change occurs at 166° and for 550 nm at 171°. According to Stam (2008), this change is due to a greater influence of multiple scattered photons. For 700 nm the zero-crossing and the change in the polarization direction vanish, as the probability of multiple scattering is too low due to the small atmospheric optical depth of only 0.037.

The moon contributes to the total flux but not to the polarized flux of the system. Therefore, the degree of linear polarization is slightly lower at all phase angles and wavelengths for the planetmoon system compared to the planet, similar to the results by Berzosa Molina et al. (2018) at 450 nm. The maximum reduction for 400 nm is found at a phase angle of 91° with a difference of −0.77%. For longer wavelengths, the maximum reduction shifts toward larger phase angles, similar to the maximum of P(α). For 550 nm the maximum reduction is found at a phase angle of 101° with a difference of −0.46%, for 700 nm it is at α = 112° with −0.25%. The phase angle dependent influence is also illustrated in the bottom panel of Fig. 2.

We find that the behavior might be converse for a real exomoon though. Let F_p, P_p be the flux and degree of polarization of the planet, respectively, and F_m, P_m the same quantities for the moon. If the planet and the moon are symmetric to the scattering plane, the direction of polarization is either perpendicular or parallel to the scattering plane. Therefore, the polarized fluxes F_pP_p and F_mP_m either have to be added or subtracted to get the polarized flux of the unresolved system. Because the planetary flux F_p is much bigger than the flux of the moon F_m, we assume that F_pP_p ≥ F_mP_m. For the model system, this assumption is fulfilled everywhere as P_m = 0. The difference between the degree of linear polarization of the unresolved system and the planet alone is therefore $$P-P_{\mathrm{p}}=\frac{F_{\mathrm{p}}{P}_{\mathrm{p}}\pm F_{\mathrm{m}}{P}_{\mathrm{m}}}{F_{\mathrm{p}}+F_{\mathrm{m}}}-P_{\mathrm{p}}=-\frac{F_{\mathrm{m}}}{F_{\mathrm{p}}+F_{\mathrm{m}}}(P_{\mathrm{p}}\mp P_{\mathrm{m}}).$$(5)

Here, the minus-plus-symbol (∓) in the last expression differentiates between cases in which the direction of polarization is identical for the planet and moon (−) and cases in which they are perpendicular (+). If P_m > P_p while the directions of polarization are parallel, the quantity P − P_p is therefore positive.

The direction of the polarization of Earth’s moon is parallel to the scattering plane for phase angles of less than the inversion phase angle α_inv, which is between α = 20° and 30°, and perpendicular for phase angles larger than α_inv (e.g., Lyot 1929; Coyne & Pellicori 1970; Shkuratov et al. 1992). If P_m > 0 is assumed, we thus expect a larger reduction of P than shown in Fig. 2 for phase angles below α_inv and at large phase angles where the planetary direction of polarization is parallel to the scattering plane. Near quadrature, where we found the maximum differences P − P_p due to a moon, the degree of linear polarization is less affected. Assuming a realistic value of P_m ≈ 10% (Lyot 1929; Coyne & Pellicori 1970), the maximum value of P − P_p at 550 nm is only −0.25% rather than −0.46% at 101°.

Even if it were possible to measure flux and polarization with a precision high enough to detect the subtle differences between the curves of the planet alone and of the planet-moon system, the following ambiguity would remain. A slightly higher flux and a slightly altered degree of linear polarization can also be explained by an adapted planetary model, for instance by a slightly higher surface albedo. Therefore, the presence of an exomoon cannot be deduced from the phase curves shown in Fig. 2 (see also Berzosa Molina et al. 2018). This ambiguity can be resolved by observing the influences of mutual events, which are investigated in the next section.

Fig. 3 Four types of mutual events simulated for a phase angle of 225°. The RGB images are based on individual maps simulated for wavelengths of 440 nm (blue), 555 nm (green), and 700 nm (red). The total flux was normalized to the geometric albedo separately at each of the wavelengths. The dashed circles indicate the position and size of the planet (blue) and moon (gray). To improve the visibility of the moon and to fit both bodies into the images for shadow and eclipse, the lunar surface albedo am was increased to 0.3 and the moon’s semimajor axis was set to 2.8 rp for this graphic.

4.2 Mutual events

The four types of mutual events that can occur in a system are shown in Fig. 3. We stick to the terminology introduced by Cabrera & Schneider (2007) and Limbach et al. (2024):

1. The moon passes between the planet and the star. It casts a shadow on the planet.

2. The moon passes between the planet and the observer, which is called a transit of the moon in front of the planet.

3. The moon is in the shadow of the planet, a lunar eclipse.

4. The planet passes between the moon and the observer. This is called an occultation of the moon by the planet.

We first calculate the influence of full occultations of an exomoon and total lunar eclipses in Sect. 4.2.1 analytically. Thereafter, the influence of transits and shadows is determined for the model system described in Sect. 3 using POLARIS simulations in Sects. 4.2.2 and 4.2.3. Finally, the differences for noncentral transits and shadows are highlighted in Sect. 4.2.4. We are interested in the change in the total flux △F and the linear polarization △P, which are $$\mathrm{\Delta }F=F_{\text{with}}-F,\mathrm{\Delta }P=P_{\text{with}}-P,$$(6) where F_with and P_with are respectively the flux and linear polarization observed during a mutual event, and F and P are the same quantities that would be observed at the same phase angle without the mutual event.

4.2.1 Full occultations and eclipses

If the moon is completely occulted by the planet or completely within the umbra of the planet, the flux change is △F = −F_m, where F_m is the lunar flux that would be measured at the respective phase angle if there was not an eclipse or occultation. Here, we neglected any radiation that is scattered toward the moon in the planetary atmosphere during a lunar eclipse (Cabrera & Schneider 2007; Limbach et al. 2024).

Because the influence of the moon vanishes during these two mutual events, the change in the degree of polarization △P during the event is the inversion of the influence described in Eq. (5). For the model system, the maximum change during these events is therefore extracted from the bottom panel of Fig. 2 by flipping the algebraic sign. We find the maximum changes to be △P = +0.77% for 400 nm if a mutual event happens at α = 91°, △P = +0.46% for 550 nm at α = 101°, and △P = +0.25% for 700 nm at α = 112°.

Assuming reflection properties similar to the Earth’s moon, the same influence that was described in the previous section is expected. Assuming P_m ≈ 10% around quadrature, the maximum △P at 550 nm is only +0.25%. The systematic error caused by the assumption of P_m = 0 is on the same order of magnitude as the actual result and △P might even be negative if P_m > P_p.

It should be noted that we calculated much smaller values for △P than reported by Berzosa Molina et al. (2018) who found maximum values for △P of +2.5% during full occultations of the moon and of +2.7% during lunar eclipses. However, using values for P_p, P_m, △F, and F_p given in Berzosa Molina et al. (2018) and Eq. (5) results in a realistic △P of approximately 0.7% to 0.8%, which is similar to the value we determined for 400 nm.

For the lunar darkening events, △P is sensitive to the surface albedo of the moon. If P_m = 0 and the albedo of the moon is ω_m > 0.2, the changes △P reach values of 3% and more.

According to Eq. (5), observations during an eclipse or occultation allow for the determination of the difference or sum of the degrees of polarization P_p ∓ P_m, because F_m = △F and F_p = F_with can be measured directly. P_p is also measured directly during a transit because the degree of polarization of the unresolved system is P = P_p when the moon is completely eclipsed or occulted. Therefore, the sign ∓ is also determined. As a consequence, both P_p and P_m can be determined at the phase angle where the eclipse or occultation is observed. In other words, such observations allow us to disentangle the planetary and the lunar polarization phase curve.

4.2.2 Central transits

We now continue with a more in-depth description of mutual events which influence the planetary flux. In this section, we start with central transits of the moon in front of the planet.

In the model system, such a central transit occurs for angles α_p + φ_ms in the interval [360° − △φ_ms, 360° + △φ_ms]. For the model system, the half width △φ_ms is 1.23° (see Appendix C). A transit in the model system therefore lasts no longer than 4.4 hours. With increasing phase angle, a decreasing part of the day side of the planet is visible to the observer. The time during which a transit causes a change in flux or linear polarization is therefore decreasing with increasing phase angle (for an analytical approach, see Cabrera & Schneider 2007; Limbach et al. 2024).

The values of △F/F and △P resulting from our simulations are visualized in Fig. 4. As △F is strongly dependent on the reflected flux phase curve of the planet as shown by Berzosa Molina et al. (2018), the quantity △F/F during a central transit is more insightful. The color maps have a resolution of 181 phase angles α_p between 0° and 180° and a resolution of 58 lunar position angles φ_ms for each α_p. The values for △F/F and △P for a single transit are found along a passage perpendicular to the abscissa through the respective diagrams as the planetary phase angle only changes about 0.2° during a transit. For the model system, transits occur 6 or 7 times during half the revolution of the planet around the star.

For phase angles below 2°, we find an obvious feature in the △F/F plots. If a transit happens at these phase angles, the shadow of the moon also falls on the planet. Disregarding this special configuration with a doubled influence, the maximum changes △F/F for phase angles up to 30° are between −8% and −12%. With increasing phase angle, central transits cause larger reductions in △F/F. During a transit at higher phase angles, the moon covers a bigger fraction of the illuminated part of the planetary disk that is now restricted to a small crescent. At a wavelength of 400 nm, △F/F reaches the biggest decrease of −32.8% for a transit at α_p = 144°. This maximum shifts to smaller phase angles for the other two wavelengths, reaching values of −32.3% at 132° at 550 nm and −33.9% at 131° at 700 nm.

The similarities between the different wavelengths are to be expected, given that the flux is mainly influenced by geometric effects, for instance, what fraction of the planetary day side is occulted. The higher decreases of △F/F during a transit toward longer wavelengths at phase angles below 60° are caused by the decrease in optical depth. The planetary brightness distribution approaches that of a diffuse reflecting sphere, so the flux is more concentrated toward the sub-stellar point on the planetary surface. The moon thus blocks a higher fraction of the planetary flux when it occults the sub-stellar point.

At large phase angles above 160°, a light ring forms around the planet as a consequence of multiple and forward scattering in the planetary atmosphere, an additional contribution to F, that decreases the influence △F/F of the lunar transit. Also, a flux change between −2% and −7% occurs before the moon transits the very small crescent on the day side of the planet starting at a phase angle of 176°. During this phase, the moon is in front of the bright atmosphere on the side opposite the crescent. We note that this light ring only occurs in a simulation when a fully spherical model of the planetary atmosphere is used.

The absolute change △P in the degree of linear polarization during central transits inhabits a slightly higher dependence on the optical depth and hence on the wavelength. Furthermore, the dependence on the planetary phase angle α_p is more complicated. Multiple different regions are found and decreases as well as increases in P are possible.

In Fig. 5, the polarized flux F_pP_p and the direction of polarization χ are shown for the resolved planet using results from the simulations in Sect. 4.1 for a wavelength of 400 nm. Similar figures for 550 nm and 700 nm can be found in Appendix E. Figure 5 is used to facilitate the following discussion.

At a wavelength of 400 nm, a slight increase of △P ≈ 0.2% is observed at phase angles up to 20° when the moon is in front of the edge of the planet as seen from the observer. In between, P returns to values close to those calculated without a mutual event. From the panel for a phase angle of 5° in the first row of Fig. 5, it is clear that the polarization directions form an almost symmetrical pattern. Consequently, P_p(5°) is small at only 0.2% at a wavelength of 400 nm and even lower for 550 nm (0.1%) and 700 nm (0.04%). As described by Berzosa Molina et al. (2018), the moon breaks this symmetry when it crosses the limb of the planet as sketched in Fig. 5.

We find that the described feature is too faint at 550 nm and 700 nm to appear in the maps for △P in Fig. 4. For phase angles below 10°, △P is below 0.1% because the polarized flux is too small even if the symmetry is broken. For larger phase angles, due to the smaller influence of multiple scattering at lower optical depths, the direction of polarization is perpendicular to the scattering plane across the planetary disk at these wavelengths which causes a different behavior.

At a wavelength of 400 nm, a direction of polarization vertical to the scattering plane throughout the planetary disk is observed for phase angles above 20°. In Fig. 5, this can be seen in the surface brightness maps of the disk-resolved polarized flux of the planet at phase angles of 30°, 75°, and 120°. The moon now causes a central peak in the degree of polarization. In addition, a slight decrease in P is observed at the beginning and the end of a transit at phase angles between 50° and 120°. These features are also present in the phase curves in Berzosa Molina et al. (2018), who find maximum absolute values at 450 nm that are similar to our results at 400 nm and 550 nm.

Because the polarization directions are parallel, the change △P depends solely on the fraction of blocked polarized flux compared to the fraction of blocked total flux. The central peak in △P during the transit occurs when the moon is in between the terminator and the edge of the planet. Here, it cannot block a high fraction of F_pP_p. Due to the higher fraction of F being blocked, P increases. At 400 nm, the largest change is an increase of △P ≈ 2.6% around phase angles of α_p ≈ 98°. With increasing wavelength, thus decreasing optical depth, we find that the maximum value in △P shifts toward larger phase angles reaching 2.7% during a transit at a planetary phase angle of 102° for a wavelength of 550 nm and 1.9% for a planetary phase angle of 116° at 700 nm. This is similar to the maximum of the planetary polarization curve that shifts toward larger phase angles with increasing wavelength as well (see Fig. 2).

Interestingly, the change in polarization during transits of the moon in front of the planet is relatively robust with respect to the actual lunar polarization phase curve (see Appendix D). Although the assumption that P_m = 0 may result in a slight underestimation of positive changes △P and a slight overestimation of negative changes △P, it does not change the general behavior shown in Fig. 4. Using a realistic maximum value of P_m = 10% (Lyot 1929; Coyne & Pellicori 1970), the expected absolute systematic errors △P_sys for the maximum changes △P given above are approximately △P_sys ≈ 0.1% for all wavelengths.

At α_p = 75°, F_pP_p is highest at the edge of the visible planetary disk. When the moon transits across this edge, P decreases because a higher fraction of F_pP_p than of F is blocked by the moon. Near the terminator, F_pP_p and F are both small. However, again a higher fraction of F_pP_p than of F is blocked and the degree of polarization decreases (for an explanation using the Stokes parameters, see Berzosa Molina et al. 2018). The biggest decrease is higher than △P = −0.1% at all wavelengths but has a slightly smaller absolute value than −0.1% at the end of transits at 700 nm as P is lower in general.

We investigate the reason why the decreases before and after the transit only occur in the mentioned phase angle range. At α_p = 30°, F_pP_p is concentrated in two regions which are on the edge of the visible planetary disk but not in the orbital plane and thus outside the areas the moon can occult. As a consequence, a decrease in P at the beginning and end of a transit does not occur. For larger phase angles of more than 120°, the decrease is not observed either. The moon still blocks a high fraction of F_pP_p when it crosses the planet’s edge. However, this fraction of F_pP_p now corresponds to a large fraction of F as well, as the scattered radiation flux of the planet is limited to a relatively small crescent. F_pP_p is even more concentrated on the edge of the planetary disk for wavelengths of 550 nm and 700 nm than for 400 nm. Consequently, a decrease in P before and after the transit is observed over a larger range of phase angles.

We find another change in the behavior when the diskintegrated direction of polarization of the radiation scattered by the planet changes at α_p = 166° for 400 nm or at α_p = 171° for 550 nm. At a wavelength of 700 nm, the direction of polarization does not change at large phase angles and consequently, the corresponding feature disappears completely. The degree of linear polarization now decreases by at most 0.75% during a transit.

As an example, a polarization map at 172° is shown in the bottom row of Fig. 5. F_pP_p is largest at the illuminated crescent and the polarization directions are vertical to the planetary edge. In the area, where the moon occults the crescent, χ is parallel to the orbital plane. For polarization vectors on the crescent in a greater distance of the orbital plane, we always find a vector in a region of approximately equal F_pP_p that lies perpendicular to the first one. Therefore, these areas do not contribute to the diskintegrated polarized flux. Around the poles of the planet, F_pP_p is too low to contribute significantly. The moon therefore covers the part of the crescent that causes a disk-integrated degree of polarization, so P decreases significantly.

This feature is much smaller at a wavelength of 550 nm because the zero-crossing of the planetary polarization phase curve occurs at a larger phase angle. As for a wavelength of 700 nm, we instead find a direct transition into the final behavior close to a planetary phase angle of 180° which happens at 178° for a wavelength of 400 nm: P is close to 0 again due to the symmetry of the system. During a transit, P now increases again. The moon again breaks the symmetry during the transit.

Fig. 4 Relative change of the total flux △F/F (top row) and absolute change of the degree of polarization △P (bottom row) during central transits of the moon in front of the planet at different phase angles αp visualized using color maps and contour lines. The upper abscissa shows the time in days since the planet was at α = 0°. On the left ordinate, the relevant angle αp + φms is shown. The angle was transformed to a time in hours on the right ordinate, where a time of 0 corresponds to the time when the planet, moon, and observer are perfectly aligned. A single transit of the moon corresponds to a passage through the color maps perpendicular to the abscissa. For phase angles above 180°, F(360° − αp, 360° − φms) = F(αp, φms).

Fig. 5 Polarization maps for different phase angles at a wavelength of 400 nm. The polarized flux FpPp is encoded as a color map. The resolution is 60 × 60 px, the planet has a diameter of 50 px. Polarization directions are shown as white lines. The larger orange circle marks the planet’s circumference. The orange dashed lines between the poles of the planet mark the terminator. The day side of the planet is to the right of it. The horizontal orange lines mark the path of the moon during a central transit (lower pair of lines), or a noncentral transit (upper pair of lines). For two positions on the edge of the planet, the circumference of the moon is shown as an orange circle.

4.2.3 Central shadows

The moon in the model system in which all bodies reside within a single plane casts a shadow on the planet if φ_ms ≈ 0. The half width of a shadow phase △φ_ms is the maximum |φ_ms| for which part of the planet is within the lunar penumbra. For the model system, we find an upper limit △φ_ms < 1.49° (see Appendix C). A shadow phase in the model system is therefore no longer than 5.4 hours. Similar to transits, the duration of a shadow phase decreases with increasing phase angle (for an analytical approach, see Cabrera & Schneider 2007; Limbach et al. 2024).

Figure 6 shows the relative flux change △F/F during a shadow phase and the absolute change in the degree of linear polarization △P. The color maps have a resolution of 181 phase angles between 0° and 180° and a resolution of 64 lunar position angles φ_ms between ±1.49°. The influence of a single shadow phase is found along an almost vertical line through the diagram. For the model system, a shadow phase occurs 6 to 7 times during half the planetary revolution around the central star. As expected, the influence of shadows cast on the planet by the moon resembles the influence of lunar transits in front of the planet. We therefore focus on the observed differences.

The maximum values of the relative flux change △F/F during a shadow phase are smaller than for transits at all wavelengths and phase angles. Because the umbra of the moon is small compared to its penumbra, the darkening of a region on the planet by the shadow of its moon causes a smaller △F than the occultation of the same region during a transit (Berzosa Molina et al. 2018). The largest reductions are −26.8% found at 146° for a wavelength of 400 nm, −25.9% at 141° for 550 nm and −26.9% at 133° for 700 nm. Again, the influence of the optical depth or the wavelength on △F/F is small.

For α_p ≤ 2°, the additional influence of a transit is visible again, this time causing a reduction in △F/F at the beginning of the eclipse that corresponds to the end of the transit. As for transits, the maximum values △F/F reached during a shadow phase increase with increasing phase angle because the illuminated crescent of the planetary surface is smaller at longer wavelengths. Starting at phase angles of α_p ≈ 140°, △F/F decreases again due to the influence of multiple and forward scattering. The small flux changes of about 2% to 5% that occur at phase angles larger than 176° after the shadow phase ended are caused by the shadow falling on the atmosphere on the side of the planet opposite to the crescent.

The changes △P in Fig. 6 show features that are related to the ones found for transits in Fig. 4. The explanation for transits applies analogously to shadow phases. As the decrease in the flux is smaller for shadow phases, the positive changes in the degree of linear polarization are also smaller, because they are mainly the consequence of a decrease in F. For a wavelength of 400 nm, the maximum increase in P is +1.8% for a shadow at a phase angle of 96°. Berzosa Molina et al. (2018) find a similar value at 450 nm. At the other wavelengths, we find +1.8% at 106° at 550 nm and +1.3% at 116° at 700 nm. As for transits, the phase angle with the maximum increase shifts toward larger phase angles with increasing wavelength.

Contrary to transits, a slight decrease in P is only observed when the shadow is on the illuminated crescent but not at the beginning of a shadow phase. A shadow never falls on the planetary night side, and its shape that appears to an observer is not circular when it first reaches the planetary surface at the terminator. Together, this prevents a decrease in P.

Similar to transits, a decrease in P of up to 0.7% is observed at a wavelength of 400 nm after the direction of polarization of the planet changes. The feature again disappears at a wavelength of 700 nm and occurs in a very small range of phase angles for a wavelength of 550 nm. These changes appear to be slightly more pronounced for central shadows than for central transits.

Fig. 6 Similar to Fig. 4, but for phases during which the moon casts a shadow on the planet. The left ordinate now shows the angle φms, which specifies the lunar position. Again, this angle was transformed to a time in hours, as shown on the right ordinate, where a time of 0 corresponds to the time when the planet, moon, and star are perfectly aligned. A single run of the shadow of the moon across the planetary surface corresponds to a passage through the color maps perpendicular to the abscissa.

4.2.4 Noncentral transits and shadows

In Sects. 4.2.2 and 4.2.3, we have shown that △P depends on the exact part of the planet that is occulted by the moon or in its shadow. Only if the moon is within the planetary orbital plane, △P takes the values found in the previous sections.

If i_b = 90° and i_m ≠ 0, the moon may also transit the planet along a noncentral path and similarly its shadow may take a noncentral path on the planetary surface. We mimic the effect of a nonzero inclination of the lunar orbit by placing the moon in a distance of r_p/2 above the planetary orbital plane in every simulation while keeping its distance to the planet d_pm constant. In the case of the model system, such a configuration is possible if i_m > 0.48°, a small value compared to the Earth’s moon with i_m ≈ 5°. For the observer, the moon now crosses the planetary disk above the planetary orbital plane, as sketched in Fig. 5.

To specify the occurrences of transiting and shadowing events, we define φ′ as the angle between the direction from the center of the planet toward the star and the direction from the center of the planet toward a projection of the center of the moon on to the planetary orbit plane. Now a noncentral transit occurs if α_p + φ′ ≈ 360° and a noncentral shadow falls on the planet if ${\phi }_{\text{ms}}^{\prime }\approx 0$.

For a wavelength of 400 nm only, we repeated the simulations for transits and shadows as described in Sects. 4.2.2 and 4.2.3, respectively, with the repositioned moon. Our results for the relative flux change △F/F and the absolute change in polarization △P for noncentral transits are shown in Fig. 7.

Due to the change in the geometry, noncentral transits last shorter than central ones. The maximum decrease of the flux is slightly smaller than for central transits. Values of −25.2% are reached for a transit at a phase angle of 141°. However, the shape of the observed △F/F curves remains unchanged compared to central transits. A moon with a smaller radius on a slightly narrower orbit around the planet can therefore cause the same changes in the observed total flux.

However, the changes △P in the degree of linear polarization behave differently for a noncentral transit, as we see from a comparison of Figs. 4 and 7. While for central transits, a maximum increase of 2.6% is found at a phase angle of 98°, the change in the degree of linear polarization is less than 0.5% at the same phase angle for a noncentral transit. A decrease in P is observed at the end and the beginning of the transit over a larger range of phase angles. Therefore, △P inhabits a sequence of a decrease, an increase, and a decrease, all with absolute values of △P between 0.1% and 0.5%, during noncentral transits at phase angles around 90°. This wave-like behavior of △P is different from the pronounced central peak found for central transits at the same phase angles, where △P reached +2.5%.

For the specific type of noncentral transits described above, the highest increase △P of about 1.5% is observed at a phase angle of 166°. This is the phase angle at which the direction of polarization changes (P is equal to zero for the model system at this phase angle). For central transits, we observe very small changes at α_p = 166°. Also, a decrease in the degree of polarization is not found at phase angles above 166°. Instead, P increases slightly during a noncentral transit at those phase angles.

The new behavior outlined in this section is explained using Fig. 5 again. For a noncentral transit at a phase angle of 75°, the moon occults parts of the edge of the visible planetary disk for a longer period of time, which causes the decrease in polarization that is connected to this to last longer as well. Additionally, the moon never covers the region in the orbital plane between the terminator and the edge. The occultation of this region is however connected to the highest increases for central transits around 100°. Therefore, smaller maximum values are observed for noncentral transits at these phase angles.

In the case of a phase angle of 166°, the polarized flux is slightly smaller in the orbital plane than where the moon crosses the planetary edge. The net polarization is zero without the transit because the equivalent region on the lower half of the planet has a similar F_pP_p and directions of polarization perpendicular to the ones on the upper half. However, during the transit, the region on the lower half of the planet contributes to the disk-integrated polarized flux. Consequently, P increases considerably. The direction of polarization is diagonal to the orbital plane during such a transit. For a similar reason, P increases for noncentral transits at higher phase angles. As shown in the panel for 172° in Fig. 5, the exomoon does not cover the part of the crescent in the orbital plane that causes the disk-integrated degree of polarization if it transits noncentrally at this phase angle. Instead, it breaks the symmetry between the upper and lower half of the visible planetary disk.

Similar to central transits and shadows, noncentral transits and shadows also cause connatural results, as is clear from our results for noncentral shadows shown in Fig. 8. The decrease in the degree of linear polarization at the beginning of a noncentral shadow phase is again much smaller than for a transit because the shadow never falls on the night side of the planet and thus never crosses the terminator. The decrease at the end of a shadow phase is observed for a larger range of phase angles. For a phase angle around 90°, the resulting changes in △P are therefore an increase during the first half of the event followed by a decrease in the second half, with similar maximum absolute values of △P during both halves.

Clearly, the behavior of △P depends on the exact path of the moon or its shadow across the visible planetary disk. Mutual events caused by a moon can thus show large variations in their influence on the linear polarization, even when two events in the same system are compared.

Fig. 7 Same as the first column of Fig. 4, but for the results of our simulations with noncentral transits.

Fig. 8 Same as the first column of Fig. 6, but for the results of our simulations with noncentral shadow phases.

5 Discussion

5.1 Resolved and unresolved observability

So far, we have only considered the flux of the planet and the moon but not the stellar flux. This resembles observations during which the planet and the star are spatially resolved. Such observations of the scattered radiation of exoplanets are rarely possible today. Until now, directly imaged planets are either massive gas giants in large separations from the host star (e.g., β Pic b, Lagrange et al. 2009) or accreting planets still within a protoplanetary disk that were discovered through detections of their thermal emission (e.g., PDS 70 b, Keppler et al. 2018). Massive gas giants in large separations from the host star are likely a promising target to discover exomoons (see also Lazzoni et al. 2022). However, high-precision polarimetry is only promising for close-in exoplanets at the moment.

Observing the scattered radiation of Earth-like exoplanets directly is a target of the Habitable Worlds Observatory, a future space telescope in the 2040s, which will be equipped with a state-of-the-art coronagraph (National Academies of Sciences, Engineering, and Medicine 2023; Vaughan et al. 2023; Tuchow et al. 2024). The HWO will be able to detect mutual events in systems analogous to the Earth and its moon up to a distance of 10 pc through the variation in flux (Limbach et al. 2024). Provided that polarimetric observations with the required precision of about 0.1% to 1% are also possible, the influence of an exomoon on the linear polarization of its exoplanet will potentially be detectable with the HWO.

On the other hand, if the stellar flux is included in the calculations, the degree of polarization caused by the planet and the moon together is less than 2 × 10⁻¹⁰ at all wavelengths of the considered model system. This is already clearly below the accuracy of a few parts per million reached by modern polarimeters (Bailey et al. 2020; Wiktorowicz et al. 2023). The influence of mutual events on the linear polarization is therefore on the order of 10⁻¹². Even the intrinsic polarization of radiation emitted by inactive stars of about 1 to 10 ppm (Kemp et al. 1987; Cotton et al. 2017) and its variability have to be considered very carefully to evaluate the possible impact of exomoons. While a particularly bright Hot Jupiter might cause a polarization of over 10 ppm (e.g., Bailey et al. 2021), it is unclear whether Hot Jupiters can even host exomoons, as they might lose their exomoons during their migration to an orbit close to the host star (Trani et al. 2020).

5.2 Usability of polarimetry for exomoon detections

To detect unknown exomoons through the observation of mutual events, the system has to be observed with a high duty cycle to prevent missing mutual events. For the model system, a mutual event only lasts between 1 and 5 hours. Therefore, such a system has to be observed with a time resolution of about 0.5 hours for the time of one lunar orbit period, that is about a month. Such observations are only possible from space (see Cabrera & Schneider 2007). In a system with higher inclinations, mutual events occur less often, which demands observations over an even longer time period.

In systems with smaller lunar orbital radii than the Earth- Moon system, mutual events occur more frequently and such systems are consequently easier to detect (see also Cabrera & Schneider 2007; Limbach et al. 2024). We note that while the time axis of our plots would have to be adapted accordingly, the changes △F/F and △P during a transit as a function of φ_ms are not affected by a smaller lunar orbital radius if the mutual illumination of the two bodies remains negligible except for combined transit and shadowing events at small phase angles. In the case of shadows, the values found for △F/F and △P would not be valid in the case of smaller lunar orbital radii.

Because the change in the polarization △P is highly dependent on the lunar orbit, the flux change △F should be measured simultaneously. To confirm that a measured change △F is caused by a moon, a prediction for future mutual events has to be made. The measurements of △P are then useful to narrow down the parameter space of the lunar orbital parameters. In Sect. 4.2.4, we showed that the path of the moon or its shadow over the planetary disk has a significant influence on the results and can be used to set limits on the inclination of the lunar orbit. After possible paths are extracted from measurements of △P, the radius of the moon and its orbital speed can be estimated from the duration of the mutual event. These constraints are sufficient to predict possible periods of time in which the next mutual event occurs. To confirm the existence of the exomoon, the system should then be observed during this period.

In some cases, it might also be possible to verify the presence of an exomoon with different approaches, for instance through transit timing or transit duration variations (e.g., Simon et al. 2007; Kipping 2009) or the orbital sampling effect (Heller 2014) in systems where the exoplanet transits the star or through observations at wavelengths where the moon outshines the planet (e.g., in the 1.4 μm water absorption band, see Limbach et al. 2024) using spectroastrometric observations (Agol et al. 2015).

This work does not take into consideration that the geometrical albedo of Earth-like exoplanets might vary significantly. An inhomogeneous surface, such as continents and oceans on Earth, can lead to a short-term variation of the polarized and total scattered radiation flux due to the planetary rotation (Song & Qu 2017; Sterzik et al. 2019; Wang et al. 2019; Groot et al. 2020). Such short periodic variations should not be mistaken for a moon. Variations in the planetary phase curves generally make it harder to detect an exomoon (see Limbach et al. 2024). Consequently, exoplanets with a less variable polarization phase curve offer a better chance to detect an exomoon. Planets covered in a dense cloud layer are therefore a reasonable expansion of our study.

6 Conclusion

We presented simulation results for the influence of an exomoon similar to the Earth’s moon on the flux and polarization phase curves of an Earth-like exoplanet at wavelengths of 400 nm, 550 nm, and 700 nm. We found that an inclined orbit of the exomoon causes considerable variations in the polarization signal observed during mutual events.

In Sect. 4.1, we first calculated the flux and linear polarization phase curves of the planet (F_p(α), P_p(α)) and the moon (F_m(α), P_m(α)) of the model presented in Sect. 3 and of the planet and moon combined (F(α), P(α)) using the threedimensional Monte Carlo radiative transfer code POLARIS (Reissl et al. 2016; Lietzow et al. 2021). For F(α) and P(α), we assumed that the planet-moon system is observed spatially resolved from the host star, while the planet and the moon remain spatially unresolved. The moon is assumed to be a diffuse reflecting sphere in the model (P_m = 0). Neglecting the influence of mutual events, the degree of linear polarization of the planet alone is therefore slightly higher compared to the two bodies combined, with a difference of less than 1% at all wavelengths, which is similar to the results by Berzosa Molina et al. (2018). However, we showed that these findings cannot be generalized. If the orientation of polarization of the disk-integrated planetary and lunar scattered radiation were identical and P_m were greater than P_p, the degree of linear polarization would be higher for the overall system of an exoplanet and its moon combined than for the planet alone. As pointed out by Berzosa Molina et al. (2018), an ambiguity remains even if the small differences could be measured. A slightly adapted planetary model can also explain a slightly higher flux and slightly changed degree of linear polarization. This ambiguity can be resolved through observations of mutual events. Our results for the influence of mutual events were presented in Sect. 4.2.

For occultations of the moon by the planet or lunar eclipses during which the moon is in the shadow of the planet, we calculated the changes in the total flux △F and linear polarization △P analytically in Sect. 4.2.1. For the model system, the biggest change in the linear polarization is an increase of △P = 0.77% for a lunar darkening event at α = 91° at a wavelength of 400 nm and +0.46% for a lunar darkening event at 550 nm. This is less than the values of around 2.5% which were previously reported by Berzosa Molina et al. (2018) for a similar system, although their flux values suggest a △P comparable to ours. Again, we found that for identical directions of polarization and P_m > P_p, the degree of linear polarization instead decreases during these mutual events. We discussed that with observations of these types of mutual events, it is possible to disentangle the linear polarization phase curves of the exoplanet and the exomoon.

The results from our polarized radiative transfer simulations of transits of the moon in front of the planet and time periods during which the moon casts a shadow on the planet were discussed in detail in Sects. 4.2.2–4.2.4. By using a Monte Carlo radiative transfer code, we were able to calculate the influence even for larger phase angles for which a completely spherical description of the planetary atmosphere is necessary. In Sect. 4.2.4, we placed the moon at a distance of r_p/2 to the planetary orbital plane to mimic the influence of an inclined exomoon orbit. This resulted in noncentral transits and shadows.

In general, the influence of a transit is higher than the influence of a shadow because of the relatively large penumbra compared to the umbra in the model system (Berzosa Molina et al. 2018). The highest relative flux changes △F/F are observed at large phase angles α_p when the moon is in the planetary orbital plane. At these phase angles, only a small crescent of the illuminated day side of the planet is visible to the observer. The moon occults a large portion of this crescent or darkens a large portion with its shadow. For transits, △F/F decreases at most by 32% to 36%. For shadows, the changes are smaller with decreases of 25% to 27%. For noncentral mutual events during which the moon has a distance of r_p/2 to the orbital plane of the planet, △F/F is smaller, reaching values of 25% for transits and 20% for shadows. From decreases in the flux alone it is hard to narrow down the space of possible exomoon properties because a smaller moon on an orbit with a smaller semimajor axis causes a very similar signal to the noncentral events.

We used radiative transfer simulations to also calculate the changes in the degree of linear polarization △P during a transit. For wavelengths of 400 nm and 550 nm the absolute value of △P is slightly higher than for 700 nm. Similar to Berzosa Molina et al. (2018), our simulations showed that increases are expected for central transits and shadows. However, we found that at large phase angles where the direction of polarization of the planet changes the degree of polarization decreases instead. The maximum increase in P is about 2.7% for wavelengths of 400 nm and 550 nm for central transits and about 1.8% for central shadows. For 700 nm, we find smaller increases of 1.9% for transits and 1.3% for shadows. Identical to the results by Berzosa Molina et al. (2018) for 450 nm, these maxima occur for transits or shadows observed at a phase angle of 96° to 106° just after quadrature, which offer the highest angular distance between the planet and the star. For noncentral transits and shadows, only a small increase of less than 0.5% is found at a wavelength of 400 nm in this phase angle range. Additionally, △P inhabits a wave-like sequence of increases and decreases that are similar in their absolute value instead of a dominant central peak. The largest increase in linear polarization during a noncentral mutual event is instead found at a phase angle of 166°, where the direction of polarization changes and is 0 without the mutual event. For noncentral transits △P = 1.5% at this phase angle, for noncentral shadows △P = 1.1%.

Because the maximum values and the shape of the linear polarization signal are dependent on the exact path that the moon or its shadow takes over the visible planetary disk seen by the observer, it can be used to narrow down possible values for the inclination of the lunar orbit. Therefore, measurements of △P during a mutual event facilitate the prediction of the timing of future mutual events, as discussed in Sect. 5.2.

If the planet-moon system were spatially unresolved from the host star, the stellar flux would dominate the total flux. Consequently, the resulting linear polarizations would be on the order of 10⁻¹². This is far below the precision of a few parts per million reached with modern polarimeters (Bailey et al. 2020; Wiktorowicz et al. 2023). However, enhanced polarimeters or coronagraphic observations might be able to measure the influence of mutual events of Earth-Moon analogs.

Acknowledgements

This research made use of NASA’s Astrophysics Data System (https://ui.adsabs.harvard.edu), Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018, 2022), Matplotlib (https://matplotlib.org/) (Hunter 2007) and Numpy (https://numpy.org/) (Harris et al. 2020). It was supported in part through highperformance computing resources available at the Kiel University Computing Center. M.M. and S.W. thank the DFG for financial support under grant WO 857/24-1. The authors thank the anonymous referee for constructive comments and suggestions.

Appendix A Moon-to-planet flux ratios

Fig. A.1 Ratios of the lunar flux to the planetary flux Fm/Fp computed according to Eq. (A.1). Fm/Fp depends on the ratio of the lunar to the planetary radius rm/rp and the ratio of the geometrical albedoes Am/Ap. Contour lines indicate equal flux ratios. Symbols indicate values for planet-moon pairs in the Solar System at a wavelength of 549 nm. The different planets are marked by different symbols; moons are indicated by their first letter(s). Circumplanetary rings are not included in the calculation. All values are from the NASA planetary fact sheets (D. R. Williams, https://nssdc.gsfc.nasa.gov/planetary/planetfact.html, as of May 7, 2024). Abbreviations of the moons: Earth: M (The Moon); Jupiter: I (Io), E (Europa), G (Ganymede), C (Callisto); Saturn: M (Mimas), E (Enceladus), TE (Tethys), D (Dione), R (Rhea), TI (Titan); Uranus: Ariel (A), Umbriel (U), Titania (T), Oberon (O), Miranda (M); Neptune: Triton (T), Nereid (N), Proteus (P); Pluto: Charon (C).

We assume that a planet and a moon are illuminated only by their central star and neglect the mutual illumination of the planet and the moon. The ratio of the flux reflected by a moon to the flux reflected by its planet is now approximated using their geometrical albedoes A_m, A_p and radii r_m, r_p, $$\frac{F_{\mathrm{m}}}{F_{\mathrm{p}}}\approx \frac{A_{\mathrm{m}}}{A_{\mathrm{p}}}\frac{r_{\mathrm{m}}^2}{r_{\mathrm{p}}^2}.$$(A.1)

Results from this equation for the moons in the Solar System at a wavelength of 550 nm are shown in Fig. A.1. The approximation is exact for a phase angle of 0°. Be aware that it does not include an additional influence of circumplanetary rings.

Appendix B Mutual illumination

While a single planet is only illuminated by its host star, a celestial body in a planet-moon system is also illuminated by the other body in the system. The overall flux of the system is made up of different photon paths which we group into different orders.

The zeroth-order flux is the flux of stellar radiation reaching the observer without interacting with the planet or the moon. This flux is (e.g., Zalewski 1995; Stam et al. 2004) $$F_{\mathrm{s}}=\frac{\pi r_{\mathrm{s}}^2}{d_{\mathrm{so}}^2}{B}_{\nu }(T),$$(B.1) where d_so is the distance between the star and the observer.

We now define the first-order flux to be the flux of stellar radiation reflected toward the observer by either the planet (F_1p) or the moon (F_1m) without interacting with the other celestial body. The first-order flux is therefore the total flux neglecting the mutual illumination if the planet and the moon are unresolved but the planet-moon system is resolved from the star. The fluxes F_1p and F_1m depend on the phase angle α_p of the planet and α_m of the moon. The flux phase curves of the two celestial bodies are thus $$F_{1\mathrm{p}}({\alpha }_{\mathrm{p}})=F_{\mathrm{p}0}{A}_{\mathrm{p}}{\psi }_{\mathrm{p}}({\alpha }_{\mathrm{p}}),F_{1\mathrm{m}}({\alpha }_{\mathrm{m}})=F_{\mathrm{m}0}{A}_{\mathrm{m}}{\psi }_{\mathrm{m}}({\alpha }_{\mathrm{m}}).$$(B.2)

Here, F_0p is defined in Eq. (4) and F_0m is defined similarly. A_p and A_m are the geometric albedoes. Also, we use the phase law ψ_p(α) = F_1p(α)/F_1p(0) of the planet and the phase law ψ_m of the moon that is defined analogously (see Lester et al. 1979).

Flux of the n-th order is scattered by either the planet or the moon n times. Therefore, second-order flux is the flux that was reflected by the planet or the moon toward the other celestial body once before being scattered toward the observer. If r_m ≪ d_pm, the flux of the stellar radiation reflected by the moon measured on the planetary surface is $$F_{1\mathrm{m},\mathrm{p}}=\pi B_{\nu }(T)\frac{r_{\mathrm{s}}^2{r}_{\mathrm{m}}^2}{d_{\mathrm{sm}}^2{d}_{\mathrm{pm}}^2}{A}_{\mathrm{m}}{\psi }_{\mathrm{m}}(\pi -{\phi }_{\mathrm{ms}}),$$(B.3) where φ_ms is defined in Fig. 1. The second-order flux of the planet F_2p is therefore calculated by replacing the stellar flux on the planetary surface $\pi p{B}_{\nu }(T)r_{\text{s}}^2/d_{\text{sp}}^2$ in Eq. (B.2) with the flux F_1m,p. Consequently, using d_po ≈ d_mo $$\begin{array}{l}{F}_{2\mathrm{p}}({\phi }_{\mathrm{ms}},{\phi }_{\mathrm{mo}})=\pi B_{\nu }(T)\frac{r_{\mathrm{s}}^2{r}_{\mathrm{p}}^2{r}_{\mathrm{m}}^2}{d_{\mathrm{sm}}^2{d}_{\mathrm{pm}}^2{d}_{\mathrm{po}}^2}{A}_{\mathrm{p}}{A}_{\mathrm{m}}{\psi }_{\mathrm{m}}(\pi -{\phi }_{\mathrm{ms}}){\psi }_{\mathrm{p}}({\phi }_{\mathrm{mo}})\\ =F_{1\mathrm{m},\max }\frac{r_{\mathrm{p}}^2}{d_{\mathrm{pm}}^2}{A}_{\mathrm{p}}{\psi }_{\mathrm{m}}(\pi -{\phi }_{\mathrm{ms}}){\psi }_{\mathrm{p}}({\phi }_{\mathrm{mo}}).\end{array}$$(B.4) where φ_mo is the angle between the direction toward the moon and the observer at the center of the planet. For a non-inclined system as in Fig. 1, we have φ_mo = φ_ms + α_p. F_1p,max and F_1m,max are the maximum values of the first-order fluxes. For the second-order flux of the moon F_2m, we get $$F_{2\mathrm{m}}({\phi }_{\mathrm{ms}},{\phi }_{\mathrm{mo}})=F_{1\mathrm{m},\max }\frac{r_{\mathrm{p}}^2}{d_{\mathrm{pm}}^2}{A}_{\mathrm{p}}{\psi }_{\mathrm{m}}({\phi }_{\mathrm{ms}}){\psi }_{\mathrm{p}}(\pi -{\phi }_{\mathrm{mo}}).$$(B.5) as the distances of the planet or the moon to the central star are comparable, d_sp ≈ d_sm. Now, we set a boundary for the secondorder flux maximum in comparison to the first-order influence of the moon, $$\frac{F_{2\mathrm{p},\max }+F_{2\mathrm{m},\max }}{F_{1\mathrm{m},\max }}<2\frac{r_{\mathrm{p}}^2}{d_{\mathrm{pm}}^2}{A}_{\mathrm{p}}.$$(13)

For the model system described in Sect. 3, the right-hand side is less than 1.9 · 10⁻⁴. In that case, the second-order flux and thus the mutual illumination are negligible. However, this cannot be generalized to all planet-moon systems. If the moon orbits its planet near the Roche limit at about 2.5r_p (e.g., Williams 2003), the value of the right-hand side rises to about 0.1. Whether the mutual illumination can be neglected when calculating the influence of an exomoon on the phase curve of an exoplanet has to be decided for each system individually.

Appendix C Half width of shadows and eclipses

For a transit, the half width is calculated from Fig. C.1, $$\mathrm{\Delta }{\phi }_{\mathrm{ms}}=\arcsin \frac{r_{\mathrm{p}}+r_{\mathrm{m}}}{d_{\mathrm{pm}}}.$$(C.1)

Here, we assumed that the distance between the observer and the planet is much larger than the distance between the planet and the moon. This is the case for Solar System-based observations of exoplanets and exomoons.

The half width of a shadow phase is calculated according to Fig. C.2. In this case, three contributions were identified. The quantity ϕ_s is the angle under which the star is seen from the planetary center. It is interpreted as describing the time the center of the moon needs to cross the stellar disk from the perspective of an observer on the planet. However, an extended moon needs slightly longer to cross the stellar disk. This effect is included by adding the angle under which the moon is seen from the planetary center, ϕ_m. Finally, ϕ_p is added to accommodate the time the shadow of the moon needs to cross the planet’s surface. From Fig. C.2, we get $$\mathrm{\Delta }{\phi }_{\mathrm{ms}}<\arcsin \frac{r_{\mathrm{s}}}{d_{\mathrm{sp}}}+\arcsin \frac{r_{\mathrm{p}}}{d_{\mathrm{pm}}}+\arcsin \frac{r_{\mathrm{m}}}{d_{\mathrm{pm}}},$$(C.2) where we used an upper limit for ϕ_p that is reached when d_sp approaches infinity (for a very high distance between the planet and the star).

Fig. C.1 Half width of a transit △φms calculated from the triangle between the direction toward the observer, the center of the planet P and the positions M1 and M2 of the center of the moon at the beginning and the end of the transit. All arrows point toward the observer and are roughly parallel. The figure is not to scale.

Fig. C.2 Half width of a shadow phase △φms estimated from the three angles ϕm, ϕp, and ϕs.

Appendix D Influence of the lunar polarization phase curve during transits

Here, we estimate the influence of the actual lunar polarization phase curve on the changes △P in the linear polarization that are observed during transits. Again, F_p, P_p are the flux and degree of linear polarization of the planet, respectively, and F_m, P_m the same quantities for the moon. We assume that the polarized flux of the planet is reduced from F_pP_p to aF_pP_p and the total flux is reduced from F_p to bF_p, where a, b ∈ [0, 1]. Furthermore, we assume that the direction of polarization is identical for the planet and moon. Assuming parameters similar to the moon of Earth, this is the case at all phase angles between the inversion phase angle of the moon of about 30° (e.g., Lyot 1929; Coyne & Pellicori 1970; Shkuratov et al. 1992) and the change in the direction of polarization at about 166° for a wavelength of 400 nm. Also, the direction of polarization is assumed to be constant across the planetary disk which is approximately true for the above-mentioned range of phase angles (see Fig. 5). This results in a change in the degree of linear polarization of $$\mathrm{\Delta }P=\frac{a{F}_{\mathrm{p}}{P}_{\mathrm{p}}+F_{\mathrm{m}}{P}_{\mathrm{m}}}{b{F}_{\mathrm{p}}+F_{\mathrm{m}}}-\frac{F_{\mathrm{p}}{P}_{\mathrm{p}}+F_{\mathrm{m}}{P}_{\mathrm{m}}}{F_{\mathrm{p}}+F_{\mathrm{m}}},$$(D.1) thus, $$\begin{array}{l}\mathrm{\Delta }P=(\frac{a}{b{F}_{\text{p}}+F_{\text{m}}}-\frac{1}{F_{\text{p}}+F_{\text{m}}})F_{\text{p}}{P}_{\text{p}}\hfill \\ +(\frac{1}{b{F}_{\text{p}}+F_{\text{m}}}-\frac{1}{F_{\text{p}}+F_{\text{m}}})F_{\text{m}}{P}_{\text{m}}.\hfill \end{array}$$(D.2)

If we let P_m = 0, the systematic error in △P is given by the second term of Eq. (D.2). If F_p/F_m ≫ 1/b and F_p/F_m ≫ 1, the term simplifies to $$\mathrm{\Delta }{P}_{\mathrm{sys}}\approx \left(\frac{1}{b}-1\right)\frac{F_{\mathrm{m}}}{F_{\mathrm{p}}}{P}_{\mathrm{m}}.$$(D.3)

In our simulations, F_p/F_m > 43 at 700 nm and even higher at wavelengths of 400 nm and 550 nm. With 1/b < 1.5, the assumption is justified. We note that the simplification might not be possible if the parameters of the planet and the moon are different, especially when the moon has a similar or higher surface albedo than the planet. In that case, the second term of Eq. (D.2) should be used instead of Eq. (D.3).

△P_sys is always positive for equal polarization directions of the planet and the moon. In our study, using a realistic maximum value of P_m = 10% (see Lyot 1929; Coyne & Pellicori 1970) leads to △P_sys ≤ 0.08% for a wavelength of 400 nm. While this may result in a slight underestimation of positive changes △P and a slight overestimation of negative changes △P, it does not change the general behavior shown in Fig. 4. Even the slight decreases in the polarization direction at the end and beginning of a transit are unaffected because 1 − b is small when the moon is only partially in front of the edge of the planet or the terminator, so 1/b is close to 1.

Appendix E Disk resolved polarized flux

In Fig. E.1 and E.2, we present results for the disk-resolved polarized flux of the planet in our model system (see Sect. 3) at wavelengths of 550 nm and 700 nm, respectively.

Fig. E.1 Same as Fig. 5, but for a wavelength of 550 nm.

Fig. E.2 Same as Fig. 5, but for a wavelength of 700 nm.

References

All Tables

All Figures

	Fig. 1 Geometry of the non-inclined baseline system. The centers of the star, the planet, and the moon as well as the observer are all coplanar. The angles αp and αm are the phase angles of the planet and the moon, respectively. φms is the angle between the direction to the star and the direction to the moon for an observer at the center of the planet.
In the text

Fig. 2 Phase curves of the planet, the moon, and the planet and moon together but without the influence of mutual events. Top: Flux F/F0 as a function of the phase angle α. For the planet and for planet and moon combined, F0 = F0p. In contrast, for the moon, F0 = F0m. In addition, we use the analytic result from Russell (1916) for a Lambertian sphere for the phase curve of the moon. Middle: Phase curves for the linear polarization P(α) for the planet and planet and moon combined, assuming Pm = 0 (i.e., of the diffuse reflecting moon surface). The stellar flux is not included in the calculation of P. Bottom: Difference P − Pp between the polarization phase curve of the unresolved system consisting of planet and moon and the planetary polarization phase curve.

In the text

Fig. 3 Four types of mutual events simulated for a phase angle of 225°. The RGB images are based on individual maps simulated for wavelengths of 440 nm (blue), 555 nm (green), and 700 nm (red). The total flux was normalized to the geometric albedo separately at each of the wavelengths. The dashed circles indicate the position and size of the planet (blue) and moon (gray). To improve the visibility of the moon and to fit both bodies into the images for shadow and eclipse, the lunar surface albedo am was increased to 0.3 and the moon’s semimajor axis was set to 2.8 rp for this graphic.

In the text

Fig. 4 Relative change of the total flux △F/F (top row) and absolute change of the degree of polarization △P (bottom row) during central transits of the moon in front of the planet at different phase angles αp visualized using color maps and contour lines. The upper abscissa shows the time in days since the planet was at α = 0°. On the left ordinate, the relevant angle αp + φms is shown. The angle was transformed to a time in hours on the right ordinate, where a time of 0 corresponds to the time when the planet, moon, and observer are perfectly aligned. A single transit of the moon corresponds to a passage through the color maps perpendicular to the abscissa. For phase angles above 180°, F(360° − αp, 360° − φms) = F(αp, φms).

In the text

Fig. 5 Polarization maps for different phase angles at a wavelength of 400 nm. The polarized flux FpPp is encoded as a color map. The resolution is 60 × 60 px, the planet has a diameter of 50 px. Polarization directions are shown as white lines. The larger orange circle marks the planet’s circumference. The orange dashed lines between the poles of the planet mark the terminator. The day side of the planet is to the right of it. The horizontal orange lines mark the path of the moon during a central transit (lower pair of lines), or a noncentral transit (upper pair of lines). For two positions on the edge of the planet, the circumference of the moon is shown as an orange circle.

In the text

Fig. 6 Similar to Fig. 4, but for phases during which the moon casts a shadow on the planet. The left ordinate now shows the angle φms, which specifies the lunar position. Again, this angle was transformed to a time in hours, as shown on the right ordinate, where a time of 0 corresponds to the time when the planet, moon, and star are perfectly aligned. A single run of the shadow of the moon across the planetary surface corresponds to a passage through the color maps perpendicular to the abscissa.

In the text

	Fig. 7 Same as the first column of Fig. 4, but for the results of our simulations with noncentral transits.
In the text

	Fig. 8 Same as the first column of Fig. 6, but for the results of our simulations with noncentral shadow phases.
In the text

Fig. A.1 Ratios of the lunar flux to the planetary flux Fm/Fp computed according to Eq. (A.1). Fm/Fp depends on the ratio of the lunar to the planetary radius rm/rp and the ratio of the geometrical albedoes Am/Ap. Contour lines indicate equal flux ratios. Symbols indicate values for planet-moon pairs in the Solar System at a wavelength of 549 nm. The different planets are marked by different symbols; moons are indicated by their first letter(s). Circumplanetary rings are not included in the calculation. All values are from the NASA planetary fact sheets (D. R. Williams, https://nssdc.gsfc.nasa.gov/planetary/planetfact.html, as of May 7, 2024). Abbreviations of the moons: Earth: M (The Moon); Jupiter: I (Io), E (Europa), G (Ganymede), C (Callisto); Saturn: M (Mimas), E (Enceladus), TE (Tethys), D (Dione), R (Rhea), TI (Titan); Uranus: Ariel (A), Umbriel (U), Titania (T), Oberon (O), Miranda (M); Neptune: Triton (T), Nereid (N), Proteus (P); Pluto: Charon (C).

In the text

	Fig. C.1 Half width of a transit △φms calculated from the triangle between the direction toward the observer, the center of the planet P and the positions M1 and M2 of the center of the moon at the beginning and the end of the transit. All arrows point toward the observer and are roughly parallel. The figure is not to scale.
In the text

	Fig. C.2 Half width of a shadow phase △φms estimated from the three angles ϕm, ϕp, and ϕs.
In the text

	Fig. E.1 Same as Fig. 5, but for a wavelength of 550 nm.
In the text

	Fig. E.2 Same as Fig. 5, but for a wavelength of 700 nm.
In the text