© ESO 2019

1 Introduction

With the discovery of the first potentially habitable planet in the solar neighborhood by Anglada-Escudé et al. (2016), the quest for remote detection of signatures of life in other worlds has become one of the most exigent problems in modern astrophysics. Remotely detectable signs of life are being carefully assessed (see the reviews of Schwieterman et al. 2018; Catling et al. 2018), and observational prospects and requirements are being mapped onto the next generation of astronomical facilities (see the review of Fujii et al. 2018). Facing colossal difficulties, further progress is needed in all areas encompassing astrobiology, theoretical concepts and frameworks, methodology, models, and instrumental capabilities.

For the time being, Earth is the only inhabited planet we know. Planet Earth is therefore a benchmark object to infer biosignatures of life as we know it today. The signature of Earth seen as an exoplanet, thus from afar, depends strongly on the local illumination and viewing geometries and cannot be derived from Earth remote-sensing observations such as taken by low-Earth-orbit satellites. One way to study Earth from afar is to observe Earthshine. Earthshine is sunlight scattered by Earth and reflected from the lunar surface back to Earth, where it can be observed with ground-based astronomical facilities (see the reviews of Arnold 2008; Vazquez & Pallé 2010). Biosignatures such as the vegetation red edge (VRE) and trace gases of biotic origin such as O₂ and CH₄ can be inferred from spectroscopic observations (Turnbull et al. 2006). Different surface sceneries of Earth are probed as the Sun–Earth–Moon phase angle α and corresponding viewing geometry changes. Temporal changes in the associated light curve contain further information about surface and atmospheric properties of Earth, mainly through spatial and temporal modulation of Earth’s albedo. For example, high-cadence photometric time series measurements from the Deep Space Climate Observatory imaging Earth from afar can be used to reconstruct planetary rotation, cloud patterns, and surface types without prior detailed knowledge of its properties (Jiang et al. 2018).

However, detection and reliable extraction of biosignatures from Earth’s reflectivity spectra remains difficult. In the case of Earthshine, the final transmission of light through Earth’s atmosphere contaminates the signatures in the Earthshine’s spectral flux through extinction by scattering and absorption and has to be carefully corrected for. Extinction by Earth’s atmosphere does not, however, usually change the degree of polarization P of light coming from an astronomical source, as both polarized and unpolarized fluxes are affected to the same proportion. Polarimetry implemented as an intrinsically differential measurement technique should therefore overcome most telluric calibration issues and allow more reliable measurements from the ground. Polarization spectra, in general, also have a higher diagnostic value, which allows us to characterize their scattering and reflecting particles and surfaces more fully than in flux spectra. Therefore, we explored spectropolarimetry of Earthshine at optical wavelengths to improve the extraction of characteristic properties of Earth contained in the sunlight that Earth reflects (Sterzik et al. 2012). The method allows us to determine the fractional contribution of clouds and ocean surfaces and can distinguish visible areas of vegetation as small as 10% by comparing two datasets from different epochs with different aspects of Earth. While the measured O₂ -A band strength and the VRE feature were fully compatible with the results from spectral polarization models for Earth-type exoplanets (Stam 2008), the shape of the polarized spectral continuum remains unexplained: it is significantly flatter in the red spectral part than expected from the models.

Relatively flat polarization spectra of Earthshine were corroborated by Takahashi et al. (2013), who presented a series of polarization spectra covering phase angles between 49° and 96° that were recorded in five consecutive nights. Miles-Páez et al. (2014) extended the wavelength regime for spectropolarimetry of Earthshine to the near-infrared, and found particularly high degrees of polarization for the H₂ O band around 1.12 μm.

Classic observations of Earthshine polarization were presented by Dollfus (1957), covering phases from around 30° to 140° in the V band. The main characteristics of the polarization phase curve of Earth were qualitatively and quantitatively explained. Dollfus (1957) found a steady increase in the fractional polarization of Earthshine from about P ≈ 2% around a phase of 30° to around 10% peak values for 100°, and a decrease for larger phase angles. Dollfus (1957) noted a wavelength dependence of depolarization caused by backscattering at the lunar surface and concluded that fractional polarization of the light scattered by Earth is around P_V ≈ 33% around quadrature.

Because Earthshine is reflected by the lunar surface, ground-based observations have to be corrected for the resulting change of the degree of polarization. Indeed, Bazzon et al. (2013) probed Earthshine polarization in four B, V, R, I photometric bandpasses for phase angles between 30° and 110°. These authors demonstrated that the reflection by different regions on the Moon with different albedos (Highlands and Mare) has a significant and distinct impact on the polarization measured in Earthshine. To mitigate this problem, they introduced a method to correct for lunar depolarization efficiency that depends on wavelength and lunar albedo. Their measurements convincingly demonstrated that the lunar surface albedo needs to be taken into account when calibrating and interpreting the absolute values of Earth’s degree of polarization derived from Earthshine, and that lunar albedo effects may be of the same order as those caused by variations of the albedo of Earth.

In the last decade, a substantial amount of theoretical work has been carried out to model and simulate the polarization spectrum of Earth. Stam (2008) employed a method to approximate an inhomogenuous surface albedo and clouds by using weighted sums of light reflected by horizontally homogenous planets with a specific surface reflection function (Fresnel or Lambert), covered by an anisotropic Rayleigh scattering atmosphere containing Mie scattering clouds with a fixed optical depth (equals to ten). These calculations covered a wide range of phase angles and a wide range of surface conditions. Results can easily be retrieved from look-up tables and enable the systematically exploration of a multiparameter problem. Karalidi & Stam (2012) extended the methodology to horizontally inhomogeneous clouds and surfaces of Earth-like exoplanets. Karalidi et al. (2011, 2012a) further incorporated effects of liquid and ice water clouds on the degree of polarization and their phase dependence. An alternative Monte Carlo approach simulates a whole Earth-type planet in full spherical geometry, and was presented by García Muñoz & Mills (2014) and García Muñoz (2015).

The three-dimensional Monte Carlo vector radiative transfer code MYSTIC (Mayer 2009; Emde et al. 2010) has been extended to fully spherical geometry, so that it allows us to simulate Earthshine polarization spectra (Emde et al. 2017). The MYSTIC code is one of the solvers included in the libRadtran package¹ (Emde et al. 2016), which is widely used for Earth remote-sensing applications. The influence of aerosols, clouds, and the potential impact of sunglint was studied to explain in detail their impact on two spectra published by Sterzik et al. (2012). They incorporated three-dimensional fields of cloud properties (cloud cover, liquid, and ice water content) from the ECMWF IFS weather forecast model² at the time of observations and a two-dimensional surface albedo map derived from MODIS satellite observations, and simulated high spatial and high spectral resolution maps of the Stokes vectors across the surface of Earth. The simulations approximate the observed polarization spectra much better than the weighted sum spectra of the horizontally homogeneous planets (Sterzik et al. 2012).

The purpose of this paper is to extend the analysis of Earthshine polarization spectra to a much larger sample of high-quality observations. In a monitoring campaign, around 50 polarization spectra of Earthshine have been recorded between April 2011 and February 2013, all of which were observed with FORS2 at the Very Large Telescope (VLT). These spectra cover phase angles of Earth between 50° and 135°. Suitable parameters derived from these spectra allow statistical analysis of the temporal and spectral variation of the degree of polarization. In this paper, we discuss 34 spectra of the highest quality. One goal of this work is to establish a set of observables from Earthshine spectropolarimetry that can be used to characterize the surface of Earth and atmosphere at the time of observations, minimizing the uncertainties caused in particular by lunar (de)polarization efficiencies. We extensively use the models of Stam (2008) to compare them with the statistical properties in our sample and infer global trends. However, we do not intend to simulate and explain specific and individual datasets because they require dedicated modeling, which is beyond the scope of this paper.

The process of acquiring and reducing our datasets is a challenging task. Essentially, this process consists of spectropolarimetric measurements of a fraction of the visible lunar disk not illuminated by the sun. Observationally, the lunar disk is an extended target that is moving at a quickly changing, nonsidereal velocity, which cannot automatically be tracked by the VLT. In addition, the Earthshine signal that is reflected by the dark fraction of the lunar disk is contaminated by the (polarized) signal of the Moonshine (the part of the lunar disk that is illuminated by the sun). This contamination increases with increasing lunar phase.

This paper is organized as follows. In Sects. 2 and 3 we explain the procedure to obtain accurate data. In Sect. 4 we elaborate on our results. We present our approach to correct the Earthshine polarization for lunar depolarization, and derive the polarization phase curve of Earth in B, V, R, I spectral bandpasses. We derive polarization color ratios, the differential polarized vegetation index and the equivalent width of the O₂ -A band for our observations, look into the short-term variability of these parameters, and compare them with literature and model results. Finally, we summarize our conclusions and outlook in Sect. 5.

2 Observations

In this section we describe the instrument and instrument settings for which we have obtained Earthshine measurements, as well as the special procedures that we adopted for data acquisition and data reduction.

2.1 Instrument and instrument settings

All our observations were obtained with the FORS2 instrument (Appenzeller et al. 1998) of the ESO VLT. A multipurpose instrument capable of imaging and low-resolution spectroscopy, FORS2 is equipped with polarimetric optics. Polarimetric optics consist of a retarder waveplate (either a λ∕2 retarder waveplate for linear polarization measurements, or a λ∕4 retarder waveplate for circular polarization measurements) followed by a Wollaston prism that splits the radiation into two beams linearly polarized in perpendicular directions separated by 22′′. In front of the retarder waveplate, a Wollaston mask with nine 22′′ slitlets prevents thesuperposition of the field of view with the beam split by the Wollaston prism, following the optical scheme of Scarrott et al. (1983).

