© ESO 2021

1 Introduction

As a solvent, a liquid phase seems more favorable for biochemical reactions than gas and solid phases (Benner et al. 2004). Therefore, the discovery of planets with a liquid ocean has crucial importance in the context of the search for extraterrestrial life. How can we find an exoplanet with a surface ocean?

With regard to non-polarimetric signatures of oceans, two types signature have been confirmed by the astronomical observations of Earth. These are (i) photometric variations caused by the rotation of a planet covered by a surface with inhomogeneous reflectances and/or colors, and (ii) the spectroscopic signature of glint (specular reflection) from an ocean.

For the first type, model calculations for an Earth-like planet expect a diurnal intensity (and color) variation because of the rotating brighter (redder) continents and darker (bluer) oceans (Ford et al. 2001; Oakley & Cash 2009; Fujii et al. 2010); however, it was also expected that the existence of clouds would make surface-type determination much more difficult than the cloud-free case (Ford et al. 2001; Oakley & Cash 2009). From an examination of the multiband light curves of Earth obtained from space by the EPOXI (Extrasolar Planet Observation and Characterization (EPOCh) + Deep Impact Extended Investigation (DIXI)) mission, a correlation between color and ocean fraction was identified (Cowan et al. 2009).

The second type of signature was identified from Earth’s spectra. Robinson et al. (2014) observed Earth from the Lunar CRater Observation and Sensing Satellite (LCROSS) and showed good agreement with the model spectrum that considered glint. Enhanced intensity contrasts between wavelengths sensitive to the surface (e.g., ~1.6 μm) and those insensitive due to strong molecular absorption (e.g., ~1.4 μm) are key features indicating the contribution of glint from the oceans.

A polarimetric signature of oceans is also expected, and thus, polarimetry has the potential to be another powerful technique. Specular reflection (glint) from a smooth liquid surface is highly polarized (~100% when the incident angle equals Brewster’s angle), as expressed by the Fresnel equations. Some researchers theoretically studied the polarization of light reflected from planets with an ocean (McCullough 2006; Stam 2008; Williams & Gaidos 2008; Zugger et al. 2010, 2011; Kopparla et al. 2018). A common conclusion was that a cloud-free planet with a full ocean coverage exhibits a very high peak polarization degree (>70%), which is significantly larger than that of planets with other surface types. However, it is commonly noted that diffuse scattering by clouds, atmospheric Rayleigh scattering, and various other effects dilute the polarimetric signature of the oceans. Most of these theoretical works calculated the intensity and polarization of the glint considering wind-driven tilts of the ocean surface based on the Cox–Munk model (Cox & Munk 1954).

One benefit of polarimetry is that the degree of polarization is virtually insensitive to telluric extinction because of the nature of the relative measure (Stam 2008) that makes polarimetry applicable both in space and on the ground. In contrast, the photometric and spectroscopic signatures seem to be severely affected by telluric extinction and its variability, which make them more suitable for space observations than ground-based observations.

Another benefit is the compatibility with high-contrast observations. Polarimetric differential imaging (PDI) is a widely utilized technique to enhance the contrast performance to detect polarized objects (planets and disks) around an unpolarized central star.The PDI technique can be used not only to detect planets but also to measure their polarization (Murakami et al. 2006; Takahashi et al. 2017). In addition, the fact that multiwavelength observation is not required to detect the polarimetric signature is favorable for high-contrast observations because high-contrast optics are very sensitive to wavelength, and therefore the effective bandwidth is currently limited to 10% of the central wavelength (N’Diaye et al. 2016; Llop-Sayson et al. 2020). This seems advantageous to the polarimetric technique, in comparison with the techniques based on the photometric (color) variation and the spectroscopic contrasts.

We review previous polarimetry of Earth. The detectability of the polarimetric signature remains uncertain because of the lack of a conclusive detection of the polarimetric signature of the oceans in the disk-integrated Earth light. The POLarization and Directionality of the Earth’s Reflectances (POLDER, Deschamps et al. 1994) instruments on Earth-orbiting satellites conducted polarimetry of Earth. They measured the polarization of reflected light from various types of media, including the ocean, from different viewing angles, while sweeping Earth with a swath of ~ 2000 km. However, when searching an exoplanetary ocean, we need a polarization measurement of the disk-integrated Earth light as the benchmark because it will be extremely difficult to spatially resolve an exoplanet in the foreseeable future. In principle, it is possible to calculate the net Earth polarization by summing up polarized intensities from different media on Earth. Wolstencroft & Breon (2005) performed such a calculation using the POLDER measurements for six types of media, including cloud-free ocean and cloudy ocean. However, they had to assume a certain global cloudfraction and uniformity of cloud distribution. We still require the observations of the Earth’s whole disk at onceto confirm the significance of the polarimetric signature of the oceans in the disk-integrated Earth light.

In this sense, it would be ideal to conduct the astronomical polarimetry of Earth from a space probe. To our knowledge, Voyager 1 and the Pioneer Venus Orbiter measured the polarization degree of the disk-integrated Earth for large phase angles (Coffeen 1979). However, with a very limited number of the measurements, it is impossible to identify the ocean signature in the net Earth polarization.

Ground-based observations of Earthshine on the Moon can be considered an alternative technique. Earthshine is Earth’s reflected light, which is back-scattered from the Moon to Earth; it appears as a faint glow on the night side of the Moon. It is commonly utilized to observe Earth as an exoplanet (a remote disk-integrated planet), as reviewed by Arnold (2008) and Pallé (2010). If the glint from the oceans provides a significant contribution to the polarization of the total reflected light on Earth, a higher Earthshine polarization degree is expected when we have a glint spot in the view from the Moon than that when we do not have it. Because there is a higher probability that a glint is seen from the Moon when we have a larger ocean fraction in the Earthshine-contributing region (the region illuminated by sunlight and viewable from the Moon), the dependence of the Earthshine polarization degree on the ocean fraction may be observed.

In the past, we measured Earthshine polarization degrees in visible wavelengths on 19 nights. Although it was implied that Earthshine from an ocean-dominant Earth surface has a higher polarization degree than that from a land-dominant surface, the difference was not statistically significant because of the large observational errors (Takahashi et al. 2012).

Sterzik et al. (2012, 2019) presented polarization degree spectra (in the visible wavelengths) of Earthshine observed with the Very Large Telescope in Chile. They showed that the polarization degrees of Earth with the Pacific in view were significantly higher than those with the Atlantic in view. As possible causes of the polarization difference, Sterzik et al. (2019) suggested two factors: (a) different cloud coverages, and (b) a larger contribution of ocean glint from the Pacific side; however, conclusive evidence regarding the ocean signature has not been presented yet.

Although most of the previous Earthshine polarimetry was conducted in the visible wavelengths (<1 μm), Zugger et al. (2011) pointed out that near-infrared polarimetry is more favorable to the search for an ocean than that at the visible wavelengths because atmospheric Rayleigh scattering is reduced in the near-infrared wavelengths. Our analysis on the phase variation of the Earthshine polarization degree spectra (Takahashi et al. 2013) suggested that Earth’s polarization degrees in the near-infrared should be more sensitive to surface properties than that in the visible wavelengths. This suggestion is consistent with the observations by Miles-Páez et al. (2014), who presented a visible-to-near-infrared (0.4–2.3 μm) polarization degree spectrum of Earthshine: the visible polarization degrees rapidly decreased with increasing wavelengths, which can be explained by atmospheric Rayleigh scattering; the near-infrared spectrum was virtually flat except at atmospheric molecular bands (H₂O and O₂), which implies that Rayleigh scattering is almost ineffective for polarization degrees in the near-infrared. With a single-night observation, it is difficult to identify a contribution from the sea glint to the polarization. Based on the previous Earthshine polarimetric results and discussions, we launched a project to perform near-infrared polarimetry for Earthshine to detect the possible dependence of the polarization degree on the ocean fraction.

This paper is organized as follows. In Sect. 2, our observations and data reduction methods are summarized. The main results are presented in Sect. 3, where we retrieve the polarimetric signature of the oceans fromthe observed data. After the examinations of possible impacts by factors other than the oceans in Sect. 4, we discuss the distinctiveness of the polarimetric signature of an Earth-like ocean and feasibility offuture polarimetric search for exoplanetary oceans in Sect. 5. We conclude this article in Sect. 6. The details of the data reduction methods and the observation-model comparison are provided in the appendices.

2 Observations and data reduction

We conducted near-infrared polarimetry for Earthshine using the Nishiharima Infrared Camera (NIC, Ishiguro et al. 2011) mounted at the Cassegrain focus (f∕12) of the 2.0 m Nayuta altazimuth telescope at the Nishi-Harima Astronomical Observatory(134.3356°E, 35.0253°N, and 449 m in altitude). The NIC is equipped with three detector arrays and two dichroic mirrors that enable simultaneous J- (central wavelength: 1.25 μm), H- (1.63 μm), and K_s-band (2.15 μm) imaging. In the imaging polarimetry mode, a rotatable half-wave plate and a polarizing beam displacer are inserted in the optical path(Takahashi et al. 2018; Takahashi 2019). A pair of ordinary and extraordinary images with a size of ~ 24″ × 69″ is obtained with a single exposure.

Observations were conducted between May 2019 and April 2020 (Table A.1). Valid data were obtained for 32 nights. The Moon was in the waxing phase for 20 nights and in the waning phase for the other 12 nights. As observed from Japan, Earthshine on the waxing Moon is usually contributed by the Eurasian and African continents and the Indian Ocean, whereas that on the waning Moon originates from the Pacific Ocean and the Americas (Fig. 1). Under our observation conditions, the ocean fraction (with consideration of the cloud distribution) in the Earthshine-contributing region ranged from ~ 15–40% for the waxing phase, and ~20–45% for the waning phase. On average, the ocean fraction is larger in the waning phase than in the waxing phase. We covered a wide range of ocean fractions (~15–45%), which allowed us to investigate the possible dependence of the Earthshine polarization degrees on the ocean fraction. The ocean and land fractions also vary on an hourly timescale because of the Earth’s rotation, as shown in Figs. 1c, d, and this enabled us to explore the hourly variations of the Earthshine polarization degrees.

To minimize the undesired effects caused by observing different lunar locations, we conducted observations according to the following procedure. On each observing night, we first pointed the Nayuta telescope toward the crater Grimaldi (selenographic coordinate: 68.6°W, 5.2°S) in the waxing phase and the crater Neper (84.5°E, 8.8°N, east of Mare Crisium) in the waning phase, after correcting the pointing error measured using a nearby star. Both craters are near the lunar edge (distances ≲2′). Then, we scanned the Moon along the RA axis until the edge of the Moon was placed near the center of the field of view (FOV). An example of the observed (and reduced) images is shown in Fig. 2. Our target locations are not on a major maria and near sites repeatedly observed in previous Earthshine photometry because they were expected to have roughly comparable albedos (Qiu et al. 2003; Pallé et al. 2004; Montañés-Rodríguez et al. 2007). Half of the FOV was reserved for the sky, which allows the sky background intensities and their positional gradients to be measured. The position angle of the instrument (ϕ_inspa) was maintained at 90° from the equatorial north, as measured counter-clockwise, so that the long side of the FOV was aligned with the RA axis. Telescope tracking was conducted in accordance with the sky motion of the Moon, which was calculated at the Jet Propulsion Laboratory (JPL) Horizons system¹. Because the tracking was not perfect, we shifted the telescope east or west with a typical interval of ~ 30 min so that the lunar edge remained near the center of the FOV. Features on the Moon were hardly recognizable in the raw images because of the dim Earthshine and strong scattered light from the day side of the Moon, though we were able to visually identify the lunar edge in most cases². Despite our efforts, the actually observed location may have varied night by night even within one phase (waxing phase or waning phase), or on an hourly timescale during a single night. Possible impacts induced by different lunar locations (namely different degrees of depolarization) are discussed in Sect. 4.1.

The exposure time for a single frame was usually 20–180 s depending on the brightnesses of the Earthshine and the sky.A series of four exposures corresponding to four different rotation angles of the half-wave plate (ϕ_hwp = 0°, 45°, 22.5°, and  67.5°) produced a set of normalized Stokes parameters q = Q∕I and u = U∕I. We call this single series a “sequence” of observations. With a typical interval of ~ 30 minutes, we observed a blank sky region 60″–90″ east or west of the observing lunar edge. The exposure time for the blank sky observations was set to be the same as that for the Earthshine observations.

After basic image processing including flat fielding and the subtraction of the sky background, the maps of normalized Stokes parameters q = Q∕I and u = U∕I were produced (Fig. 2). The values of q and u in a region of ~16″ × 8″ are extracted and averaged. The polarization degree (fractional polarization, P) and polarization position angle (Θ) are converted from q and u with a correction of positive bias (Plaszczynski et al. 2014). We confirm that the derived Θ is almost always perpendicular to the scattering plane (the plane that includes the Sun, Earth, and Moon), as shown in Fig. 3, and this supports the fact that we successfully extracted the polarization by the reflection of sunlight by Earth. Details on the data reduction are presented in Appendix B.

Fig. 1 Views of cloud-free Earth from the Moon at different observation times. At each panel, the illuminated hemisphere is the Earthshine contributing region. These images were created with the Earth and Moon Viewer3 developed by John Walker.

3 Results

3.1 Nightly means

All observed polarization degrees (P) as nightly means are summarized in Table A.2 and illustrated in Fig. 4. The only previous near-infrared polarimetry for Earthshine (Miles-Páez et al. 2014) reported a P of ~ 3–5% at α ~ 100°, which approximately agrees with our results. The observed P increased with the increasing Sun-Earth-Moon phase angle (α), and it reached its peak of ~4% or larger at an α between 120° and 150°. The overall shape of the P phase curve and α for the peak P agree with theoretical predictions by Williams & Gaidos (2008), Zugger et al. (2011), and Kopparla et al. (2018), who calculated the polarization degree of an ocean planet in the near-infrared wavelengths (or considering no contribution from atmosphericRayleigh scattering).

Our primary focus is the possible dependence of the polarization degree on the ocean fraction. The ocean fraction in the Earthshine-contributing region is expressed by the darkness of the plot colors in Fig. 4. This set of figures provides an interesting impression that data points with a larger ocean fraction tend to have a larger P than those with a smaller ocean fraction at a similar α.

We performed the following analysis to illustrate the possible dependence of P on the ocean fraction in a more quantitative manner. In general, the polarization degree of reflected light depends on both properties of the reflecting body and phase angle. We fit a curved line to all data points (except some outliers) in Fig. 4 (see Appendix B.4 for details of the fitting). The fit curve, denoted by P_mean, corresponds to the polarization phase curve (phase angle dependence of polarization degrees) for the typical scene combination on the Earthshine contributing region. The contrast of the observed P to P_mean for the same α represents the extent to which P deviates from the polarization degree of the typical Earth scene, and it clarifies the discussion on the dependence of P on the actual Earth scene because the phase-angle dependence is suppressed.

Figure 5 (top row) displays P∕P_mean plotted against the ocean fraction. We see a clear positive correlation of P∕P_mean with the ocean fraction for all J, H, and K_s bands. In other words, P tends to be larger when we have a larger ocean fraction if α is fixed. This is probably attributed to the greater contribution from highly polarized sea glint. We also plot P∕P_mean against land fraction and cloud fraction (Fig. 5, middle and bottom rows). In Fig. 5 (middle row), P∕P_mean appears to be negatively correlated with the land fraction, though the correlation is less clear than that for the ocean fraction. In Fig. 5 (bottom row), no clear correlation of P∕P_mean is found with the cloud fraction.

We classified clouds into three types based on the cloud top height (h_top). Following Lamb & Verlinde (2011), we defined clouds for h_top ≥ 7 km as high clouds, those for 2 km ≤ h_top < 7 km as middle clouds, and those for h_top < 2 km as low clouds. Figure 6 explores the possible dependences of P∕P_mean on fractions of high clouds (top row), middle clouds (middle row), and low clouds (bottom row). Although it is interesting that high clouds and the other types appear to have opposite dependences, none of the correlations are as strong as that for the ocean fraction (Fig. 5, top row).

The stronger correlation with the ocean fraction (f_o) than that with the land fraction (f_l) or cloud fraction (f_c) can be explained as follows. The three types of scenes can be divided into two groups: (i) oceans as a strong polarizer (because of specular reflection), and (ii) lands and clouds as weak polarizers (because of multiple scattering). Because polarimetric effects from the lands and the clouds are (very) roughly similar, the net polarizationof Earth should be largely determined by the ratio of f_o to (f_l +f_c), regardless of the specific values of f_l and f_c. Hence, the strong correlation of P∕P_mean with f_o was observed. Once f_o is given, (f_l + f_c) is automatically fixed since we have a relation of f_o + f_l + f_c = 1. In contrast, even if f_c (or f_l) is given, the ratio of f_o to f_l (f_c) should have a significant impact on the net Earth polarization. This is probably the reason why we observed a weaker correlation with f_c (f_l) than that with f_o. The weaker correlation with f_c does not deny a major role of clouds in the net Earth polarization. The type of surface covered by clouds is important.

Fig. 2 Intensity (left), Stokes q0 (middle), and u0 (right) images observed on 2020 January 3. North is right, and east is up. The FOV is ~ 19″ × 64″ (smaller than original FOV because of trimming). The parallelograms are the sampling regions. The values of q0 and u0 on the sky are scattered because of division of ~0 by ~ 0.