The FORS2 instrument allows us to directly measure the (wavelength dependent) quantities P_Q = Q∕I and Ū ∕I, and thus the fractional polarization P(λ) (or degree of polarization, referred to in this work as “polarization”, in percent) defined as $$P=\sqrt{(P_Q{}^2+P_U{}^2)}.$$(1)

The angle ofpolarization (ϕ) may be obtained from $$\tan (2\cdot \varphi )=U/Q.$$(2)

As our slit direction (and therefore the FORS2 instrument rotation) was consistently oriented in the east-west direction, we rotated the measured quantities Q∕I and U∕I by 90° into the instrument reference system to yield an equatorial reference system pointing to the north celestial pole. In the following, we calculate and discuss both quantities P and ϕ.

Most of our observations were obtained with grism 300V, which spans the optical range between 3600 and 9200 Å. An order separating filter may be inserted in the optical beam to cut the radiation shortward 4200 Å and prevents second-order contamination at λ≳6500Å. We obtained observations with and without an order separating filter and confirm that polarization spectra are not affected by contamination from second order (Patat et al. 2010). We also used grism 600 I that covers the spectral range 6700–9300 Å. We adopted a 2′′ slit width for a spectral resolution of 220 and 750 with grism 300V and grism 600I, respectively.

The performance and general calibration techniques of the polarimetric modes of the FORS2 instrument are well known and documented. In particular, the variation of the instrument polarization over the field of view, the chromatism of the retarder waveplate and the cross-talk between different Stokes parameters has been characterized in detail (see, e.g., Bagnulo et al. 2009, 2011). The systematic effects are typically on the order of 0.1% or less in the field center, which is practically negligible as contribution to systematic errors in our context.

2.2 Data acquisition

The FORS2 detector consists of two chips separated by a 4′′ gap. As a general rule, we aimed at pointing to the Moon with chip 1 containing Earthshine and chip 2 background sky, hence with instrument position angle on sky (relative to the north celestial pole) equal to 90° when pointing to the waxing moon or 270° when pointing to the waning moon. Acquisition consisted of a first pointing to the center of the Moon followed by a 15 arcmin offset either to the east or west, and a final tuning using more direct imaging at the lunar limb and through-slit images. An example of a typical acquisition image is shown in Fig. 1, on the left, obtained immediately before the dataset corresponding to ID F.6 (Table A.1). The lunar limb is clearly seen, and five 22′′ long slitlets on chip 1 are superposed to scale on the image to indicate their actual position because through-slit images do not allow us to reliably identify the actual position on the Moon. We always tried to position the limb in the gap between both chips to avoid observing the immediate vicinity of the limb. The four slitlets of chip 1 that cover the background are not shown in the image.

Science observations were generally obtained with the retarder waveplate set at all position angles 0, 22.5, 45, …, 337.5°, in order to apply the so-called “beam swapping technique” (see, e.g., Bagnulo et al. 2009). Basically, at each position angle of the retarder waveplate, the light flux is measured in perpendicular polarization states. The measurement is repeated after a 45° rotation of the λ∕2 retarder waveplate, thus with swapped polarization states in the two beams split by the Wollaston prism. The combination of the fluxes measured in the various beams allows us to remove most of the instrumental effects and calculate the reduced Stokes parameters P_X = X∕I (where X represents the various Stokes parameters Q, U, and I is the unpolarized intensity) with an uncertainty ideally given by ${(\text{S/N})}^{-1}$. Practically, the actual uncertainty of our Earthshine measurements is not determined by the photon noise, but is dominated by systematic instrumental effects and by the fact that the reflection properties of the lunar surface are only approximately known (see Sect. 4.2.2). Uncertainties were also statistically checked through inspections of the so-called null profiles, as discussed by Bagnulo et al. (2009).

2.3 Observations and Earth viewing geometry

In Table A.1 we list all parameters that characterize a given cycle of our polarimetric observations of Earthshine: a unique identifier (ID), the date and time when the observation cycle started, the airmass at this time, and the instrumental setup with the grism used. We also list the exposure times of individual exposures. The total exposure time is obtained by multiplying this number by the total number of exposures at individual settings of the retarder plate, which is also specified in the table. The phase angle α is defined as the angle between Sun – Earth – Moon. The viewing sceneries of Earth refer either to observations of the region to the west of the place of the observations in Chile, i.e., the Pacific (“P”, observed after sunset), or to observations of the region to the east, i.e., the Atlantic (“A”, observed before sunrise). As Earth rotates around its axis, its phase angle changes slowly. Also global weather patterns are changing slowly, but continuously. We expect these patterns to modulate the observed polarization steadily, butnot abruptly (see Sect. 4.6). We indicate observations that belong to the same sequence of observationsby the same capital letter for their ID in Tables A.1 and 1, while distinguishing individual observations with a numerical suffix. Observations a few weeks later, however, are likely influenced by completely different cloud coverage maps, which may have a strong effect on the overall polarization signal. Distinct observing sequences are therefore separated by larger horizontal spaces in the table and have different letters in their IDs.

A few examples of the appearance of Earth during the observations for representative observing epochs are shown in Fig. 2 as true color RGB composite images. The images were generated by MYSTIC radiative transfer simulations at wavelengths of 645 nm (red), 555 nm (green), and 469 nm (blue) as described in Emde et al. (2017). As input we used three-dimensional cloud field data from the ECMWF model closest to the date and time of the observations and land surface albedo data derived from MODIS (Schaaf et al. 2002). The ocean surface is simulated using a polarized bidirectional reflectance distribution function (Tsang et al. 1985; Mishchenko & Travis 1997). The aspect of the surface of Earth is seen from the lunar center.

Fig. 1 Left panel: acquisition image of the lunar limb that contains Earthshine observed on 2013-02-04UT08:10:50. Right panel: map of the apparent albedo from Velikodsky et al. (2011) of the same region that has been scaled and rotated. Corresponding FORS2 slitlets are superposed on both images.

3 Data reduction

The FORS pipeline (Izzo et al. 2010) handles spectropolarimeric data of point sources, but is not designed for spectropolarimetry of extended objects. Therefore, various dedicated procedures were developed and applied to our Earthshine data.

3.1 Frame preprocessing

All science frames were preprocessed with the ESO FORS pipeline (vers. 4.8.7.) to remove bias and perform a two-dimensional wavelength calibration. Among the various final pipeline products we used only the frames mapped in wavelength and corrected for field distorsion, but not flat fielded.

3.2 Flat fielding

In contrast to what is normally expected for polarimetric data obtained with the beam-swapping technique, flat fielding is important in the reduction of spectropolarimetry of Earthshine. The reason is that Moonshine background is not constant but must be interpolated in chip 2 and then (linearly) extrapolated to chip1 to be correctly substracted from Earthshine observations (see Sect. 3.3). Flat-field spectropolarimetric images may be obtained using either the calibration unit or twilight sky. Twilight sky data should be more suitable than data obtained with the internal calibration unit because twilight sky images follow the same optical path as science data. Therefore, these images should better represent large-scale spatial gradients in the system response. Ideally, flat-field data should be obtained through a continuously rotating retarder waveplate to measure a totally unpolarized image. This option is not implemented by the instrument software; therefore, our master flat fields were obtained adding up images obtained at various retarder waveplate position angles. Data were smoothed in wavelength, and the master flatfield was normalized along the direction perpendicular to the dispersion, by dividing the entire image by a one-dimensional spectrum obtained as average of its central 50 raws. We note that with a suitable combination of the frames obtained at various positions of the retarder waveplate (e.g., all positions 0°, 22.5°, 45°, …, 337.5°) we could in theory obtain a totally unpolarized master flat field, but in practice this did not work because the sky polarization changes quickly during twilight. This did not appear to be a major problem because the total flux needs to be calibrated to correct for field distortion; this flux is obtained by adding the two beams of the master flat field.

Sky flat fields were available for datasets with IDs D.x and observations later than December 2012, whereas screen flat fields were available for all datasets. We used sky flat fields whenever available. In order to minimize the effect of using different type of flat fields, we constructed stacked and smoothed sky flat fields from those available, and applied the correction to all screen flat fields. We compared the extracted spectra for pure sky- and screen-corrected flat-field calibration and did not find any significant differences in the extracted spectra. As we are not aware of changes of instrument configuration happening during April 2011 to December 2012 that could affect the optical beam, we are confident that this procedure ensures an optimal flat-field correction even for those datasets without corresponding twilight sky calibration data available.

3.3 Background subtraction

Background subtraction is a crucial step in the entire data reduction process. Its accuracy relies on the assumption that (1) field distortions have been correctly removed by two-dimensional wavelength calibration and flatfielding, (2) that Moonshine intensity decreases linearly with the distance from the terminator, and (3) that Moonshine polarization is constant within the 7 arcmin range of the PMOS slit. For each wavelength bin, at the spatial position y, the Earthshine polarization spectrum is then obtained as $$P_X{}^{(\text{ES})}(\lambda ,y)=\frac{P_X{}^{(\text{tot})}(\lambda ,y)F^{(\text{tot})}(\lambda ,y)-P_X{}^{(\text{bkg})}(\lambda )F^{(\text{bkg})}(\lambda ,y)}{F^{(\text{tot})}(\lambda ,y)-F^{(\text{bkg})}(\lambda ,y)},$$(3)