Fig. 3 Observed position angles of Earthshine polarization as plotted against position angles normal to scattering plane.

Fig. 4 Earthshine polarization degrees (P) in J (left), H (middle), and Ks (right) bands, plotted against Sun–Earth–Moon phase angle (α). The ocean fraction was calculated with concentrated weighting (see Appendix B.5 for details on the derivation of the fractions). The dashed lines represent the polarization phase curve for the typical Earth scene. They are derived from fitting a curved line to all data points (with some exceptions described below) with free parameters w (single scattering albedo) and s (scaling factor). The crosses represent data points corresponding to |Θ − N| > 15° (where Θ denotes the position angle of polarization and N denotes the position angle normal to the scattering plane) or α < 50°, which were excluded from the fitting (see Appendix B.4 for details on the fitting). These figures give an impression that P for a larger ocean fraction (plots with a darker color) tends to be larger than those for a smaller ocean fraction at a similar α, which suggests that the contribution from the sea glint (specular reflection) enhances the polarization degree of Earth.

Fig. 5 Dependence of polarization degree on ocean (top), land (middle), and cloud (bottom) fractions (in J, H, and Ks bands from left to right). Each polarization degree (P) in Fig. 4 is divided by typical polarization (Pmean: dashed line in Fig. 4) at the phase angle, and then plotted against the ocean, land, or cloud fraction (fo, fl, or fc, respectively). The filled and open plots correspond to observations in the waxing and waning phases, respectively. The dashed lines are regression lines of the form af + b with a correlation coefficient r. The crosses correspond to those in Fig. 4 and were excluded from the linear regression. The fractions was calculated with concentrated weighting.

Fig. 6 Dependence of polarization degree on high- (top), middle- (middle), and low-cloud (bottom) fractions (in J, H, and Ks bands from left to right). The legends are the same as those in Fig. 5.

3.2 Hourly variations

Fractions of scene types on the Earthshine contributing region vary on an hourly timescale corresponding to Earth’s rotation (see Figs. 1c, d for 2020 January 3). Hence, it is possible to observe the hourly variation of P in accordance with the scene transition. We investigated the time variation of P on six dates, on which we made a valid observation for more than two hours on a single night. Time-resolved P values from all six dates are divided by P_mean, as obtained from Fig. 4, and plotted against the ocean, land, and cloud fractions in Fig. 7. Similarly to what is seen in Fig. 5, a clear positive correlation of P∕P_mean with the ocean fraction is deduced again.

Time-series P on the six dates is presented in Figs. 8–9, with scene fractions and the observed position angle of polarization. Among the six dates, we observed significant variations in P on three dates (2019 December 18, 2020 January 3, and 2020 March 2) in all three bands (Fig. 8, left). The ratios of peak-to-peak variation (ΔP) to the averaged polarization degree ($\overline{P}$) range from ~0.2 to ~1.4 (Table 1). The position angle of the polarization was almost constant and it was confined to be perpendicular to the scattering plane for all six dates and all three bands (Figs. 8–9, right).

We attempted to reproduce the observed hourly variations in P (including the non-variations) based on scene fractions at the time, by referring to a model of planetary reflected light (Williams & Gaidos 2008). The detailed description of the model is provided in Appendix C. We see an excellent (2020 January 3 and 2020 March 2) or fairly good (2019 November 21, 2019 December 19, and 2020 April 29) agreement between the observed P and modeled P, except for on 2019 December 18 (Figs. 8–9). The disagreement on 2019 December 18 may be attributed to insufficient time resolution of the referred data of cloud distribution (see Appendix C). We discuss other possible causes (such as lunar depolarization, telluric polarization, and artificial polarization) in Sect. 4. We observe a resemblance between the time-variation of the ocean fraction and that of the modeled P (Fig. 8); this indicates that mainly the ocean fraction controls the hourly variation in the polarization degree of Earth.

Fig. 7 Same figure as Fig. 5, except the data source is the time-resolved data on 2019 November 21, 2019 December 18, 2019 December 19, 2020 January 3, 2020 March 2, and 2020 April 29. The data on the same date are connected by lines.

4 Impacts of misleading factors

4.1 Depolarization at the lunar surface

The polarization of Earthshine is not the same as the polarization of Earth as observed from outside the planet. The light from Earth is depolarized when it is back-scattered from the lunar surface (Dollfus 1957). The depolarizing factor (or polarization efficiency, ϵ) of the Moon at the near-infrared wavelengths is not well known. It is known that ϵ depends on surface albedo (Dollfus 1957; Bazzon et al. 2013), and a medium with a higher albedo has a lower ϵ. Bazzon et al. (2013) derived an empirical formula (Eq. (9) of that paper) of ϵ as a function of albedo and the wavelength, which is valid in the visible wavelengths. When we extend the formula to near-infrared wavelengths (1.2–2.2 μm) with typical highland albedos (0.15–0.25 in visible wavelengths), ϵ ~ 0.2–0.3 is deduced. Hence, the observed Earthshine polarization degree of ~ 4% at the peak (as shown in Fig. 4) probably corresponds to Earth’s polarization degree of ~ 13–20%.

We always show the observational results in (unconverted) Earthshine polarization degrees because there is a considerableuncertainty in the conversion from the Earthshine polarization degree to the Earth’s polarization degree. Furthermore, this is why we avoided relying on the absolute value of the Earthshine polarization degree in our discussion of the ocean signatures in Sect. 3. Instead, we discuss it in a relative manner (i.e., using P∕P_mean and $\Delta P/\overline{P}$) because the dependence on ϵ disappears as long as ϵ is constant. Although we believe the impact from ϵ is minimized in this manner, ϵ varies if we observe different lunar locations with different albedos, and thus may cause an undesired impact on our discussion. Below, we examine its impact on our discussion with respect to the nightly means and the hourly variations of Earthshine P. We note that we do not need to consider the phase dependence of ϵ because the phase angle of the depolarizing back-scattering on the Moon (the Earth–Moon–Earth angle) is always zero regardless of the lunar phase.

4.1.1 Impact on nightly-mean P

In Sect. 3.1, we treat the combined dataset of Earthshine P from both the waxing and waning lunar phases, and we found a correlation of P with the ocean fraction, as shown in Figs. 4 and 5. For Earthshine observations, we must point to the opposite (western and eastern) sides of the Moon between the waxing and waning phases because the opposite sides are illuminated by sunlight. If two different lunar locations for the waxing and waning phases have different albedos, an apparent difference of Earthshine P may be induced by different ϵ. Indeed, the near-infrared spectro-polarimetry of Earthshine by Miles-Páez et al. (2014) showed a contrast with a factor of 1.8 ± 0.3 between polarization degrees observed on two separate lunar locations. Because the ocean fraction is tied to the waxing-or-waning phases (on average, the ocean fraction is larger for the waning phase than for the waxing phase), there is a potential risk that the effect from different ϵ values may be wrongly interpreted as the effect from the oceans.

Figure 10 displays the same (P, α) dataset as Fig. 4; however, it distinguishes the waxing and waning phases by plot styles. It is barely recognizable that P is likely to be higher for the waning phase than for the waxing phase. However, the dependence of P on the ocean fraction as shown in Fig. 4 seems more obvious than the difference between the waxing and waning phases as shown in Fig. 10.

If the Earthshine P were severely affected by asignificant difference in ϵ between twodifferent lunar locations corresponding to the waxing- and waning-phase observations, the difference in Fig. 10 between the two phases should be more distinct. The indistinct waxing-or-waning difference is easily understood by accepting that Earthshine P is affected by the ocean fraction. Although the ocean fraction is larger in the waning phase than in the waxing phase on average, it is often similar between the two phases depending on the date and time. During our observations, the ocean fraction was ~ 15–40% for the waxing phase and ~20–45% for the waning phase. There is a large overlap in the ocean fractions.

Furthermore, Fig. 5 (top row) supports our interpretation. We consider two cases where Earthshine P is affected by two different ϵ values corresponding the waxing and waning phases.

First, we assume that Earthshine P is affected by different depolarizing factors, but it is not affected by ocean fractions. In this case, the data plots in Fig. 5 (top row) should be split into two levels, rather than showing a linear correlation. For instance, if ϵ is smaller (i.e., more depolarizing) in the waxing phase than in the waning phase, the plots should be distributed on a lower level (a smaller P∕P_mean) for the waxing phase and in a higher level (a larger P∕P_mean) for the waning phase.

Second, we assume that Earthshine P is affected by both of the different depolarizing factors and the ocean fractions. In this case, we can draw two separate regression lines for the waxing and waning phases in Fig. 5 (top row). In reality, however, the data points from both phases appear to roughly fall on a single regression line.

Based on the above discussions, we are convinced that the retrieved linear correlations in Fig. 5 (top row) do not result from the difference (if any) in the lunar depolarizing factors between waxing- and waning-phase observations, but are instead caused by the Earth’s oceans.

In this section, we discuss a possible impact of the difference in ϵ between the waxing and waning phases. Within one or the same phase (waxing phase or waning phase), we invested our best effort to observe the same lunar location as described in Sect. 2. However, the actually observed location may be slightly different from night to night because of the limited pointing accuracy, and this can cause night-to-night differences in ϵ. In contrast to the difference between the waxing and waning phases, the night-to-night differences in ϵ (caused by telescope pointing) is not coupled with Earth’s ocean fraction. Therefore, it is unlikely that the night-to-night differences cause the linear correlation shown in Fig. 5 (top row). Nonetheless, the deviations from the regression line in Fig. 5 (top row) may be caused in part by the night-to-night differences in ϵ.

Fig. 8 Time-series polarization degrees (left column), scene fractions (middle column), and polarization position angles (right column) for dates when significant hourly variation of polarization degree is detected. Left: Polarization degrees: Squares, circles, and triangles represent data in the J, H, and Ks bands, respectively. Open plots represent data points with |Θ − N| > 15°. The dashed line is the model curve for Earthshine polarization calculated based on the ocean, land, and cloud fractions. The applied lunar polarization efficiency (depolarizing factor, ϵ) is shown in the inset. Middle: Scene fractions: Ocean, land, and cloud fractions are exhibited as solid, dashed, and dotted lines, respectively (left y-axes). Fractions were calculated with concentrated weighting (see Appendix B.5). The crosses correspond to lunar elevation (right y-axes). Right: Polarization position angles: Solid lines in the right panels show the position angle normal to the scattering plane.

Fig. 9 Same as Fig. 8, but for dates when significant hourly variation of the polarization degree is not detected.

Fig. 10 Same as Fig. 4, except plot styles are distinguished by waxing-or-waning lunar phases. The filled and open plots represent observations in the waxing and waning phases, respectively.

4.1.2 Impact on hourly variation of P

Because different lunar locations may have different albedos, imperfect telescope tracking can lead to a false hourly variation in P that does not correspond to any of Earth’s properties. We attempt to estimate the possible variation of P caused by the shift of the observation location on the Moon. Although our target locations are not on the major maria, we occasionally recognize a dark patch on the reduced lunar images. From the visual inspection of the intensity (I) images on the three dates when a P variation is detected, we approximately determined the intensity contrasts between a dark region and the surrounding typical region (I_dark∕I_typ) in addition to the maximum area ratio of the dark region within the sampling region (S_dark∕S_smpl). Then, assuming that the albedo is proportional to the observed intensity, we estimate the possible highest albedo contrast for the night by $A_{\text{min}}/\overline{A}\cong 1-\left(S_{\text{dark}}/S_{\text{smpl}}\right)\left(1-I_{\text{dark}}/I_{\text{typ}}\right)$, where A denotes the effective albedo of the sampling region (A_min is the minimum and Ā is the mean)⁴. From the difference between Eq. (9) in Bazzon et al. (2013) for A_min and Ā, we have $\log ({ϵ}_{\text{min}}/\overline{ϵ})=-0.61\log (\overline{A}/A_{\text{min}})$. Although the albedo in Eq. (9) in Bazzon et al. (2013) is at wavelength of 602 nm, we assumed that the ratio of two albedos has a negligible wavelength dependence. Then, we estimate the relative variation of P by $\Delta P/\overline{P}=1-{ϵ}_{\text{min}}/\overline{ϵ}$. The results from the calculations are listed in Table 1.

For comparison, we determine $\Delta P/\overline{P}$ in the observed values. Time-series P (Fig. 8; left column) is fit by a linear function. We take the difference of P at the two ends of the fit line as ΔP and the average as $\overline{P}$. This derivation aims to avoid the overestimation of ΔP caused by a single extreme data point. The derived $\Delta P/\overline{P}$ in the observed values is summarized in Table 1.

The $\Delta P/\overline{P}$ estimated based on the variation in lunar depolarization is ~0.1 at the maximum. That variation is significantly smaller than the observed $\Delta P/\overline{P}$, which ranges from ~0.2 to ~1.4. Therefore, the variation in lunar depolarization cannot explain the observed variation of P.

Even if the depolarization variation caused by the tracking error contributes part of the hourly variations in the observed Earthshine P, it is very unlikely that the tracking error, which is independent of the scene fractions, can create clear correlations with the ocean fraction as shown in Fig. 7.

4.2 Telluric and telescope polarization

The term “telluric” in this article refers to the Earth’s atmosphere on the path from a celestial body to a ground-based observer. Telluric effects should be eliminated because we are only interested in the Earth’s properties as observed from outside the planet. Although it is usually assumed that telluric extinction does not polarize celestial light because of isotropy, extremely precise polarimetry by Bailey et al. (2008) indicated that telluric airborne dust can polarize celestial light because of the dichroic extinction caused by the dust. Nonetheless, the observed maximum polarization caused by this effect is as small as ~ 5 × 10⁻⁵, which is much smaller thanour observed P and its variations. Based on very sensitive solar polarimetry, Kemp et al. (1987) identified telluric polarization attributed to double scattering by aerosols and molecules in the Earth’s atmosphere. However, the measured polarization degree by this effect was ~8 × 10⁻⁶ at the maximum. Although other telluric polarizing sources may exist, we believe that a ~ 5 × 10⁻⁵ polarization degree, measured from the Canary Islands under a relatively strong effect from the Saharan dust (Bailey et al. 2008), provides a good upper limit to telluric polarization.

In addition, both of the above-mentioned telluric polarizing effects tend to be larger for a larger airmass (i.e., a lower elevation) (Kemp et al. 1987; Bailey et al. 2008). However, all our observed variations in P exhibited the opposite transition: on 2019 December 18, P increased with time while the Moon ascended; on 2020 January 3 and 2020 March 2, P decreased with time while the Moon descended (Fig. 8; left and middle columns). Therefore, we are convinced that polarization caused by telluric effects does not significantly affect our observations.

When the telescope pointing elevation is below 22°, part of the light beam incident onto the primary mirror is blocked by the enclosure wall. This breaks symmetry with respect to the telescope optical axis and can induce significant telescope polarization. The degree of the telescope polarization should increase as the pointing elevation decreases. However, P varied in the opposite sense, as described above. Thus, we exclude this effect from the causes of the observed variation in P.

5 Implications

5.1 Distinctiveness of polarimetric signature

One of the issues we should address is whether it is possible to distinguish between planets with an ocean and those without an ocean based on the observations of rotational variations in polarization. Comparison with near-infrared polarization of the Solar System objects would help the discussion; however, a lack of previous near-infrared polarimetry forces us to rely on results at the visible wavelengths.

Anti-correlation between the albedo and polarization degree of the reflected light is known as the Umov effect (Hapke 2005). The polarimetry of the integrated disk of the Moon showed different P at the peak of phase curves between the waxing and waning phases (Lyot 1929; Coyne & Pellicori 1970): the waning Moon was more polarized than the waxing Moon. The western (in selenographic coordinates) part of the Moon, illuminated in the waning phase, has a larger fraction of maria, and therefore it has a larger polarization than the eastern part. According to past observations (Coyne & Pellicori 1970), P at the effective wavelength of 534 nm was 10.9% at its peak in the waning phase, whereas it was 8.1% in the waxing phase. This implies that we will obtain $\Delta P/\overline{P}\cong 0.3$ when we observe a rotation of the Moon from outside the Earth-Moon system (above the lunar equator). These previous observations by Coyne & Pellicori (1970) were performed at several different wavelengths between 336 nm and 534 nm; the corresponding $\Delta P/\overline{P}$ values at different wavelengths do not exhibit an obvious wavelength dependence. Therefore, we expect that the near-infrared $\Delta P/\overline{P}$ of the Moon will not differ significantly from that at visible wavelengths (i.e., ~0.3).

Asteroid (4) Vesta is the only asteroid known to exhibit a convincing rotational variation in the polarization degree (Cellino et al. 2016) owing to the inhomogeneous albedo distribution. In previous visible polarimetry for Vesta, $\Delta P/\overline{P}$ was in the range 0.06–0.24 (Degewij et al. 1979; Broglia & Manara 1989; Lupishko et al. 1999; Wiktorowicz & Nofi 2015). These values were observed at phase angles of less than 20° when the polarization is negative (parallel to the scattering plane). It is not certain whether $\Delta P/\overline{P}$ in the negativepolarization regime is similar to that near the peak of positive polarization.