where $P_X{}^{(\text{tot})}(\lambda ,y)$ is the polarization measured at wavelength λ and positiony on the CCD (e.g., measured in arsecs from the edge of the CCD that is outside the lunar disk); F^(tot)(λ, y) is the total measured intensity; $P_X{}^{(\text{bkg})}(\lambda )$ is the polarization of the background at wavelength λ and average over the y range in chip 2; and F^(bkg)(λ, y) is obtained as linear extrapolation of the flux measured in chip 2. The background subtraction method also works in absorption band regions (like those around the O₂ bands), although overall fewer photons are available as compared to the continuum and a corresponding reduction of the S/N achieved.

In the ideal case, $P_X{}^{(\text{ES})}(\lambda ,y)$ should actually be independent of y (i.e., along the slit direction), but this is rarely the case. This is either because instrument distortions were not perfectly corrected, or the background does not change exactly linearly with distance. However, the largest contribution is due to lunar albedo changes along the slit (see Fig. 1), and its effect on the lunar polarization efficiency, as discussed in Sect. 4.2.2. Indeed, a clear anticorrelation between intensity and fractional polarization exists along each slit. Therefore, we always average the degree of polarization along the slit.

Fig. 2 Aspect of the Earth for eight representative observing epochs. True color RGB composite images of Earth as seen from the Moon, simulated using the MYSTIC radiative transfer model with cloud data from the ECMWF forecast model and MODISsurface albedo. The epochs correspond to the phase angles during observations.

4 Results

4.1 Spectral characteristics

We extracted polarization spectra for all observations listed in Table A.1 consistently applying the method described in Sect. 3. Figure 3 shows the polarization spectra separated in four groups: viewing sceneries pacific and atlantic, for viewing phase angles α above and below 90°. Following Sect. 2.2, the fractional polarization P is shown versus wavelength, ranging from 4200–9200 Å. The formal error of the polarization spectra mainly depends on the net signal obtained (dark moon minus background extrapolation). Summing up the signal along the slit direction, and for all retarder settings, we typically deal with total counts above 10⁶ per wavelength bin, allowing a formal statistical accuracy of below 1o for all polarization spectra discussed in this work. Spectra obtained with the 600I sometimes exhibit residual “fringing patters” redward of 8000 Å possibly caused by variability of telluric sky lines. These patterns could not be removed by our flat-fielding procedure.

Most polarization spectra are smoothly decreasing with increasing wavelength (except for molecular bands, see below). The spectra for datasets with IDs B.x appear to have a broad “bump” in the blue spectral range around 5000 Å. We verified that the flat-field procedure did not introduce artifacts, but we cannot exclude some systematic errors for these datasets in this spectral regime.

The different shapes of the individual polarization spectra as shown in Fig. 3 hamper their systematic comparison. But the main shape factors that distinguish the spectra can easily be identified are (a) the absolute value of the polarization at a certain wavelength; (b) the slope of the polarization across certain wavelength bands, in particular in the blue and the red; and (c) the shape and strength of the O₂ -A band (7600 Å) and to a lesser extent the much weaker O₂-B band (6700 Å).

Certain simplifications of the description of polarization spectra have been introduced before. Bagnulo et al. (2015a) suggested normalizing different polarization spectra at a certain wavelength (in their case at 5500 Å) to allow for the comparison of spectra of different asteroids observed at different phase angles. In their case, such a normalization removes the phase angle dependence of their observations and hence enables the grouping of the observations in different object classes that are otherwise hard to find.

For our purpose, we simplify the description of the detailed spectral shape of individual observations by introducing mean polarization values for different spectral bandpasses. The midpoints of the spectral bandpasses are derived from the usual astronomical (Johnson) filter system, which is centered around 4450 Å(B), 5510 Å(V), 6580 Å(R), and 8060 Å(I). A single observation with grism 300V allows us to extract the four mean polarization values and allows for an accurate differential measurement of the polarization values in these bands simultaneously. In practice, and allowing for inclusion of data from other sources, we defined 200–600 Å wide bandpasses across which we average the measured degree of polarization. The bandpasses chosen are P_B: 4350–4550 Å, P_V: 5450–5650 Å, P_R: 6450–6650 Å, and P_I: 8050–8650 Å. The choice of a wide bandpass for P_I is motivated to extend the usual I-bandpass to enable future comparison with data obtained with the POLDER satellite instrument (Deschamps et al. 1994), which is the only instrument that provides polarization measurement from space. The POLDER instrument measures polarization only in its reddest channel (centered around 865 nm). The choice of "spectral bands" allows a rather accurate determination of the mean polarization value within the band because the spectral slope across these relatively narrow bands is almost constant, and the average value essentially eliminates the residual, statistical variations within the band. The errors in the polarization values are thus given by their standard deviation measured in their passbands, and are typically very small (< 0.1o).

Table A.1 lists the values of degree of polarization and their statistical errors for all spectra observed. The spectral coverage of the 600I grism only allows determination of P_I. For comparison, Table A.1 includes values determined from various sources: low-resolution spectra of Takahashi et al. (2013) cover the wavelength range of 4500–8500 Å and allow us to determine the polarization values in all four bandpasses defined above. Their observations consistently contained the African and Asian continents in Earthshine. The values reported in Bazzon et al. (2013) refer to measurements of polarization in specific bandpasses, focusing on two regions on the Moon having distinct albedos (lunar highlands and mare). They used the standard Bessell B, V, R, I filter set. Unlike our definition of P_I that spans a region from 8050-8650 Å, their I is centered around 8000 Å. We estimate the value of polarization P_I, assuming linear extrapolation from their V, R, and I colors to a wavelength of 8550 Å. This should reduce bias in this band to compare with our data.

Within the same bandpasses defined above, we can calculate the angle of polarization ϕ_B, ϕ_V, ϕ_R, and ϕ_I. For a Rayleigh scattering atmosphere, this angle should coincide with the normal of the scattering plane on Earth, which is defined by the Sun, observer (on Earth), and the Moon, with a reference direction of the north celestial pole. We calculate this angle Φ using Eq. (5) from Bagnulo et al. (2006) for each geometrical configuration of the Sun-Earth-Moon system, and list it in Table A.1. As expected, the values of ϕ_B, ϕ_V, ϕ_R, and ϕ_I have at most only a weak dependance on the wavelength, and coincide quite well with Φ as calculated. A similar behavior was found in Takahashi et al. (2013) and serves as an important sanity check of the quantities P_Q (λ) and P_U (λ) measured.

Fig. 3 Spectra of the fractional polarization of Earthshine as observed with FORS. For better visibility, the spectra are grouped as follows. Left panels: spectra of the Pacific in the Earthshine (to the west of the observing site), right panels: spectra of the Atlantic (to the east). Upper panels: spectra taken at phase angles α < 90°, lower panels: those at α > 90°. Individual observations are identified by colors specified in the legend, and each spectrum is identified with its observation ID and its actual phase angle in Table A.1.

Fig. 4 Measured fractional polarization of Earthshine as a function of the phase angle (Sun – Earth – Moon). The different symbols refer to different viewing geometries during the observations: the Pacific or Atlantic ocean. The bandpasses used correspond approximately to the traditional B, V, R, and I bandpasses (see text). The lines connect at least three independent observation cycles that belong to a consecutive sequence of observations within a given observing run (see Table A.1).

4.2 Polarization phase curves

Astronomical objects (with or without any atmosphere) exhibit a characteristic variation of their polarization with phase angle (see, e.g., Kolokolova et al. 2015). Ever since Dollfus (1957) it is known that the observed polarization of Earthshine follows a characteristic phase curve. Figure 4 shows our measured polarization values at phase angles from 50° to 140° (see Table A.1). The four passbands (B, V, R, and I) have been indicated by blue, violet, green, and red lines, respectively. The two different types of symbols distinguish the viewing geometry of the Earth during the observations: left half circles are observations of the Atlantic (east to the observational site in Chile), while right half circles are of the Pacific (west to the observational site). This distinction is important as it indicates both the differentglobal sceneries and the fact that a different lunar limb region was observed at the time. One set of polarization measurements in the four bands at a given phase angle and with a given viewing geometry (Pacific or Atlantic) represents a single Earthshine polarization spectrum.

In order to simplify the interpretation of the cluster of points in Fig. 4, we connected at least three independent observation cycles that belong to a continuous sequence of observations in a given band and observing epoch, by lines using a second-order polynomial for a least-squares fit procedure. The lines highlight the overall shape and grouping of polarization values. Individual observations that belong to one sequence of observationswithin one observing run lie in general close to the connecting line. The polarization generally follows a smooth curve as a function of the phase angle. The polarization reaches its highest value roughly around α = 80°–100° (i.e., around quadrature) and decreases toward lower and higher phase angles within our observational range. There is a tendency for longer wavelengths to have the polarization maximum at larger phases angles.

There isa significant dispersion of the polarization when comparing distinct observation epochs in all bands. The variation of the polarization values appears to be largest around quadrature. For example, observations of the Pacific side around quadrature exhibit maximum polarization values in the blue of more than 16% (B.2, B.3, and G.1), while observations at a similar geometry of the Atlantic side are about 3% lower (F.4). In general, the polarization appears to be lower for observations on the Atlantic side than on the Pacific side. This was already noted by Sterzik et al. (2012) for datasets A.3 and B.4 and interpreted in terms of different cloud fractions at the time of the observations; higher polarization corresponds to a lower cloud coverage fraction. A further important difference between the Pacific and Atlantic observations is the highly polarized sun glint on the ocean, which is mostly visible on the Pacific side and only partly on the Atlantic side (Emde et al. 2017).