The previous observations of the Moon and Vesta suggest that for airless rocky bodies $\Delta P/\overline{P}$ are likely to be ~ 0.3 or less. For comparison, the near-infrared $\Delta P/\overline{P}$ of the Earth was 0.2–1.4 when P was highly variable (Table 1). Hence, we believe that the $\Delta P/\overline{P}$ of a planet with an Earth-like ocean fraction can be significantly larger than that of airless rocky planets.

For small icy bodies, we notice that some satellites exhibit a large difference in polarization depending on the central latitude (Rosenbush 2002; Ejeta et al. 2013). An analysis of the previous polarimetry of Jupiter’s satellite Callisto showed a $\Delta P/\overline{P}$ of ~ 0.3–0.5 at visible wavelengths (Rosenbush 2002). Spectro-polarimetry results of Saturn’s satellite Iapetus for both its leading and trailing hemispheres correspond to a $\Delta P/\overline{P}$ of ~ 0.8–1.5 at a wavelength ~900 nm (Ejeta et al. 2013), which is comparable to values from our Earthshine polarimetry in the near-infrared. These results suggest that the surface of the icy planetary bodies can have an extraordinarily distinctive albedo contrast that causes a large $\Delta P/\overline{P}$ comparable to that of a planet with a partial ocean. This should be considered when we interpret the polarization of planets near the outer edge of the habitable zone or beyond.

5.2 Feasibility estimate

The comprehensive feasibility evaluation of the polarimetric technique is beyond the scope of this work, though we briefly discuss it by referring to our previous work (Takahashi et al. 2017), in which we demonstrated the feasibility of the ground-based detection of a near-infrared spectro-polarimetric feature of water vapor in an exoplanetary atmosphere. We showed that the feature with a strength of ΔP_feature≅10% and continuum level of P_cont≅10% was detectable for 5–14 known exoplanets with a total exposure time of 15 hours using a 40-m class telescope such as the Extremely Large Telescope (ELT). In the estimate, we assumed that the high-contrast instrument suppressed stellar light down to 10⁻⁸–10⁻⁹, and its total throughput was 10%.

In the current case for ocean detection, the target signature is the polarization time variation of ΔP_time≅10% with a mean polarization level of $\overline{P}\cong 10\%$ (converted from the Earthshine polarization of ~2.5% near the quadrature phase assuming a lunar depolarization factor⁵ of ~ 0.25 and $\Delta P_{\text{time}}/\overline{P}\cong 1$), which is similar to the previous case for water vapor detection (i.e., ΔP_feature≅ΔP_time and $P_{\text{cont}}\cong \overline{P}$). Although the target signature in the previous work was a spectro-polarimetric feature with a feature width of Δλ≅ 0.05 μm, in the current work it is a broad-band signature. Hence, we can set Δλ≅ 0.15 μm assuming the coronagraph bandwidth to be ~10% of the H-band central wavelength (~ 1.6 μm). This reducesthe required exposure time to 15 × 0.05∕0.15 = 5 h. In the meantime, planets in the habitable zone around M-type stars, which are the main target of the in-development ground-based extremely large telescopes, are likely to be tidally locked (Kasting et al. 1993). Thus, a comparison of polarization between the waxing and waning near-quadrature phases will be effective for searching an inhomogeneously distributed ocean in the star-facing hemisphere, as long as the system is not face-on. In this case, the five-hour exposure time is sufficient to compare the two orbital phases.

Although forthcoming extremely large ground-based telescopes will not be optimized for polarimetry, some envisioned high-contrast instruments – namely, the Planetary Camera and Spectrograph (PCS, Kasper et al. 2021) for ELT, and the Planetary Systems Imager (PSI, Fitzgerald et al. 2019) for the Thirty Meter Telescope (TMT) – will have imaging polarimetry capabilities. It is worth seriously considering a search for an exoplanetary ocean using these instruments.

Although we demonstrated an estimate of the feasibility of ground-based detection, it does not imply that this technique cannot be implemented by space telescopes. The space-based time-series polarimetry for habitable-zone planets orbiting G-type stars is also worth considering when detecting the rotational variability of polarization caused by the existence of a partial ocean.

6 Conclusions

Our near-infrared polarimetry of Earthshine indicated the polarimetric signature of Earth’s oceans: we found a clear positive correlation of P with the ocean fraction on the Earthshine-contributing region (Figs. 4, 5 and 7); furthermore, we observed hourly variations of P in accordance with the rotational transition of the ocean fraction (Fig. 8). Although our simple model reproduced the observed hourly variation of P (Figs. 8 and 9) fairly well, modeling in a more sophisticated manner and inputting more appropriate Earth scene data (hourly time-resolved cloud maps) may resolve the exceptional observation-model disagreement and confirm the indicated ocean signature. The observed relative variation, $\Delta P/\overline{P}$, reached as large as ~0.2–1.4. An effective observation is estimated to be possible using a 40-m class ground-based telescope with a 5-h exposure. Therefore, we propose near-infrared polarimetry as a prospective technique for the detection of an exoplanetary ocean.

Acknowledgements

This work was supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI Grant Numbers 15K21296, 17K05390, and 21K03648; Tokubetsu Kenkyu Joseikin (2019–2021), funded by University of Hyogo; and the Optical and Near-Infrared Astronomy Inter-University Cooperation Program, funded by Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. This work was awarded the Best Presentation Award from the IAU Symposium 360 Science Organizing Committee held in Hiroshima, Japan 2021.

Note added in proof. Model simulations by Trees & Stam (2019) showed P phase curves of ocean planets with realistic ocean surfaces and Earth-like atmospheres. Our observations (Fig. 4) seem consistent with their simulations (at λ = 865 nm, the longest wavelength in their calculations) for partly cloudy ocean planets with f_c = 0.25 and 0.50.

Appendix A Additional tables

Appendix B Data reduction details

B.1 Image processing

The NIC’s raw images suffer from column-pattern (stripe-looking) noises whose strengths and spatial patterns differ by the detector’s quadrant and vary every frame. In the polarimetry mode, only a twin of ~ 150 × 430 pixels is used in the entire 1024 × 1024 pixels. This configuration allows us to read the pattern strength outside the polarimetric images and to subtract it from the whole frame. In the procedure, all constant offset counts (bias and dark counts) are subtracted. This procedure leaves virtually uniform counts (near-zero) outside the polarimetry windows (except some abnormal counts on bad pixels and occasional periodic noises described below), which guarantees the reliability of this method for subtracting the offset counts (see Fig. 4 in Takahashi et al. 2018).

After the subtraction of the column-pattern noises, periodic noises with an amplitude of a few analog-to-digital units (ADU) and a wavelength of several tens of pixels are often recognizable. They have very similar patterns between the left-half region and right-half region in a frame. Because astronomical images fall within the right-half region, we computed the fast Fourier transform for each column in the left-half region and extracted the strongest three frequencies (a constant component is included if its power is strong). The noise image reproduced using these frequencies was subtracted from the right-half region.

Next, flat-fielding was conducted. Data acquisition and image processing required to generate the flat frame is summarized as follows. (1) Illuminated screen on the enclosure wall was observed on 2018 May 7. (2) Data acquisition was conducted for eight different angles of the Cassegrain instrumental rotator (ϕ_insrot) ranging from − 135° to 180° with 45° steps. (3) For each single ϕ_insrot, we made 20 sequences of polarimetric observations. In total, we obtained 8 × 20 × 4 = 640 frames. (4) After column-pattern subtraction, all frames were averaged and normalized such that the mean count over the frame was unity.

Because the screen was illuminated by lamps with non-zero phase angles, the reflected light may be polarized. To generate a sensitivity map responding to unpolarized light, we obtained and combined data with symmetrically distributed angles of ϕ_insrot and ϕ_hwp. After flat-fielding, we fixed abnormal pixels using the Image Reduction and Analysis Facility (IRAF, Tody 1986, 1993) fixpix and cosmicray tasks.

Earthshine signals were veiled by strong scattered light from the day side of the Moon. The scattered light included in the sky background gets stronger as the position comes nearer to the day side of the Moon. The appropriate subtraction of the sky background intensity is necessary to measure the Earthshine intensities and derive polarization parameters.

We conducted the subtraction of the sky background in two steps. The sky background includes telluric emission and scattered Moonshine (light from the day side of the Moon). The scattered Moonshine has a spatial gradient in its intensities; it is stronger for a region nearer to the day side of the Moon. In the first step, blank sky frames (60″–90″ away from the target position) were used. For each single Earthshine frame, we selected a sky frame acquired with the same ϕ_hwp at the nearest time (unless the weather was too different). The selected sky frame was subtracted from the corresponding Earthshine frame. This procedure removes the sky background intensities superimposed upon an Earthshine image fairly well (Fig. B.1). However, because sky background intensities and their positional gradients on a blank sky region are not perfectly identical to those at the corresponding Earthshine region⁶, some residual counts must exist after the first background subtraction.