4.2.1 Comparison with other observations of Earthshine polarization

How do the measurements relate to observations published before? For comparison, Fig. 5 shows values from the literature as listed in Table A.1. Our Earthshine measurement can most directly be compared to those of Takahashi et al. (2013) and Bazzon et al. (2013). Bazzon et al. (2013) reported polarization values in certain passbands. From the spectra of Takahashi et al. (2013) we derived the passband values with the same procedure as for our own spectra. We also connected their data by lines following least-squares fits, as all observations fall within a few days. In general, their polarization values follow the same trends as ours. The data that belong to one continuous set of observations show little scatter around the fitted lines. However, compared amongst each other and with our data, the scatter in the Earthshine polarization data of different authors obtained at different epochs but at similar phase angles reaches around 2%.

An important source of scatter of the polarization data was clearly identified in Bazzon et al. (2013), who measured Earthshine polarization in two different regions on the lunar surface (on Highlands and in Mare), quasi-simultaneously. The different surface types covering these regions (their albedos are 0.21 and 0.11, respectively) appear to yield significantly different lunar depolarization factors, and hence a relative difference of the measured polarization of the Earthshine of up to 30% in the blue (B) spectral band. In the next section, we investigate the effect of lunar depolarization on our data.

4.2.2 Correction for lunar depolarization

Earthshine is reflected by the lunar surface. Depending on the composition and structure of the local lunar surface, the reflection changes the state of polarization of the Earthshine. The lunar depolarization factor or polarization efficiency ϵ is defined as $$ϵ(\lambda )=P^{\text{out}}(\lambda )/P^{\text{in}}(\lambda ),$$(4)

where Pⁱⁿ is the polarization of the light that is incident on the moon and P^out the polarization of the reflected light. We note that to measure ϵ, one could illuminate lunar (analog) surface samples with light with a known polarization state and measure the polarization of the reflected light for a wide range of illumination and viewing geometries and wavelengths. As far as we know, such measurements have not been done yet.

Bazzon et al. (2013) introduced a method to correct Earthshine polarimetry using the knowledge of the lunar albedo at the location where the Earthshine is observed, assuming that ϵ depends on the wavelength and the lunar surface albedo. In addition, it is assumed that the phase angle dependence of ϵ is negligible because the angle between the lit part of Earth and the observer, as seen from the moon, is small and always around 1° ± 0.5°. We derived the spectral dependence of the polarization efficiencies on albedo through an analysis of lunar samples by Hapke et al. (1993). Finally, we obtained lunar albedos from an extrapolation of absolute lunar albedos maps by Velikodsky et al. (2011) to backscatter angles of 1°. The polarization efficiency, log ϵ, as a function of the lunar albedo at 603 nm, loga₆₀₃, and the wavelength, logλ, as derived by Bazzon et al. (2013) is then $$\log ϵ(\lambda ,a_{603})=-0.61\log a_{603}-0.291\log \lambda [\mu \text{m}]-\mathrm{0.955.}$$(5)

We used the same approach to derive the polarization efficiency for our observations. The acquisition images obtained immediately before the spectropolarimetric measurements were used to identify their location on the lunar surface, and we matched these images with the extrapolated 1-deg albedo maps by Velikodsky et al. (2011). This can be seen in Fig. 1 where we scaled and rotated the albedo map to best match the orientation of the acquisition image. The determination of the rotation angle α_Moon with this procedure is affected by field distortion, the finite spatial resolution of the images, and the sometimes low contrast of the acquisition images. We estimated the accuracy of the rotation is limited to ± 1°, while the accuracy to locate and extract the correct albedo values at the position of the slitlets is accurate to no more than ± 15–20′′ on the moon. We therefore sampled the albedos a₆₀₃ scanning theslit mask ± 5 pixels (corresponding to ±18′′ on the moon) around its suspected center position in and perpendicular to the slit direction. We then derived a mean albedo ${\overline{a}}_{603}$ for all 11 × 11 raster slit images and estimated the error of the albedo measurement using the standard deviation within these samples. In Table 1 we list the rotation angle α_Moon (counted counterclockwise from the north) and the mean albedo ${\overline{a}}_{603}$ and its standard deviation derived from these positions for each observation.

We plotted the lunar depolarization corrected polarizations for all four passbands in Fig. 6. As in Fig. 4, different colors indicate different bandpasses, and the two different symbols distinguish the observations of the Pacific (P) and the Atlantic (A) side of the Earth. The errors derived from applying minimum and maximum polarization efficiencies due to the scatter in the lunar albedos across the region where the Earthshine is measured are shown with vertical bars. For each passband and orientation of the Earth (“P” or “A”), the data appears to be scattered around a mean phase curve. To fit these phase curves, we used a modified Rayleigh function of the form (Korokhin & Velikodsky 2005) $$P(\alpha )=\frac{{({\sin }^2(\alpha -\Delta \alpha ))}^W}{1+{\cos }^2(\alpha -\Delta \alpha )+\text{dePol}}.$$(6)

Parameters W, Δα, and dePol characterize each curve with respect to its width and phase shift with respect to α = 90°. Parameter dePol relates to the maximum polarization. In Fig. 6, we show the best fits with the 1σ confidence intervals indicated with gray bands. The fit parameters are listed in Table 2. The upper and lower errors are obtainedby fitting the data with the upper and lower polarization values obtained by assuming different polarization efficiencies ϵ due to systematic albedo variations as listed in Table 1. The datasets over the Pacific and the Atlantic are distinguishable by different fit parameters that are statistically significant at the 1–3σ level. It is interesting to note that the differences between the sets are more significant in the redder (R and I) spectral passbands. Another interesting feature in the curves is the apparent crossing of the polarization curves for the different bands near 120°: planet Earth becomes “white” in polarization.

The scatter of the observed polarization that is visible in Fig. 6 is caused both by intrinsic scatter resulting from the relative uncertainties of the lunar albedos and the associated uncertainty in the lunar polarization efficiency, and by the sensitivity of the polarization to the atmospheric and surface properties of the Earth as seen from the Moon at the time of the observations. In particular data points outside the error regime of lunar depolarization are likely caused by differences in the atmospheric and surface patterns, as they influence the Earth’s polarization significantly. A comparison with dedicated models should help to understand the causes of these deviations.

Fig. 5 Polarization of Earthshine as a function of phase angle α. The different symbols refer to different observations, including observations by Takahashi et al. (2013) and Bazzon et al. (2013). The bandpasses correspond approximately to B, V, R, and I bandpasses, as described in the text.

Fig. 6 Polarization of Earth after correction for lunar depolarization effects.

4.2.3 Comparison with models

In this section we compare our measurements with idealized models of polarization spectra of Earth-like planets by Stam (2008). The results of these numerical models are very useful for a qualitative comparison with our data. They are available in tabulated form with sufficient spectral resolution for the full phase curve and for a variety of parameters that characterize the surface and atmosphericproperties such as the fraction of free ocean and land surfaces (including vegetated areas), as well as cloud cover. The tabulated data are limited in that they apply to horizontally homogeneous model planets, that the reflection by land surfaces is only Lambertian (i.e., isotropic and depolarizing), and while the reflection by ocean surfaces is described by Fresnel reflection, thus including polarization, there are no waves on the ocean. The glint due to the reflection of the direct, i.e., nonscattered, sunlight on the water is thus described by a delta function, not by an extended region, as it would be on a wind–ruffled surface (Zugger et al. 2010, 2011). The tabulated data further include a single type of (liquid water) cloud only that has an optical thickness of ten. Despite the limitations, these tabulated data can serve as a first approximation to explain the most prominent observational features of the Earthshine and may allow us to identify the main physical mechanisms that impact the Earthshine spectra. These data have already been used in Sterzik et al. (2012). By taking weighted sums of the tabulated data, we mimic horizontally inhomogeneous planets (Stam 2008).

In Fig. 7 we compare Earth’s polarization values, which have been individually corrected by polarization efficiencies as described above, with three representative models of Stam (2008). We consider cloud coverage fractions of 30, 40, and 50%, and the remainder of the disk is covered by cloud-free ocean. As expected, in every bandpass, the polarization decreases with increasing cloud coverage because at most phase angles the polarization from the clouds is lower than that from the atmospheric gas above the dark ocean. The polarization decreases with increasing wavelength because the scattering by the gas decreases with increasing wavelength. The scattering by the clouds thus becomes more prominent with increasing wavelength. In the red bandpass, P_Q changes sign at phase angles largerthan 90° with 40 and 50% cloud fraction³, because at those phase angles the polarized flux reflected by the clouds has a direction parallel to the reference plane and the angle of polarization ϕ defined in Eq. (2) changes by 90° (see Karalidi et al. 2012a, for models at a large range of cloud coverages and phase angles). We note that our observationsdo not include the phase angle range in which the rainbow due to the scattering of light in the spherical water cloud droplets is expected: this peak in polarization would appear around α = 40°, where the increase in polarization toward the peak starts below α = 50° (Karalidi et al. 2012a), just the smallest phase angle in our observations. The lines for the lower cloud fraction and for the higher cloud fraction bracket the observations for the B band and for most observations in the V band. The lines for the R band and in particular the I bands, however, are too low to bracket the observations. We conclude that models of the Earth with a pure ocean surface and liquid water clouds bracket the observed polarization phase curves in the blue spectral region. Differences in the observed cloud coverage fractions and in particular cloud types (optical thickness, cloud particle size, and possibly thermodynamic phase) for the different observing epochs likely cause the different polarization fractions observed during the course of one observation epoch covering at most a few days. The red R and in particular I spectral bands remain a challenge to fit with a weighted sum of horizontally homogeneous cloudy and ocean-covered Earth polarization models.