Fig. B.1 Intensity profiles before and after subtraction of sky background. The solid lines indicate profiles along the detector y-axis before the subtraction (top line), after the first subtraction (middle), and second subtraction (bottom). The top line is vertically offset by −1000 ADU for visibility. These profiles are averages of the central 20 x-positions. The left half area (y < ~ 200 pix) is the Moon (Earthshine), and the other half is sky. The dashed line is the fit line for the sky region. The corresponding image data after the second subtraction is shown in Fig. 2 (left).

The second background subtraction was executed. We fit a linear function with respect to the y-axis (RA) to the sky area on an Earthshine frame after the first background subtraction. The fit function was extrapolated toward the Moon and subtracted from the whole frame⁷ (Fig. B.1). Linear extrapolation and subtraction is common in previous Earthshine data reduction (e.g., Qiu et al. 2003; Hamdani et al. 2006; Bazzon et al. 2013). Frames after this second background subtraction are the final processed intensity images that are ready for photo-polarimetry. A sample is displayed in Fig. 2 (left). The sky region in the background-subtracted images is fairly flat (Fig. B.1). In these images, we did not recognize any complicated pattern thatimplies stray light in the instrument.

B.2 Derivation of Stokes parameters

A set of normalized Stokes q and u was derived from a sequence of observationswith four different ϕ_hwp. At least two different paths exist to proceed from a set of 2D intensity maps to representing two values of q and u. In the “0D” method, we first averaged the counts over the defined sampling region on the Earthshine intensity map (reduced data dimension from 2D to 0D), and then, q and u were calculated. In the “2D” method, we generated 2D maps of q and u, and then the averaged values of q and u were extracted from the sampling region. Below, we describe these two methods in more detail and compare them.

In the 0D method, we define the sampling region as the region enclosed by a parallelogram with widths of 100 pixels (~ 16″) along the x-axis (DEC) and 50 pixels (~ 8″) along the y-axis (DEC), as shown in Fig. 2 (left). A buffer with a distance of 20 pixels (~ 3″) along the y-axis was taken between the automatically detected lunar edge and the sampling region. We computed the average and standard deviation of sampled intensity values with 3σ clipping, after compressing⁸ the 2D data into 1D data that are a function of y.

Normalized Stokes parameters defined in the instrumental coordinates q₀ and u₀ are derived as follows: $$\begin{array}{lll}{R}_q\hfill & =\hfill & \sqrt{\frac{I_{\text{e},0^{°}}}{I_{\text{o},0^{°}}}\frac{I_{\text{o},{45}^{°}}}{I_{\text{e},{45}^{°}}}},\hfill \\ R_u\hfill & =\hfill & \sqrt{\frac{I_{\text{e},{22.5}^{°}}}{I_{\text{o},{22.5}^{°}}}\frac{I_{\text{o},{67.5}^{°}}}{I_{\text{e},{67.5}^{°}}}},\hfill \\ q_0\hfill & =\hfill & \frac{1-R_q}{1+R_q},\hfill \\ u_0\hfill & =\hfill & \frac{1-R_u}{1+R_u},\hfill \end{array}$$

where I denotes the intensity averaged over the sampling region. The first subscript of I specifies where the I value comes from, namely, ordinary (o) or extraordinary (e) rays. The second subscript is ϕ_hwp, with which the corresponding I frame was observed. The double rationing in Eqs. (B.1) and (B.2) removes any multiplicative effect common between ordinary and extraordinary rays (such as isotropic telluric extinction) and different optical throughputs between the paths of ordinary and extraordinary rays.

Because no significant instrumental (including telescope) polarization was detected from the observations of unpolarized standard stars (Takahashi et al. 2018), we did not perform the subtraction of the instrumental polarization. Errors in a single set of q₀ and u₀ were estimated from the standard deviations of I following the rule of error propagation. The values of q₀ and u₀ derived from the 0D method are shown in Fig. B.2 as open plots for the J-band data on 2020 January 3.

In the 2D method, the calculations of Eqs. (B.1)–(B.4) were executed as 2D image processing rather than the calculations of the averaged scaler values (over the sampling region). Misalignment and different spatial resolutions between the I frames may produce spurious polarization. Therefore, before the image processing of Eqs. (B.1)–(B.4), we aligned the I frames by shifting them along the y-axis based on the detected y position of the lunar edge, and we smoothed these images using a median filter with a window size of 10 × 10 pixels (1.6″ × 1.6″, which corresponds to a typical seeing size at the observatory). Typical q₀ and u₀ images are displayed in Fig. 2 (middle and right).

Sampling and averaging values from q₀ and u₀ images were taken in exactly the same manner as for the I images. Errors in a single set of q₀ and u₀ were estimated as standard deviations in sampled values after 2D-to-1D compression. The values of q₀ and u₀ derived from the 2D method are shown as filled plots in Fig. B.2.

A quick glance at Fig. B.2 tells us that results from the 0D and 2D methods are almost identical. This is true for data on other dates except for limited sets with very low signals. If the misalignment and/or resolution mismatch between images in the 2D method caused significant errors, such errors should be minimized in the 0D method. The coincidence of the results from the two methods implies such errors were insignificant.

Estimated errors from the 0D method (shown as pale bars in Fig. B.2) seem overestimated because they are considerably larger than the scattering of values in a series of several (q₀ u₀) sets. This overestimation may be attributed to true distribution of intensities on the Moon (i.e., reflectivity distribution, which should not be counted as an error) and/or certain types of systematic errors removed by the double rationing in Eqs. (B.1)–(B.2). Errors from the 2D method (shown as thick bars in Fig. B.2) may be underestimated because they appear to be slightly smaller than local scattering. This underestimation possibly stems from the smoothing and compressing before calculating standard deviations.

Fig. B.2 Time-series Stokes q0 (circles) and u0 (triangles) on 2020 January 3 for J band. The open and filled plots correspond to the results from the 0D and 2D methods (see text), respectively.

For further data reduction and discussion, we used q₀ and u₀ derived from the 2D method. The results from the two methods did not differ significantly; however, the 2D maps of q₀ and u₀ were useful to investigate possible systematic errors in the spatial distributions of q₀ and u₀. For example, if we had found a gradient in q₀ or u₀ along the y-axis, it could have been attributed to the imperfect sky background subtraction and/or the effect of glancing views of the lunar edge. Fortunately,we did not find such a gradient.

Errors derived from the 2D method may be underestimated as described. We did not use these errors for error estimates in degrees and position angles of polarization. Instead, we used them as weights when averaging multiple sets of q₀ and u₀.

B.3 Derivation of polarization degrees and position angles

The polarization degree (P₀) and polarization position angle (Θ₀) were derived by $$\begin{array}{lll}\hfill & \hfill & P_0=\sqrt{{\overline{q}}_0^2+{\overline{u}}_0^2},\hfill \\ \hfill & \hfill & \tan (2{\Theta }_0)={\overline{u}}_0/{\overline{q}}_0,\hfill \end{array}$$

where ${\overline{q}}_0$ and ū₀ are averages of multiple sets of q₀ and u₀ calculated with 3σ clipping and weighting based on estimated errors described in the previous subsection. When we derive P₀ and Θ₀ as nightly means, ${\overline{q}}_0$ and ū₀ are averages of all data on the night. When we investigate hourly variations of P₀ and Θ₀, ${\overline{q}}_0$ and ū₀ are averages for every 30-minute bin.

Errors in P₀ and Θ₀ were estimated by $$\begin{array}{lll}{\sigma }_P\hfill & =\hfill & \frac{\sqrt{q^2{\sigma }_q^2+u^2{\sigma }_u^2}}{P_0},\hfill \\ {\sigma }_{\Theta }\hfill & =\hfill & \frac{{180}^{°}}{2\pi }\frac{{\sigma }_P}{P_0},\hfill \end{array}$$

where σ_q and σ_u are the standard deviations of q₀ and u₀, respectively.

When σ_P is relatively large, P₀ may be positively biased. Therefore, we applied the following correction developed by Plaszczynski et al. (2014): $$P=P_0-{\sigma }_P^2\frac{1-{\text{e}}^{-P_0^2/{\sigma }_P^2}}{2P_0}.$$(B.9)

Because Θ₀ is defined in the instrumental coordinates, we converted it to equatorial coordinates using $$\Theta =-({\Theta }_0-{\varphi }_{\text{inspa}}),$$(B.10)

where Θ denotes the position angle of polarization measured counter-clockwise from the equatorial north.

Although the possible occurrence of depolarization and a Θ offset as instrumental effects was suggested in a previous performance evaluation (Takahashi 2019), we did not attemptto correct these effects. The polarization efficiency (depolarizing factor) may be down to ~ 0.9, and the Θ offset may be up to ~ 1°. The absenceof such correction will not affect the essence of the discussion in this work.

B.4 Derivation of the mean phase curves of polarization

The mean phase curves of the polarization degree (P_mean) in Figs. 4 and B.3 were retrieved based on Eq. (13.11) in Hapke (2005)⁹. This equation considers polarized specularly reflected light and unpolarized multiply scattered light, with an intention of explaining the positive branch (polarized perpendicularly to the scattering plane) of the polarization phase curves of solid planetary bodies such as the Moon and asteroids. Although Earth differs from these solid bodies in various respects, the overall shape of Earth’s polarization phase curve (Dollfus 1957; Takahashi et al. 2012, 2013; Bazzon et al. 2013; Sterzik et al. 2019) is similar to the positive branch of polarization phase curves of the Moon (Lyot 1929; Coyne & Pellicori 1970) and other solid bodies (as summarized in Ito et al. 2018). These polarization phase curves have their peak polarization degree at an α (phase angle) between 90° and 150°, and it approaches zero when α goes toward either 0° or 180°. As described in the main text (Sect. 3.1), the purpose of this fitting is to extract the polarization phase curve of the typical Earth scene. Any function reproducing the general phase angle dependence of the polarization degree is sufficient for this purpose.

In the equation, the phase curve is determined by refractive index, n, and single scattering albedo, w. We fit the equation by scanning w and a scaling factor s, whereas n was fixed to be 1.3, the value for water. If Θ differs from N by more than 15°, the corresponding P data points were excluded from the fitting. Data points with α < 50° were also excluded from the fitting because some theoretical works for Earth-like planets expect a rainbow feature (enhanced polarization by water droplets) at α ~ 20°–50° (Bailey 2007; Stam 2008; Zugger et al. 2010), but the applied equation does not consider such an effect. Very recently, therainbow feature was detected by Earthshine polarimetry in the visible wavelengths at phase angles of ~ 30°–45° (Sterzik et al. 2020). The excluded points are shown as crosses in Figs. 4 and B.3. The obtained P_mean curves are displayed as dashed lines in those figures. The sole purpose of this fitting procedure is to draw the mean phase curves. The derived parameters such as w and s do not provide insightful information.

B.5 Derivation of scene fractions

The ocean, land, and cloud fractions in the Earthshine-contributing region were derived based on data from the Moderate Resolution Imaging Spectroradiometer (MODIS) aboard the Aqua and Terra satellites. Namely, data product MCD12C1 (Friedl & Sulla-Menashe. 2015) for the year 2018 was referred to for the surface type (oceans or lands) classification, and product MOD08_D3 (Platnick et al. 2015) was used as the reference of the cloud distribution for each observation date. Moreover, we retrieved the cloud’s top height from MOD08_D3 and used it for the classification of the clouds (Fig. 6). The time resolution ofMOD08_D3 is one day, and therefore, we could not consider the exact cloud distribution at a specific time and its hourly movement.

The surface dataset (MCD12C1) is a collection of grid data with a grid size of 0.05° (longitude) × 0.05° (latitude). Each grid is classified as one of 17 surface types. In this work, we defined the “water bodies” type as oceans, and all other types as lands. The grid size of the cloud dataset (MOD08_D3) is 1° × 1°. A single cloud grid corresponds to 20 × 20 surface grids (hereafter, “sub-grids”). Within these 400 sub-grids, we counted the number of ocean sub-grids (n_o) and that of land sub-grids (n_l). For each cloud grid, we retrieved the value of “cloud_fraction_mean” and treated it as the local cloud fraction (f_{c, local}). Considering the cloud fraction, the effective counts of ocean, land, and cloud sub-grids within a single cloud grid are derived by $$\begin{array}{lll}{n}_{\text{o},\text{eff}}\hfill & =\hfill & (1-f_{\text{c},\text{local}})n_{\text{o}},\hfill \\ n_{\text{l},\text{eff}}\hfill & =\hfill & (1-f_{\text{c},\text{local}})n_{\text{l}},\text{and}\hfill \\ n_{\text{c},\text{eff}}\hfill & =\hfill & f_{\text{c},\text{local}}(n_{\text{o}}+n_{\text{l}}),\hfill \end{array}$$

respectively.

In essence, we summed n_{o, eff}, n_{l, eff}, and n_{c, eff} for all cloud grids in the Earthshine-contributing region, and then we derived the final fractions, f_o, f_l, and f_c. We have a relation f_o + f_l + f_c = 1. When summing n_{x, eff}, we applied grid-to-grid weighting by ${\left(\cos ({\phi }_{\text{ss}})\cos ({\phi }_{\text{sl}})\right)}^d$ after the correction of different grid area sizes¹⁰, where φ_ss and φ_sl represent the angular separations of the grid from the subsolar and sublunar points around Earth’s center, respectively. This function yields the heaviest weight at the glint point (midpoint between subsolar and sublunar points).

When power d is larger, the function becomes steeper, which makes the weights more concentrated around the glint point. We applied two different d values of 1 and 10. When d = 1 is employed,the weighting function considers only tilts of a grid as viewed from the Sun and Moon. We refer to this method as “normal weighting”.

From a remote viewpoint, Earth’s specular reflection is seen from a limited region ~ 30° (in angular separation about Earth’s center) wide around the glint point (Williams & Gaidos 2008). There is a possibility that Earthshine polarization is more sensitive to fractions in a limited region around the glint point than those in the Earthshine-contributing region as a whole because the former is more directly related to whether the sea glint is seen from the Moon. Considering this possibility, we applied d = 10, which makes a full width at half maximum (FWHM) of the weighting function of ~ 30° around the glint point. We refer to this calculation method as “concentrated weighting”.

The two versions of scene fractions are listed in Table A.1 where values in parentheses are from concentrated weighting. The values in the table are the nightly means, that is, averages over the observation mid-time ±2 hours. Figures B.3 and B.4 are drawn with the fractions by normal weighting, whereas Figs. 4 and 5 are by concentrated weighting. Both versions indicate common trends: the ratio P∕P_mean has a clear positive correlation with the ocean fraction, a less clear negative correlation with the land fraction, and no clear correlation with the cloud fraction. Slopes (a) of regression lines for the ocean fraction are shallower in Fig. 5 than in Fig. B.4 because the fractions are more widely distributed in the concentrated weighting version. We note that P∕P_mean values in the both figures are identical.

Fig. B.3 Same as Fig. 4, except ocean fraction was calculated with normal weighting.

Fig. B.4 Same as Fig. 5, except ocean, land, and cloud fractions were calculated with normal weighting.

Although the two different methods did not make any essential difference in the examination of the observational results as nightly means, we found that fractions calculated by concentrated weighting are more favorable to explain the observed hourly variations in P. This point is described in Appendix C.

Appendix C Observation-model comparison

Among the long-time datasets on the six dates, those on 2019 December 18, 2020 January 3, and 2020 March 2 exhibited significant hourly variations (Fig. 8; left column), whereas those on the other dates did not (Fig. 9; left column). The variations in P on 2020 January 3 and 2020 March 2 appear synchronized with the ocean fraction, whereas those on the other dates are difficult to explain in such a simple manner.

We hereby attempt to reproduce the observed hourly variations (including non-variations) based on scene fractions at the time. A model of planetary reflected light (Williams & Gaidos 2008) is referred to. Equation (6) of Williams & Gaidos (2008) is simplified as $$F_{\text{Earth}}=f_{\text{c}}{A}_{\text{c}}+f_{\text{l}}{A}_{\text{l}}+\frac{p_{\text{wav}}}{f_{\Omega }}{f}_{\text{o}}{A}_{\text{o}},$$(C.1)

where F_Earth denotes the flux of Earth’s reflected light normalized by the incident flux, f_x denote the disk area fractions of each scene type, and A_x denote the albedos of each scene type. Subscripts c, l, and o represent clouds, lands, and oceans, respectively. Originally, another term for surface ice was included in Williams & Gaidos (2008), which is merged into lands in this work. Specularly reflected light is scattered into a solid angle f_Ω2π (str). The parameter p_wav is the probability that waves are properly oriented for sending the specularly reflected light into the observer.

We input the fractions calculated with the concentrated weighting (shown in Figs. 8–9, middle columns) into f_x. Although calculations of Eq. (C.1) were performed for every grid on the Earth surface model in Williams & Gaidos (2008), we simply calculated itfor a single set of fractions representing the properties of the Earthshine-contributing region at a specific time. The ocean albedo A_o is calculated using Fresnel equations with a refractive index (n) of 1.3; hence, it depends on phase angle α. Following Williams & Gaidos (2008), we apply A_c = 0.6 and A_l = 0.3, which are assumed to be independent of α, and f_Ω = 10⁻⁵. Although p_wav is not well known, p_wav = 3 × 10⁻⁵ was set by trial and error.