The relatively flat continuum of the polarization spectra for red wavelengths as compared to Stam (2008) has already been noted in Sterzik et al. (2012). This has triggered the development and application of more sophisticated models, in particular for the treatment of clouds. A dedicated model was explained and presented in detail in Emde et al. (2017). Their simulation was actually designed to explain observations ID A.3 and ID B.4. In order to compare the simulation result with measurements, we derived polarization values directly from their spectra and list these in Table 1. As can be seen in Fig. 7, their polarization values tend to fit the observations in all four B, V, R, and I bands.

Fig. 7 Comparison of the fractional polarization of Earth in dependency of phase angle (Sun – Earth – Moon) with models in B, V, R, and I bandpasses. The error bars indicate the effects of uncertain lunar albedo on polarization. The phase curves for three representative models of Stam (2008) are overplotted with different line styles. Values extracted from the simulated spectra of Emde et al. (2017) for observation IDs A.3 and B.4 are indicated by stars.

4.3 Polarization color ratios

The comparison of Earthshine polarization with models of Earth polarization is hampered by the uncertainties induced by the rather inaccurate knowledge of the lunar (de-)polarization factor. The method of Bazzon et al. (2013) applied in Sect. 4.2.2 is only an approximation. The rather uncertain determination of the exact slit position on the Moon introduces an additional error budget that we estimated above. But in particular the rather unknown scattering parameters of (different) lunar soil across the Moon may introduce even more systematic errors, and Eq. (5) holds only approximately. The absolute polarization derived for Earth from Earthshine may thus be uncertain by a few percent with relative errors as large as 10%.

However, Eq. (5) suggests only a rather weak dependence of the lunar polarization efficiency ϵ on λ. A factor ten difference in wavelength results in a factor of less than two difference in ϵ, and the relative error introduced when considering polarization ratios such as P_B ∕P_V (or P_R ∕P_I) is expected to be only on the order of 4–5%. Thus, even if Eq. (5) may just be an approximation, polarization ratios, in particular pertaining to adjacent wavelength bands, should be rather robust quantities to compare with models. Such differential quantities could thus serve as generic observables of Earthshine that allow a more reliable comparison with models because they are less sensitive to the lunar depolarization factor and its uncertainties.

In Fig. 8, we show the polarization ratios P_B ∕P_V and P_R ∕P_I for our Earthshine polarization data as functions of phase angle α, using different (filled) symbols for the observations covering the Pacific side and Atlantic side. Errors on these ratios are within a few percent, and thus on the order of the symbol sizes. The polarization ratios for the data that are corrected for the lunar depolarization are indeed similar to the data that are uncorrected, thus confirming our choice for using polarization ratios. Figure 8 also shows model computations from Stam (2008) and Emde et al. (2017). Spectra with a relative shallow decrease of polarization with increasing wavelength have relatively small polarization color ratios, while flat spectra have P_B ∕P_V = 1 and P_R ∕P_I = 1. For both color ratios P_B∕P_V and P_R ∕P_I, there is a general tendency for an anticorrelation with phase angle α: the smaller α, the steeper the polarization spectra. However, for α > 100° all spectra appear to be flatter, and a few spectra from the “P” sample even exhibit slightly increasing slopes in the red spectral range.

It is interesting to compare the Earthshine (resp. polarization efficiency corrected) polarization color ratios with those for Moonshine. As explained above, each Earthshine observation is associated with a Moonshine observation, which is captured in the part of the detector that sees a region of sky adjacent to the lunar limb and that is used for background subtraction. The Moonshine consists of scattered moonlight and can be used to extract polarization spectra in the same way as we did for the Earthshine data. In Fig. 8 we included the Moonshine polarization color ratios, indicated with nonfilled symbols. Over a wide range of phase angles, the slopes of Moonshine spectra, and thus their polarization color ratios, are distinct from those of the Earthshine/Earth spectra. While the absolute value of the local lunar polarization depends significantly on the local lunar albedo, there is no significant difference in the color ratios of the Moonshine originating from the west- or the eastside of the Moon. For the phase angle range considered in this paper, the polarization color ratios of the Moonshine appear to be rather independent of phase angle. This behavior has been noted before by, for example, Gehrels et al. (1964) and Shkuratov & Opanasenko (1992).

The polarization color ratios of the Earthshine can also be compared to model simulations. As before, we use a representative set of Earth polarization models from Stam (2008). Conceptually, the simplest case is a model composed of a Rayleigh-scattering atmosphere with an ocean surface below. The ocean albedo is zero for all wavelengths, and while at the shortest wavelengths the atmosphere is optically thick enough for multiple scattering that lowers the polarization, with increasing wavelength, the polarization increases to its single scattering value of nearly 0.9, and becomes largely wavelength independent; the Fresnel reflection lowers the polarization slightly when compared to a black Lambertian surface. This mechanism is largely independent of α, thus the expected polarization color ratios are close to one for all phase angles, as can be seen in Fig. 8. A Lambertian reflecting surface with an albedo of 0.2 steepens the polarization spectra considerably, as can be seen from the polarization ratios; as with increasing wavelength, more light reaches the surface and adds more unpolarized light to the Earthshine (see Stam 2008). The anticorrelation of the measured polarization color ratios with the phase angle is actually approximated by such a model for the blue color ratio P_B ∕P_V, but in particular for the red color ratio P_R∕P_I the model ismuch too steep, as seen in Fig. 8.

The various dependencies of the polarization color ratios on wavelength and phase are even more complicated when clouds are considered. It is well known that clouds with their rich and complex macro- and microphysical properties have manifold impact on polarization spectra. With a 20% cloud coverage fraction in the models of Stam (2008; combined with 80% clear ocean surface), the polarization color ratios increase to values qualitatively compatible with our observations for the blue ratio P_B ∕P_V, but they are off-scale for the red ratios P_R∕P_I. Increasing the cloud coverage to a more realistic 50% (with 50% clear ocean surface) steepens the polarization spectra even more and appearsto become incompatible with both the observed blue and red polarization ratios. For comparison, we also include the polarization color ratios derived from the simulations of Emde et al. (2017) to fit datasets ID A.3 and B.4. As expected, their more detailed cloud and surface properties, as derived from remote-sensing data at the time of the observations, appear to flatten the spectra and lower the polarization ratios and agree more quantitatively with the observations, but do not match the observationswithin the errors.

Fig. 8 Polarization color ratios PB∕PV (left panel) and PR∕PI (right panel) as functions of α. Earthshine data are indicated by filled symbols. Blue left half-moon symbols indicate the Pacific side and green right half-moon symbols the Atlantic side. Values derived from Moonshine spectra are indicated by nonfilled (open) symbols for the respective sides of the Moon from which the moonlight is scattered. Models of Stam (2008) are overplotted with different linestyles. Values derived from model spectra of Emde et al. (2017) are indicated by stars.

4.4 Polarization vegetation index

The previous section concentrated on the shape of the continuum polarization spectra. In order to study gaseous absorption bands and other variations in the spectra, we normalize the spectra by subtracting a continuum, following the same procedure described in Sterzik et al. (2012): a fourth-order polynomial is fitted to a spectral range between 5300 and 8900 Å. Spectral regions that contain band and/or variations above 1.5 sigma above the mean values are excluded from the fit. The resulting fit is subtracted from the original spectrum, and the residuals then represent the spectral variations on the spectrum. Figure 9 shows a representative set of observations in the interesting spectral region between 6400 and 8000 Å. This region contains the molecular absorption bands O₂ -B (around 6900 Å), H₂O (around 7200 Å), and O₂-A (around 7600 Å), as well as the VRE, the spectral region where the reflectivity of Earth’s vegetation sharply increases (Tinetti et al. 2006a).

The VRE shows up in the polarization spectra mainly because of the added unpolarized flux that is reflected by the vegetation. In order to have a closer look at the VRE, we eliminated effects of the molecular absorption bands as much as possible. We therefore defined two bands, blueward and redward of the VRE (Tinetti et al. 2006a) in the polarization spectra: between 6750 and 6850 Å and between 7480 and 7780 Å, but avoiding the region affected by the O₂ -A band between 7580 and 7680 Å.

Fig. 9 Selection of normalized spectra between 6400 and 8000 Å. The horizontally dotted lines indicate the regions across which the vegetation index ΔPVI, associated with the VRE: blueward (between 6750 and 6850 Å; green dotted line) and redward (between 7480 and 7780 Å, but avoiding the region affected by the O2 -A between 7580 and 7680 Å; red dotted line). The difference between the green and red dotted lines defines the amount of ΔPVI. The derivation of EW(O2-A) is explained in the text.

4.4.1 Analysis of ΔPVI

We averaged the normalized polarization over all wavelengths across both regions, and call that difference the “differential polarization vegetation index” (ΔPVI). Errors of the ΔPVI are defined by the standard deviation in the bands. Negative values of ΔPVI indicate a lower continuum in the blue than in the red as a result of a sharp increase of the albedo of vegetated surfaces in the red.

Figure 10 shows all values of ΔPVI from the observations listed in Table 1 as functions of the phase angle α. As before, different symbols indicate different sceneries (blue, left-half circles indicate the Pacific and green, right-half circles the Atlantic).

As canbe seen in the Fig. 10, values of ΔPVI are scattered, and there is only a weak correlation with phase angle. But the values of ΔPVI for the Atlantic side tend to be more negative than those for the Pacific side. We quantify the correlation by a formal linear regression analysis of the parameters and their errors. The population is sparse and has outliers and errors. Therefore, we apply maximum-likelihood estimator techniques described in Kelly (2007). The 1σ confidence bands around the linear regression have been shaded in gray for the two subsamples. Both populations become distinct with increasing significance for larger phase angles. While for small phase angles both distributions overlap, the means of the two distribution are distinct by more than 2σ for phase angles around 110°. As outliers exist in both populations that are compatible with the mean values of the other population, the significance of the different linear regressions is lowered. However, a formal KS-test gives a probability of only 0.143% for the two populations “A” and “P” being drawn from the same underlying ΔPVI distribution. In this sense, both distributions are statistically different.

By design, ΔPVI is supposed to be sensitive to the amount of visible surface vegetation and caused by the steep albedo change across these wavelength bands. In order to compare the observations with models, we derived the ΔPVI parameter using polarization model spectra of Stam (2008); the models consider only one type of vegetation, i.e., deciduous forest. In order to avoid methodological biases, we apply the same procedure to remove the continuum and to derive ΔPVI as done for the observations. It is evident from the curves referring to the models in Fig. 10 that ΔPVI depends sensitively on the amount of visible vegetation at most phase angles. As an ensemble, the observations of the Atlantic side are bracketed by models containing up to 3–10% surface vegetation. Observations of the Pacific have a larger scatter, but models without vegetation better describe their distribution.

We include in Fig. 10 the values derived from the two simulations of Emde et al. (2017) as listed in Table 1. Their models contain 3 and 10% visible vegetation, and the value of ΔPVI = −4.49o corresponding to the high vegetation case for observation A.3 falls outside the range shown in the figure.

4.4.2 Relation to NDVI

The parameter ΔPVI is sensitive to the difference in polarization, which is supposed to be induced by the abrupt increase of the surface albedo of Earth’s vegetation in the near-infrared. The observed range of ΔPVI values for the “A” and “P” subsamples can be explained with models that contain different fractions of surface covered by vegetation. According to these models, already small fractions of vegetation enhance the signal. Our measurements appear to be sensitive to vegetation fractions larger than 3% and support the earlier claim in Sterzik et al. (2012) that “A” and “P” sceneries can be distinguished by their different amount of visible vegetation contained in Earthshine.

The presence of outliers in both (“A” and “P”) subsamples of the ΔPVI distribution likely washes out the statistically differences of the ensemble averages, in particular for smaller phase angles. As the amount of visible vegetation cover appears to sensitively change the ΔPVI parameter, individual observations are expected to depend sensitively on the actual presence of vegetation in a given underlying scenery.

We approximate the actual (or true) vegetation cover of individual observations by deriving the “normalized differential vegetation index” (NDVI). Based on the visible parts of Earth, and their MODIS surface albedo for the R and I bands, we calculate NDV I = (I_858nm − I_645nm)∕(I_858nm + I_645nm) to measure the proportion of each image pixel covered with vegetation (cf. Fig. 2) and derive an average vegetation cover of the total illuminated area for each observing epoch.

Figure 11 shows the relation of ΔPVI with the actual vegetation cover observed. In general, observations with a larger actual vegetation cover fraction have a lower value of ΔPVI. The quantities are anti-correlated, but outliers exist. It is interesting to note that even the “P” sample sometimes contains observed sceneries with a significant vegetation cover, and vice versa. ΔPVI appears to be sensitive to detect these cases.

We conclude the discussion of potential effects of vegetation on our data by noting that the series of observations (B.1 to B.5) contained a broad “bump” in their polarization spectra around 5000 Å (see Fig. 3). We cautioned earlier that we cannot 100% exclude artifacts caused by the flat-fielding procedure. But interestingly, the only observations of the Pacific ocean during northern summer are datasets B.1–B.5. We speculate that this bump could also be due to surface reflection: if the albedo increases toward 5000 Å (green vegetation), P decreases there. A possible absorber, which might only be present seasonally on the ocean, could be algae (chlorophyll B). However, the efficiency of this mechanism needs to be modeled before any further conclusion can be drawn.

Fig. 10 ΔPVI as a function of the phase angle with the symbols indicating the viewing sceneries (Pacific or Atlantic). The green and blue lines indicate the linear regression fit from all the observations of either scenery. Confidence bands of 1σ around the linear regression are shaded in gray for the two subsamples. While the Pacific side shows a very weak correlation with the phase angle, the population from the Atlantic exhibits decreasing ΔPVI with increasing phase angles. Both samples are statistically distinct according to a two-sample KS-test. Red lines indicate the ΔPVIs as derived from the models from Stam (2008). The star symbol indicates the ΔPVI as derived from the model for the B.4 from Emde et al. (2017); the value derived from the model for the A.3 dataset is −4.49, and falls outside the plot.

4.5 O₂-A band strength

The spectral resolution of our polarization spectra also allows us to determine the strength of the polarization signal in specific absorption band regions, most prominently in the O₂ -A band region (around 7600 Å). The polarization in this band is often higher than in the adjacent continuum; sometimes it is flat and it is rarely lower than in the continuum. The variability of the polarization in the band can readily be seen for the sample of spectra shown in Fig. 9. Various processes determine the polarization in an absorption band as compared to that in the continuum. The absorption of light by, for example, O₂ decreases the amount of multiple scattered light and there is usually a lower degree of polarization than the singly scattered light. This process yields a band polarization higher than that in the continuum. Absorption of light by O₂ also limits the amount of light that is scattered at low altitudes in the atmosphere. Indeed, the stronger the absorption, the lower the altitude from which the Earthshine originates. In a vertically inhomogeneous atmosphere, different types of particles at different altitudes can yield a polarization that varies across the band. In particular, the cloud top altitude influences the band strength. Finally, with increasing gaseous absorption, less light that has been reflected by the surface reaches the top of the atmosphere.Because a reflecting surface usually increases the amount of unpolarized flux, increasing absorption increases the polarization of the Earthshine. For more detailed explanations and sample computations for whole planet signals, see Fauchez et al. (2017).

To quantify the behavior of this feature, we introduced a quantity called equivalent width (EW), which is frequently used in stellar spectroscopy. In this case, it measures the area of polarization (P_λ) over wavelength (λ) integrated over a specific spectral region (from λ₀ to λ₁) normalized to its continuum value P_c as follows: $$\text{EW}(\lambda )={\displaystyle {\int }_{{\lambda }_0}^{{\lambda }_1}(1-P_{\lambda }/P_{\text{c}})}\text{\hspace{0.17em}d}\lambda .$$(7)

In practice, we numerically integrate P_λ across the passband λ₀ = 7580 Å to λ₁ = 7680 Å divided by the adjacent continuum level P_c. We determine P_c by a second-order polynomial fit of two 1000 Å wide regions red- and blueward of the band. Negative values of EW indicate a band in “emission”. We estimate the error in EW and also determine its values for both regions in the continuum. These values indicate the intrinsic error of the integration over a flat spectral region not affected by the spectral band. The high quality of the polarization spectra and the very good sampling of the O₂ -A band region lead to small errors for the determination of EW, typically less than 1%. The EW is largely independent of the spectral resolution used, and its values obtained by measurements with different grisms can be compared directly.

Figure 12 shows the values of EW(O₂-A) for all the observations listed in Table 1 as functions of the phase angle. As in Fig. 10, different symbols indicate different sceneries of the Earth. Evidently, there is a considerable scatter of the EWs, and only a marginal correlation within the range of phase angles for the Pacific and Atlantic sides. Again, we apply a robust least-squares fitting procedure to both samples as explained above, and show the two regression lines in Fig. 12 for the blue (Pacific) and green (Atlantic) sets. The gray areas correspond to a ± 1σ confidence interval around the optimal regression curves. Both regression curves are distinct by ≈ 2σ from each other, and offset over the full range of phase angles covered. As in the case of ΔPVI, outliers exist for both samples that are compatible with the other populations. A formal two-sample KS-test gives a very low probability (2.4 × 10⁻⁵) that both samples are drawn from the same underlying population.

In order to increase our understanding of the behavior of the EW(O₂-A), we extracted the same quantity for a set of models of Stam (2008). Although these model spectra are at a significantly lower spectral resolution than our observations, they allow us to derive the EW in the same way. Results for different model planets (all with the same cloud top altitude) are plotted in Fig. 12 with different line styles. Qualitatively, the models suggest an anticorrelation between the EW(O₂-A) and the phase angle with increasing cloud coverage fraction. Interestingly, cloud coverage fractions larger than 50% introduce significant variations in the EW(O₂-A) values.

We also determine the EW(O₂-A) for the two simulations of Emde et al. (2017; see Table 1). They show that the introduction of cloud layers at different altitudes and with different optical thicknesses and droplet sizes have significant effects on the appearance and strength of the O₂ -A band in polarization spectra. For a suitable choice of cloud parameters it may thus not be surprising that the EW(O₂-A) derived from their models correspond better to the observations.