The polarized component within F_Earth can be expressed by $$P_{\text{Earth}}{F}_{\text{Earth}}=P_{\text{c}}{f}_{\text{c}}{A}_{\text{c}}+P_{\text{l}}{f}_{\text{l}}{A}_{\text{l}}+P_{\text{o}}\frac{p_{\text{wav}}}{f_{\Omega }}{f}_{\text{o}}{A}_{\text{o}},$$(C.2)

where P_Earth denotes Earth’s polarization degree, and P_x denote the polarization degrees of each scene type with the same subscripts as in Eq. (C.1).

Similarly to A_o, the ocean polarization P_o is obtained from Fresnel equations with n = 1.3, which give a peak polarization of 100% at a phase α of ~ 105°. Although scattered light from the clouds and lands was treated as completely unpolarized in Williams & Gaidos (2008), we introduce a small polarization to meet the observed polarization. P_l and P_c are expressed by scaled Fresnel equations with n values of 1.4and 1.3 and the scaling parameters of 0.1 and 0.03, respectively. They obtain peak polarization degrees of 10% and 3% at α of ~ 110° and ~ 105° for lands and clouds, respectively. Finally, P_Earth is derived by dividing Eq. (C.2) by Eq. (C.1).

Before we compare the modeled P_Earth with observed lunar Earthshine polarization, we must implement the depolarization at the lunar surface. As described in Sect. 4.1, the depolarizing factor (polarization efficiency, ϵ) is not known, especially for the near-infrared region. When we extend Eq. (9) in Bazzon et al. (2013) to the near-infrared wavelengths (1.2–2.2 μm) with typical highland albedos (0.15–0.25 in visible wavelengths), ϵ ~ 0.2–0.3 is deduced. Because we may have observed different lunar locations for different dates, we select one of 0.2, 0.25, and 0.3 as ϵ to match the model to the observedP separately for each date. The model Earthshine polarization is presented as dashed lines in the left column of Figs. 8–9 with the selected ϵ. We see a resemblance between the time-variation of the modeled P and that of f_o, which signifies that it is mainly f_o that controls the variation of the modeled P.

For 2020 January 3 and 2020 March 2, we see excellent agreement between the observed hourly variation of P and the modeled P (Fig. 8). For 2019 November 21, 2019 December 19, and 2020 April 29, the insignificant variations in the observed P appear to be consistent with a small variation in the modeled P (Fig. 9). Large variations in P can be explained by large variations inthe ocean fraction (Δf_o > 25%).

For 2019 December 18, the model fails to reproduce the observed P (Fig. 8). A continuous increase in P was observed, whereas the modeled P peaks at UT~17 hour in accordance with the ocean fraction. An explanation for this disagreement is the possible inaccuracy of the referred cloud distribution. As described in Appendix B.5, the time resolution of the MODIS cloud distribution is one day (limited by the FOV of the instrument and orbital frequency of the satellite); thus, MODIS data may not reflect the exact cloud distribution at specific times. We found a mass of cloud near the glint point at UT 18–19 hour on the date in the visualized MODIS data. Thus, there is a possibility that a small shift (~ 1000 km) of cloud locations may lead to significant differences in the calculated cloud fractions for the Earthshine-contributing region. Although it is better to refer to an hourly time-resolved cloud distribution to test this possibility, we leave this for future work. More sophisticated modeling and alternative cloud data sources found in Tinetti et al. (2006), Montañés-Rodriguez et al. (2005), and Montañés-Rodríguez et al. (2006) may help improve our model.

Throughout the discussion of the hourly variation in P, we have always been referred to the scene fractions calculated by concentrated weighting instead of those by normal weighting. This is because we found that models based on the former fractions match better with the observed P than the latter. The model from normal weighting cannot reproduce the large variation on 2020 March 2 (variation of the ocean fraction is too small), though it explains the observed significant variation on 2020 January 3, and the insignificant variations on the other three dates fairly well. For 2019 December 18, we noticed similar discrepancies in the models from both weighting methods.

References

All Tables

All Figures

	Fig. 1 Views of cloud-free Earth from the Moon at different observation times. At each panel, the illuminated hemisphere is the Earthshine contributing region. These images were created with the Earth and Moon Viewer3 developed by John Walker.
In the text

	Fig. 2 Intensity (left), Stokes q0 (middle), and u0 (right) images observed on 2020 January 3. North is right, and east is up. The FOV is ~ 19″ × 64″ (smaller than original FOV because of trimming). The parallelograms are the sampling regions. The values of q0 and u0 on the sky are scattered because of division of ~0 by ~ 0.
In the text

	Fig. 3 Observed position angles of Earthshine polarization as plotted against position angles normal to scattering plane.
In the text

Fig. 4 Earthshine polarization degrees (P) in J (left), H (middle), and Ks (right) bands, plotted against Sun–Earth–Moon phase angle (α). The ocean fraction was calculated with concentrated weighting (see Appendix B.5 for details on the derivation of the fractions). The dashed lines represent the polarization phase curve for the typical Earth scene. They are derived from fitting a curved line to all data points (with some exceptions described below) with free parameters w (single scattering albedo) and s (scaling factor). The crosses represent data points corresponding to |Θ − N| > 15° (where Θ denotes the position angle of polarization and N denotes the position angle normal to the scattering plane) or α < 50°, which were excluded from the fitting (see Appendix B.4 for details on the fitting). These figures give an impression that P for a larger ocean fraction (plots with a darker color) tends to be larger than those for a smaller ocean fraction at a similar α, which suggests that the contribution from the sea glint (specular reflection) enhances the polarization degree of Earth.

In the text

Fig. 5 Dependence of polarization degree on ocean (top), land (middle), and cloud (bottom) fractions (in J, H, and Ks bands from left to right). Each polarization degree (P) in Fig. 4 is divided by typical polarization (Pmean: dashed line in Fig. 4) at the phase angle, and then plotted against the ocean, land, or cloud fraction (fo, fl, or fc, respectively). The filled and open plots correspond to observations in the waxing and waning phases, respectively. The dashed lines are regression lines of the form af + b with a correlation coefficient r. The crosses correspond to those in Fig. 4 and were excluded from the linear regression. The fractions was calculated with concentrated weighting.

In the text

	Fig. 6 Dependence of polarization degree on high- (top), middle- (middle), and low-cloud (bottom) fractions (in J, H, and Ks bands from left to right). The legends are the same as those in Fig. 5.
In the text

	Fig. 7 Same figure as Fig. 5, except the data source is the time-resolved data on 2019 November 21, 2019 December 18, 2019 December 19, 2020 January 3, 2020 March 2, and 2020 April 29. The data on the same date are connected by lines.
In the text

Fig. 8 Time-series polarization degrees (left column), scene fractions (middle column), and polarization position angles (right column) for dates when significant hourly variation of polarization degree is detected. Left: Polarization degrees: Squares, circles, and triangles represent data in the J, H, and Ks bands, respectively. Open plots represent data points with |Θ − N| > 15°. The dashed line is the model curve for Earthshine polarization calculated based on the ocean, land, and cloud fractions. The applied lunar polarization efficiency (depolarizing factor, ϵ) is shown in the inset. Middle: Scene fractions: Ocean, land, and cloud fractions are exhibited as solid, dashed, and dotted lines, respectively (left y-axes). Fractions were calculated with concentrated weighting (see Appendix B.5). The crosses correspond to lunar elevation (right y-axes). Right: Polarization position angles: Solid lines in the right panels show the position angle normal to the scattering plane.

In the text

	Fig. 9 Same as Fig. 8, but for dates when significant hourly variation of the polarization degree is not detected.
In the text

	Fig. 10 Same as Fig. 4, except plot styles are distinguished by waxing-or-waning lunar phases. The filled and open plots represent observations in the waxing and waning phases, respectively.
In the text

Fig. B.1 Intensity profiles before and after subtraction of sky background. The solid lines indicate profiles along the detector y-axis before the subtraction (top line), after the first subtraction (middle), and second subtraction (bottom). The top line is vertically offset by −1000 ADU for visibility. These profiles are averages of the central 20 x-positions. The left half area (y < ~ 200 pix) is the Moon (Earthshine), and the other half is sky. The dashed line is the fit line for the sky region. The corresponding image data after the second subtraction is shown in Fig. 2 (left).

In the text

	Fig. B.2 Time-series Stokes q0 (circles) and u0 (triangles) on 2020 January 3 for J band. The open and filled plots correspond to the results from the 0D and 2D methods (see text), respectively.
In the text

	Fig. B.3 Same as Fig. 4, except ocean fraction was calculated with normal weighting.
In the text

	Fig. B.4 Same as Fig. 5, except ocean, land, and cloud fractions were calculated with normal weighting.
In the text