Next we focus on observations of the O₂-A band obtained with higher spectral resolution. The appearance of this band in polarization depends on the fraction of (usually highly polarized) single scattered light to that of the (usually low polarized) multiple scattered light, which increases with increasing absorption. It also depends on the vertical distribution of scattering particles, including cloud particles, in the atmosphere, and it depends on the surface albedo. Examples of this band in polarization can be seen in Fig. 13 for those Earthshine spectra observed using grism 600I. Resolved fine structure in the band is clearly visible and some parts of the band show polarization to more than 15% above the continuum value of 20%, while others do not show enhanced polarization, or even slightly reduced polarization. Increasing the spectral resolution would further enhance the contrast in the band (see Stam et al. 1999, for example). More examples with high spectral resolution have been shown in Emde et al. (2017). The new observational contribution in this work, however, is the large variation of this line independent of phase and not directly correlated to “A” or “P” sceneries. Fauchez et al. (2017) have investigated the influence of the planet surface albedo, cloud optical thickness, altitude of the cloud deck, and the O₂ mixing ratio on the polarization in the O₂-A band. They have found that these parameters may not be easily discriminated in the higher optical depth regimes that can be probed with high spectral resolution. The cloud deck altitude and horizontal distribution across the region of the Earth that contributes strongest to the observed signal may be decisive factors for the appearance of the band, but more detailed simulations of Earth polarization measurements are needed to confirm this.

Fig. 11 ΔPVI as a function of the actual vegetation cover observed. The actual vegetation cover has been calculated from NDVI values for each image pixel. Different symbols indicate the viewing sceneries (Pacific or Atlantic). The black line shows a linear regression fit with all the observations. Confidence bands of 1σ round the linear regression are shaded in gray. The anticorrelation of ΔPVI with the actual vegetation cover is strong, although outliers exist.

Fig. 12 Equivalent width (EW) in the O2-A band as a function of phase angle α. The symbols indicate the different sceneries (Pacific and Atlantic). Full lines and gray areas indicate the linear regression and its errors. Both regression curves are distinct by ≈ 2σ from each other, and offset over the full range of phase angles covered. A very low probability found in a two-sample KS-test indicates that both samples are not drawn from the same underlying population. Red lines with different line styles refer to models of Stam (2008) and star symbols refer to the two models of Emde et al. (2017).

4.6 Short-term variability

In this section we try to disentangle the following three main effects that have an impact on P^E on a timescale of a few hours: (a) continuous phase angle changes due to the movement of Earth and Moon, (b) Earth rotates and makes different parts visible at different times, and (c) changes in cloudiness are introduced by changing weather patterns.

As described in Sect. 2.2 and in Table A.1, polarimetric spectra were usually acquired with 16 different settings of the retarder waveplates. Four settings are sufficient to reliably derive the Stokes parameters Q, U, and thus P, albeit with correspondingly less S/N. Therefore, we can increase the temporal resolution of the observations to about 15 min, which is the typical duration of an observing cycle with four retarder settings. The spectra are then equally processed as described in Sects. 2 and 3, and the same parameters (with their statistical and systematic errors) determined as describedin Sect. 4.

In particular the datasets observed on 2012-12-19, 2013-02-03, 2013-02-19, and 2013-02-22 allow a continuous monitoring of the polarization of Earth over about three hours, with the sampling time of 15 min as defined by a full observation cycle. We use these data to investigate the variation of P^E in the four passbands, and of ΔPVI and EW(O₂-A) on this timescale.

Figure 14 shows the variation of these parameters over time. The four panels correspond to four distinct observations dates. The colored half-moon symbols refer as before to the values of P^E defined for passbands B, V, R, I. The errors indicated correspond to the (dominant) systematic errors associated with the uncertainty of the lunar albedo, corrected in the same way as discussed in Sect. 4.2.2. Dashed lines indicate the global fit solution using the empirically determined parameters from Table 2 for the modified Rayleigh function (Eq. (6)) of the phase function (cp. Fig. 6). In Fig. 14 we show only the subsets corresponding to the exact time (resp. phase angle) when the datasets were observed. Typically, the phase angle α changes within three hours by not more than 1°, and the change of P^E expected for this change in α is less then 1%. The straight dashed lines indicate the expected change of P^E only due to changes of α. This change is rather low within 3 h. Within the errors, many of the measured values of P^E in fact are fully compatible with a slow change expected (see, e.g., datasets 2012-12-19 and 2013-02-03). The measurements are not expected to follow the empirical fit exactly. The fits correspond to an average sampling of Earth’s global appearance, and the offsets may thus just hint to deviations of the actually observed scenery from the averaged, global, scenery. The datasets from 2013-02-19, and in particular 2013-02-22, show more variation, in particular in the I band. For both cases, the intermediate measurements were observed with another grism (600I) than those preceding and following (300V). While an offset that is at least partially caused by systematic effects due to the different instrument setup is not excluded, larger amplitude trends and variations are also noted in the other bands, in particular around 01:30UT.

In the same Fig. 14, we also plot the values of ΔPVI and EW(O₂-A) with different symbols, together with their formal errors, in a common scale. Apparently, both parameters show some variability, but are relatively constant in particular for 2012-12-19 and 2013-02-03. Measurements of 2013-02-03 correspond to the “A” sample, and ΔPVI is lower compared to the others. Within the 2013-02-19 datasets, ΔPVI shows excursions to rather low ΔPVI values, indicative for the temporal appearance of a free surface covered by vegetation around UT01:30. This needs to be confirmed if this interpretation holds and a corresponding surface or scenery comes into sight. In the same dataset we also noted variation of EW(O₂-A) of around 10 Å, fluctuating during the duration of the observation sequence. This dataset might be in particular affected by changing sceneries in combination with changing cloud patterns during rotation of Earth.

Fig. 13 O2-A band regionbetween 7550 and 7700 Å for eight high-spectral resolution polarization spectra observed with grism 600I. Spectra are corrected for lunar depolarization, i.e., they should correspond to true polarization values of Earth at these epochs. The fine structure seen in the absorption band is real.

Fig. 14 Short-term variability of polarization features of selected observations in timescales of a few hours. Polarization in the B, V, R, I bands is shown in the usual colors (blue, violet, green, red). The EW(O2-A) and ΔPVI values are plotted with black squares and diamonds in a common scale indicated on the right axis.

5 Conclusions

Earthshine observed with the VLT naturally allows us to divide the sample into two groups. One group contains Earthshine from waxing Earth and contains major contributions from the Atlantic Ocean, the Amazonas region, Europe, Africa, and Antarctica (here called sample “A”). The other group contains Earthshine for the waning Earth and probes in particular the Pacific ocean and has no or little visible land surface (sample “P”). Both observational sceneries are partially covered by variable, and possibly systematically different, cloud patches. They represent different views of planet Earth. Our statistical analysis focused on finding observational properties that may distinguish the two groups. This should constrain the impact of typical surface and atmosphere characteristics on observables extracted from polarization spectra and build up an empirical basis for comparison with theoretical scattering models of Earth-like (exo)planets.

We extracted the fractional polarization in four characteristic wavelength bands from individual Earthshine spectra and constructed Earth phase curves in relation to their phase angle α in Sect. 4.2. Overall, polarization at a certain phase angle is consistent with the values reported in the literature. But there is considerable spread among datasets from different authors, and among our own different observing runsn. We observe variability of polarization spectra on timescales ranging from minutes to data that are separated by weeks to months.

We reduced the largest systematic uncertainty and derive the Earth polarization phase curve after correcting the observed Earthshine data for lunar depolarization efficiency. Notwithstanding, fluctuation around a mean phase curve, which can be well approximated by a modified Rayleigh function over the phase angle region covered by our data, exists at a relative level of 20% around the maximum degree of polarization (i.e., ±6% at P_max = 36%). This spread can, at least for the blue spectral band, be bracketed by models of Earth that have clouds covering 40% of its surface and the other 60% coveredby ocean. Generic theoretical polarization models cannot fully explain the red spectra range, in particular at larger phase angles where P tends to be relatively high and spectrally flat.

We introduced polarization color ratios in the blue, P_B∕P_V, and in the red, P_R∕P_I which largely reduce uncertainties caused by lunar depolarization. Comparing these ratios with simple models demonstrates their sensitivity to cloud and surface properties. Using cloud properties from Earth remote-sensing data at the time of the observations, the Earth models from Emde et al. (2017) appear to be fully compatible with two observations in the blue color ration P_B ∕P_V. In the red, there is less scattering by the gaseous atmosphere, and thus more light reaches the surface. The color ratio in the red, P_R ∕P_I appears difficult to explain with the current models. The explanation of the relatively shallow slopes of Earth’s polarization spectra in the red appears to require more advanced surface reflection in the models.

The difference between the “A” and “P” subsamples with respect to their different polarization phase functions, and relatively higher P_max values (see Fig. 6, and Table 2) for the “P” sample may indicate a slightly lower average cloud coverage when sampling this hemisphere. But also differences of the mean cloud optical thickness and/or altitudes can contribute to systematic differences in the polarization phase functions. Detailed models combined with Earth remote-sensing data tailored to individual observations should help to discover the correct explanation.

Reliable phase curves of planet Earth are sparse, respectively absent. So far, only few polarization values in three bands B, R, I have been derived from the satellite-borne POLDER instrument by Wolstencroft & Bréon (2005) for a phase angle of 90°, extrapolating to whole Earth cloud coverages. Their values for an assumed 55% global cloud coverage are 22.6%(P_B), 8.6%(P_R), and 7.3%(P_I), which is lower than those derived by us, and lower than the models by Stam (2008). But we note the potential advantages of monitoring Earth’s polarization directly from suitable satellites, which allows cross-calibration between Earthshine and direct satellite polarization data, and should help to further validate theoretical models. Because a satellite in a low Earth orbit does not provide an instantaneous view of the whole Earth, instead sampling small regions of the Earth along its orbit, a suitable satellite would preferable be far away, such as in a geo-stationary orbit or even in a lunar orbit (for a description of the advantages of the latter position, see Karalidi et al. 2012b; Hoeijmakers et al. 2016).

The insufficient phase coverage of our data did not allow us to detect signatures of enhanced polarization due to rainbows and/or cloud bows at a scale of global Earth. For water clouds, the primary rainbow is expected for scattering angles around 139°, thus at a phase angle around 41° with a smallspectral dependence. Bailey (2007) and Karalidi et al. (2012a) predicted a factor of two to three enhanced polarization around that phase angle, even with cloud coverage fractions as low as 30%. Unfortunately, the lowest phase probed by our observations was 50°, where no enhanced polarization would be expected (Karalidi et al. 2012a). It will be interesting to probe those small phase angles in the future, albeit the difficult viewing geometry implies very short observing times for Earthshine; the dark region on the moon across which the Earthshine would be observed is indeed very narrow at those phase angles.

Signatures of the glint, i.e., direct sunlight reflected on the (visible) ocean surface can, in principle, occur at most viewing geometries and phase angles, and may be more easily detectable in Earthshine than in rainbows. Williams & Gaidos (2008), Zugger et al. (2010), and Robinson et al. (2010) modeled glint properties on distant ocean-covered planets. Enhanced light curve modulation and/or somewhat higher degrees of polarization signals are predicted. Detection of sunglint on the Earth by the satellite LCROSS was reported in Robinson et al. (2014), adding specular reflection on the ocean surface of the Earth to the inventory of potential diagnostics for a habitable world. Emde et al. (2017) have simulated polarization spectra at phase angles of 87° and 102° corresponding to the observations presented in Sterzik et al. (2012). As planetary surfaces they included a two-dimensional Lambertian surface and an ocean surface, respectively. Their results show that the degree of polarization is up to 10% higher for the ocean surface and they attribute this change to the highly polarized sun glint. In relative terms, the increase of P in the red part of the spectrum is largest, up to 10%. The spread observed in our data may at least be partially affected by variable conditions for glint polarization and we refer to detailed modeling in a forthcoming paper.

We found that the actual vegetation coverage for a given observing scenery, as measured through NDVI, correlates with the ΔPVI parameter derived from the polarization spectra around bandpasses characteristic for the VRE. The ΔPVI thus allows us to probe contributions of surface vegetation on Earth rather sensitively and robustly. The positive detection of a VRE – even at different wavelengths – had been proposed as a potential biosignature for other planets (Tinetti et al. 2006b), but its interpretation as biosignature appears to be problematic as certain types of mineral reflectance could mimic the albedo slope and strength (Seager et al. 2005). Further characterization and potential application of VRE polarization detection is matter of active research (see, e.g., Sterzik et al. 2010; Martin et al. 2010, 2016; Berdyugina et al. 2015; Bagnulo et al. 2015b).

Understanding the O₂-A (and B) absorption band is essential to substantiate its relevance as potential biomarker in exoplanet atmospheres. Rodler & López-Morales (2014) and Lovis et al. (2017) assessed the feasibility to detect this band for the closest exoplanets with current and future high-resolution (flux) spectrographs at the VLT and the Extremely Large Telescope. Our observations stress the potential of spectropolarimetry as a double-differential technique to enhance contrast between the (unpolarized) host star and the highly polarized planet exhibiting a large dynamics in the fractional polarization across this line. Another advantage of using spectropolarimetry to detect O₂ absorption in an exoplanetary atmosphere is that the polarization across the band is independent of the absorption by O₂ in the Earth’s atmosphere, except that the latter absorption lowers the observable fluxes and S/N. Future feasibility studies should be directed to assess the merit employing high-resolution spectropolarimetry to observe this band in exoplanets.

Acknowledgements

Based on observations collected at the European Southern Observatory under ESO programs P87.C-0040 and P90.C-0096. We appreciate the comments of an anonymous referee that helped to improve the paper.

Appendix A: Additional table

References

All Tables

All Figures

	Fig. 1 Left panel: acquisition image of the lunar limb that contains Earthshine observed on 2013-02-04UT08:10:50. Right panel: map of the apparent albedo from Velikodsky et al. (2011) of the same region that has been scaled and rotated. Corresponding FORS2 slitlets are superposed on both images.
In the text

	Fig. 2 Aspect of the Earth for eight representative observing epochs. True color RGB composite images of Earth as seen from the Moon, simulated using the MYSTIC radiative transfer model with cloud data from the ECMWF forecast model and MODISsurface albedo. The epochs correspond to the phase angles during observations.
In the text

Fig. 3 Spectra of the fractional polarization of Earthshine as observed with FORS. For better visibility, the spectra are grouped as follows. Left panels: spectra of the Pacific in the Earthshine (to the west of the observing site), right panels: spectra of the Atlantic (to the east). Upper panels: spectra taken at phase angles α < 90°, lower panels: those at α > 90°. Individual observations are identified by colors specified in the legend, and each spectrum is identified with its observation ID and its actual phase angle in Table A.1.

In the text

Fig. 4 Measured fractional polarization of Earthshine as a function of the phase angle (Sun – Earth – Moon). The different symbols refer to different viewing geometries during the observations: the Pacific or Atlantic ocean. The bandpasses used correspond approximately to the traditional B, V, R, and I bandpasses (see text). The lines connect at least three independent observation cycles that belong to a consecutive sequence of observations within a given observing run (see Table A.1).

In the text

	Fig. 5 Polarization of Earthshine as a function of phase angle α. The different symbols refer to different observations, including observations by Takahashi et al. (2013) and Bazzon et al. (2013). The bandpasses correspond approximately to B, V, R, and I bandpasses, as described in the text.
In the text

	Fig. 6 Polarization of Earth after correction for lunar depolarization effects.
In the text

Fig. 7 Comparison of the fractional polarization of Earth in dependency of phase angle (Sun – Earth – Moon) with models in B, V, R, and I bandpasses. The error bars indicate the effects of uncertain lunar albedo on polarization. The phase curves for three representative models of Stam (2008) are overplotted with different line styles. Values extracted from the simulated spectra of Emde et al. (2017) for observation IDs A.3 and B.4 are indicated by stars.

In the text

Fig. 8 Polarization color ratios PB∕PV (left panel) and PR∕PI (right panel) as functions of α. Earthshine data are indicated by filled symbols. Blue left half-moon symbols indicate the Pacific side and green right half-moon symbols the Atlantic side. Values derived from Moonshine spectra are indicated by nonfilled (open) symbols for the respective sides of the Moon from which the moonlight is scattered. Models of Stam (2008) are overplotted with different linestyles. Values derived from model spectra of Emde et al. (2017) are indicated by stars.

In the text

Fig. 9 Selection of normalized spectra between 6400 and 8000 Å. The horizontally dotted lines indicate the regions across which the vegetation index ΔPVI, associated with the VRE: blueward (between 6750 and 6850 Å; green dotted line) and redward (between 7480 and 7780 Å, but avoiding the region affected by the O2 -A between 7580 and 7680 Å; red dotted line). The difference between the green and red dotted lines defines the amount of ΔPVI. The derivation of EW(O2-A) is explained in the text.

In the text

Fig. 10 ΔPVI as a function of the phase angle with the symbols indicating the viewing sceneries (Pacific or Atlantic). The green and blue lines indicate the linear regression fit from all the observations of either scenery. Confidence bands of 1σ around the linear regression are shaded in gray for the two subsamples. While the Pacific side shows a very weak correlation with the phase angle, the population from the Atlantic exhibits decreasing ΔPVI with increasing phase angles. Both samples are statistically distinct according to a two-sample KS-test. Red lines indicate the ΔPVIs as derived from the models from Stam (2008). The star symbol indicates the ΔPVI as derived from the model for the B.4 from Emde et al. (2017); the value derived from the model for the A.3 dataset is −4.49, and falls outside the plot.

In the text

Fig. 11 ΔPVI as a function of the actual vegetation cover observed. The actual vegetation cover has been calculated from NDVI values for each image pixel. Different symbols indicate the viewing sceneries (Pacific or Atlantic). The black line shows a linear regression fit with all the observations. Confidence bands of 1σ round the linear regression are shaded in gray. The anticorrelation of ΔPVI with the actual vegetation cover is strong, although outliers exist.

In the text

Fig. 12 Equivalent width (EW) in the O2-A band as a function of phase angle α. The symbols indicate the different sceneries (Pacific and Atlantic). Full lines and gray areas indicate the linear regression and its errors. Both regression curves are distinct by ≈ 2σ from each other, and offset over the full range of phase angles covered. A very low probability found in a two-sample KS-test indicates that both samples are not drawn from the same underlying population. Red lines with different line styles refer to models of Stam (2008) and star symbols refer to the two models of Emde et al. (2017).

In the text

	Fig. 13 O2-A band regionbetween 7550 and 7700 Å for eight high-spectral resolution polarization spectra observed with grism 600I. Spectra are corrected for lunar depolarization, i.e., they should correspond to true polarization values of Earth at these epochs. The fine structure seen in the absorption band is real.
In the text

	Fig. 14 Short-term variability of polarization features of selected observations in timescales of a few hours. Polarization in the B, V, R, I bands is shown in the usual colors (blue, violet, green, red). The EW(O2-A) and ΔPVI values are plotted with black squares and diamonds in a common scale indicated on the right axis.
In the text