© ESO 2020

1 Introduction

The presence of liquid water is key for a planet to be habitable (Kasting & Catling 2003), and the remote detection of water on exoplanets will be a rosetta stone in the search for biosignatures and life on other worlds.

On present-day Earth, the hydrosphere contains about 0.023% of Earth’s total mass. Thereof, 97% resides as oceans in the liquid phase, while ≈2% is in the solid phase as ice in glaciers. Only a mass fraction of about 10⁻⁵ is dispersed in the atmosphere as clouds and water vapour. Remote sensing techniques applied to Earth observations take advantage of specific signatures that the liquid water phase imprints on light by reflection on ocean surfaces or by cloud scattering.

Observations of polarisation as a function of the phase angle (or scattering angle) allow a detailed diagnosis of the physical and chemical properties of liquid droplets in a planet’s atmosphere, even for small quantities. Spherical particles produce a strongly polarised peak known as “rainbow scattering” caused by refraction and internal reflection of light in droplets. The rainbow scattering angle β is about 139° for water, and its contrast is much higher in polarisation than in flux. The mechanism is commonly described by Mie scattering and mainly depends on particle size distribution and index of refraction. Rainbows seen from above are often dubbed “cloudbows” and may appear less colourful because of the smaller particle sizes involved.

In the Solar System, the polarised phase curve of planet Venus provides an opportunity to reliably retrieve its atmospheric main constituents (in this case a concentrated solution of sulfuric acid) and sizes (around 1 μm) from ground-based observations, and demonstrates the power of measuring and interpreting phase angles around the rainbow scattering angle (Hansen & Hovenier 1974).

Evidently, this physical mechanism and its detectability by polarisation has applications beyond the Solar System. Bailey (2007) explored the effects of rainbow polarisation in the search for habitable planets. He estimated a disc-integrated rainbow peak between 12 and 15% (fractional polarisation) for Earth if seen as an exoplanet. In order to prepare the quest to find and interpret observations of a “second Earth”, radiative transfer models of Earth-type model planets with clouds have been calculated that include proper treatment of scattering and polarisation. Stam (2008) simulated polarisation spectra over the entire phase curve in order to investigate the complex interaction of light with different surface types underneath an atmosphere which itself hosts clouds. This work was extended by using more realistic and patchy clouds over inhomogeneous surfaces in Karalidi & Stam (2012). Calculations in Karalidi et al. (2012) corroborated the existence of a robust rainbow signature for Earth-like planets, even if clouds are composed from mixtures of liquid droplets and ice crystals.

However, whether or not the rainbow is visible to putative distant observers of planet Earth is unclear. It is also unclear whether or not these observers would be able to derive chemical and physical properties of its atmosphere by remote sensing, as we can do for Venus.

So far, polarisation measurements of the rainbow of Earth from above have only been performed by airborne instruments (Alexandrov et al. 2012) or satellites. Using spaceborne polarisation measurements with POLDER onboard the PARASOL satellite, Bréon & Goloub (1998) provided quantitative polarisation maps of Earth in polarisation. Single scattering by cloud droplets located at the cloud top cause multiple cloudbow features. The inversion of the observed polarised reflectance maps in three spectral bands (0.44, 0.67 and 0.86 μm) resulted in narrow droplet distributions (v_eff ≈ 0.02) and droplet effective radii between 8 and 10 μm. Knowledge of the effective sizes, shapes, and mass concentrations of cloud droplets is highly relevant for the energy budget in the atmosphere of Earth. Droplet sizes relate to cloud albedos, distinguish clean and polluted clouds (Peng et al. 2002), and represent a sensitive input parameter for global climate simulations; see e.g. Gultepe & Isaac (2004).

Furthermore, the spatial and temporal distribution of the thermodynamical phases of clouds and the contributions of mixed-phase and ice clouds for climate sensitivity in the global climate model are highly relevant, but unfortunately not well constrained (Lohmann & Neubauer 2018). In particular, disc-integrated whole-Earth observations are scarce and hard to obtain. Goloub et al. (2000) inferred the presence of an ice phase in the upper cloud layers from POLDER measurements. Wang et al. (2019) used a specific set of PARASOL data to re-construct the diurnal variation of disc-integrated polarisation in its three spectral bands at a phase angle of 55°, but these do not cover the rainbow angle. Their values are generally consistent with expectations from simple models, but are very limited in terms of viewing geometry and therefore do not allow further constraint of cloud or surface properties.

In principle, disc-integrated flux and polarisation information from the whole Earth can be obtained with Earthshine measurements. Earthshine is sunlight scattered by Earth and indirectly reflected from the lunar surface back to Earth, where it can be observed from the ground. This mimics observations of Earth as an exoplanet. Earthshine polarisation observations have been pioneered by Dollfus (1957), who determined Earth’s fractional polarisation for several phase angles sampled between 30° and 140°. While his measurements nicely describe a Rayleigh scattering atmosphere with the expected peak polarisation around quadrature, it is interesting to note that no rainbow feature is present, likely because of insufficient phase sampling.

Modern polarisation measurements of Earthshine allow us to infer new and more detailed insight into Earth as seen as an exoplanet. Using spectropolarimetry in optical wavelengths, Sterzik et al. (2012) derived fractional contributions of clouds and ocean surfaces and was able to distinguish visible areas of vegetation as small as 10% in two different aspect geometries of Earth, both close to quadrature. Takahashi et al. (2013), Bazzon et al. (2013), and Sterzik et al. (2019, hereafter Paper I) extended the phase coverage of Earthshine observations to phase angles between 30° and 140°. All three works are - within the expected variations caused by different sceneries and cloud coverage – consistent with each other and also agree broadly with the models, except at longer wavelengths, where the observed polarisation was always higher than model predictions.

The impact of ice clouds and the potential visibility of the ocean glint on the polarisation spectra of Earth was recently explained in more detail through 3D vector radiative transfer calculations by Emde et al. (2017). Cloud composition, stratification, and patchyness all have significant and observable impact on the strength and spectral shape of Earth’s polarisation spectra, and must be properly included in the simulations. Direct reflection of sunlight either directly on the surface of the ocean, or on cloud decks, also contributes to an enhanced polarisation, and can be expected ubiquitously (Trees & Stam 2019). These simulations of realistic Earth models included the treatment of mixtures of liquid and ice droplets, and predict the existence of a robust cloudbow feature in the polarisation phase curve for adequate choices of size parameters. A framework for simulations of comprehensive models of Earth’s disc-integrated Stokes vectors has also been put forward by García Muñoz (2015).

However, until now, no indication of a cloudbow feature has been observed in integrated whole-disc observations of Earth. Here, we report a clear detection of the cloudbow through polarisation observations of Earthshine for the first time. The measurements allow a relatively robust retrieval of the refractory index and characteristic sizes of the scattering agent, in this case liquid water. The detection of enhanced polarisation caused by cloudbows on other planets may become a realistic possibility and would contribute to characterisation of the composition of exo-atmospheres and help to constrain bio-markers in other worlds.

2 Observations

We used the FORS2 instrument mounted at the Cassegrain focus of the Antu telescope of the VLT during four suitable monthly observing time windows closely following the new moon from late October 2019 to late January 2020. The FORS2 instrument with its wollaston prism and rotating retarder waveplate enables direct measurement of the (wavelength-dependent) quantities P_Q = Q∕I and P_U = U∕I, and thus the fractional polarisation P(λ) (or degree of polarisation, referred to in this work as simply “polarisation”, as a percentage) defined as $$P=\sqrt{\left(P_Q{}^2+P_U{}^2\right)}.$$(1)

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

Observations were obtained with grism 300V, which spans the optical range between 3600 and 9200 Å, resulting in a spectral resolution of ≈220 with a 2′′ slit width. The slit was always oriented to cross the eastern part of the lunar limb. Half of the available detector therefore contains the Earthshine signal on the lunar surface, while the other half contains empty sky. This allows to linearly extrapolate and substract the sky background that otherwise contaminates Earthshine, as explained in Hamdani et al. (2006).

The observation dates were chosen to be able to observe Earthshine at phase angles α of the Sun–Earth–Moon (S–E–M) geometry that allow sampling of the suspected cloudbow feature around α ≈ 40° (α relates to the scattering angle β as α = 180° − β). Observations in this geometrical configuration are challenging, and are only possible during a short interval of time of only one evening twilight per lunar cycle. Observations shortly after evening twilight in Chile sample large portions of the illuminated Pacific ocean. We chose four observing epochs, listed in Table A.1. To maintain consistency with the labels used in Paper I, we denote the new datasets as H.x, I.x, J.x and K.x, respectively. The dates of observations cover phase angles α of 33, 35, 37, 42, 44, 50, 53, 54, 55, 56, 57, 62, 63 and 65 degrees. As we show below, these phases allow us to sample the cloudbow feature roughly around its maximum expected polarisation, albeit on different observing dates. We did not have access to even smaller phase angles for which cloudbow polarisation is expected to cease.

In order to largely suppress systematic uncertainties of the calibration of dual-beam retarder polarimeters, we used the so-called “beam swapping technique”, and observed the target (the lunar limb) with the retarder waveplate orientation consecutively set to position angles 0, 22.5, 45,..., 337.5°. The number of distinct retarder settings is indicated in Table A.1 by the parameter N_cyc and is usually equal to 16. However, occasionally, in particular at the lowest phase angles, observing time was not sufficient to conclude a full cycle of 16 retarder settings. In these cases, only eight or four angles of distinct retarder settings were performed. Albeit with lower signal-to-noise ratio (S/N), this is still sufficient to suppress most calibration errors, and “null polarisation” spectral profiles are always healthy with statistical errors around zero (for their definition see Bagnulo et al. 2009).

The spectra (and therefore the derived values for the fractional polarisation P) for the observations at a phase angle of α = 33° (dataset “K.1”) and α = 35° (dataset “ J.1”) are the most uncertain compared to all other spectra. For these spectra, the total time observing Earthshine during twilight was small, and in particular dataset J.1 may suffer from additional systematic errors due to the acquisition of only two retarder positions to derive the Stokes parameters. In principle, two distinct retarder positions are sufficient to determine the two Stokes parameters Q and U independently, but any systematic error introduced during the observation may then not be cancelled out by the beam-swapping techniques. We estimate these errors by subdividing a different observing sequence (J.3, which has N_cyc = 8) into four subsets of two settings only, and determine their polarisation spectra correspondingly. These subcycles are indicated by J.3-1 to J.3-4 in Table A.1. Typical errors of about 0.5% are introduced in the polarisation, but these errors appear still lower than those introduced by the uncertainties on the lunar depolarisation. These have an effect of 1–2%, as can be seen in Table A.2. We therefore consider the J.1 dataset in the analysis with its higher systematic errors. Although both spectra J.1 and K.1 have comparatively low S/N, we include them in the subsequent analysis (but take into account their higher statistical errors).

We refer to Paper I for a full description of the data-reduction techniques, method of background subtraction, and correction for lunar depolarisation. The last two steps are crucial for a proper analysis. Due to a viewing geometry close to the new moon, the contamination with moonshine background is even lower than in the datasets presented in Paper I. Linear background extrapolation and subtraction works well and gives reliable results. Here we also consider data from Paper I that refer to the “Pacific” viewing geometry, and therefore we reprocessed all datasets consistently with minor improvements. We report all values (and their errors) including updated values for datasets B.x, E.x, and G.x in Tables A.1 and A.2. Differences between the old and new data reductions are typically below 1% for Earthshine polarisation, and always within the systematic errors for Earth polarisation.

3 Results

3.1 Polarisation spectra

Earthshine polarisation spectra are shown in Fig. 1 for our four new observing epochs. Spectra are similar and tend to decrease with increasing wavelength and to increase for increasing phase-angle in the blue spectral range; they exhibit similar polarisation values in the red part of the spectrum, independent of phase-angle. Lower phase angles (37° and 42°) show polarisation that is only weakly dependent on the wavelength in contrast to higher phase angles. The polarisation spectra for α = 33° (J.1) and 35° (K.1) are spectrally rather flat, albeit noisier due to low S/N; because of a very low flux signal of Earthshine as compared to twilight sky background,their spectra are particularly uncertain around both ends of the available wavelength coverage, and in the region of the O₂ -A line around 7700 Å (which is omitted in Fig. 1 for these two cases).

All spectra have been observed with the grism 300V, and have a spectral resolution of approximately 220 (which corresponds to a resolving power of ≈30 Å), given a slit width of 2′′. In order to better analyse the spectra, we derived the degree of polarisation of Earthshine P^ES for the four bandpasses B, V, R and I, averaged over the wavelengths P_B: 4350–4550 Å, P_V : 5450–5650 Å, P_R : 6450–6650 Å, and P_I : 8050–8650 Å. The simultaneous measurement of Stokes Q and U also allows us to calculate the angle of polarisation ϕ. This can be compared to the angle Φ between the normal of the scattering plane and the celestial north pole. As shown in Paper I, we find very good agreement between these values, which are relatively independent of wavelength. For cases J.1 and K.1, the differences are larger and of the order of 10°. This may be due to residuals of twilight sky that might not have been fully removed and contaminates the spectra.

3.2 Polarisation phase curves

The fractional polarisation of Earthshine P^ES extracted from the spectra for the four bandpasses (as listed in Table A.1) is plotted as a function of phase angle α in Fig. 2. Different colours indicate different bands B, V, R and I. Symbol sizes are always larger than the statistical errors calculated from the spectra, except for datasets J.1 and K.1. For phase angles larger than 45°, polarisation P is always higher in blue than red, and is increasing for increasing α in all bands. Even for a narrow range of phase angles of a few degrees there is significant variation of P (up to 2% in the blue and up to 1% in the red) when observed at different epochs. All these findings are qualitatively consistent with observations at higher phase angles reported in Paper I.

P^ES reaches minimum values in all bands for α around 45°, ≈5% in B and ≈2% in R. For further decreasing α, P increases again for all bands. Observations between α ≈ 35° and 37° appear to sample a local peak, which reaches $P_B^{\text{ES}}$ = 6% in B, and around 3–4% in R. Observations at α ≈ 33° are distinctively lower, between 2–4% for all bands; the (uncertain) value in B is lowest.

The values of polarisation corrected for lunar depolarisation P^E in the four bands are listed in Table A.2. Two types of errors are given in the table: statistical and systematic. Statistical errors (in brackets) are calculated from null profiles. The systematic errors of P^E come from uncertainties in determining the lunar albedo at the observing slit position. As in Paper I, the relation of a wavelength-dependant polarisation efficiency as a function of the lunar albedo is assumed (Bazzon et al. 2013). We calculate the combined errors of P^E as the root of squares of the sum of statistical and systematic errors.

P^E and their errors are plotted in Fig. 3 as a function of the phase angle α. While the lunar depolarisation correction mainly re-scales the fractional polarisation P^ES by a factor of ≈5∕2, it also introduces a slight modification of the spectral shape due to its weak wavelength dependence.

The general shape of the phase curves is similar for the polarisation of Earthshine P^ES (Fig. 2) and Earth P^E (Fig. 3) in all four bandpasses. In order to help to guide the eye to follow the general trends of the phase curve in each band, the full lines in Fig. 3 show a low-order polynomial fit to the data in each bandpass. The local minimum around α ≈ 45°, the increase towards the local maximum at α ≈ 35°, and the decline for α ≈ 33° are evident.

Following the spectral shape, the fractional polarisation for α ≥ 35° is always ordered according to$P_B^{\text{E}}>P_V^{\text{E}}>P_R^{\text{E}}>P_I^{\text{E}}$. Interestingly, for α = 33° this order is reversed, and polarisation slightly increases with increasing wavelength: $P_I^{\text{E}}~P_R^{\text{E}}~P_V^{\text{E}}>P_B^{\text{E}}$.

Fig. 1 Spectra of the fractional polarisation of Earthshine obtained with FORS2 in the 2019-10 to 2020-01 observing campaigns. Individual observations are identified by colours specified in the legend, where each spectrum is denoted by its observation ID and its actual phase angle α in Table A.1.

Fig. 2 Fractional polarisation of Earthshine PES as a function of the phase angle (Sun–Earth–Moon). The bandpasses correspond approximately to traditional B, V, R, and I bands (see text). Lines connect observation cycles that belong to sequences of observations during consecutive nights, while different line styles differentiate observations separated by months (see Table A.1).

Fig. 3 Fractional polarisation of Earth PE after correction for lunar depolarisation effects as a function of the phase angle (Sun–Earth–Moon). Colours indicate the values for different bandpasses B, V, R, and I. The lines fit a low-order polynomial to the data in each band, and help to guide the eye to follow the general trends of the phase curves between 35° and 70°.

4 Comparison with models

The local minimum of fractional polarisation in the phase curve of Earth observed with high accuracy around 45° and the local maximum around 35°–37° are characteristic for a primary rainbow. Their quantitative explanation requires multiple scattering radiative-transfer simulations, which include optical properties of cloud droplets calculated following Mie theory, as for example employed for the interpretation of the polarisation phase curve of Venus (Hansen & Hovenier 1974). In the following, we use increasingly complex models to study if and how the observed phase curves allow us to retrieve the main physical and chemical properties of the atmosphere observed.

4.1 Heuristic model

First, we introduce a simple model with two physical components. For the Rayleigh component, we assume a three parameter model following Korokhin & Velikodsky (2005), which was already applied to fit the phase curve in Paper I; see Sterzik et al. (2019): $$P_{\text{Ray}}(\alpha ,\lambda )=\frac{({\sin }^2{(\alpha -\Delta \alpha (\lambda ))}^{W(\lambda )}}{1+{\cos }^2(\alpha -\Delta \alpha (\lambda ))+dePol(\lambda )}.$$(3)

The three parameters Δα, W, and dePol of Eq. (3) parameterise a possible phase shift, skewness, and depolarisation to be accounted for in realistic atmospheres. All three parameters may depend on wavelength, and may therefore be different for each observing bandpass. The parameter dePol, in particular, reflects which maximum polarisation can be attained close to quadrature in the presence of a non-zero surface albedo. Multiple scattering, for example in opaque clouds, also tends to increase dePol and reduce the fractional polarisation at all phase angles.

Rayleigh polarisation in Eq. (3) is small and approaches zero for both low and high phase angles. On the other hand, a Mie component adds to total polarisation specifically at low phase angles, while it is negligible at high phase angles. We may therefore add-up both components independently, and approximate the total polarisation using relation Eq. (4): $$P_{\text{R+M}}(\alpha ,\lambda )=P_{\text{Ray}}(\alpha ,\lambda )+c(\lambda )\cdot P_{\text{Mie}}(\alpha ,\lambda ,\langle r_{\text{eff}}\rangle ,\langle n_{\text{re}}\rangle ),$$(4)

where parameter c(λ) ensures the relative weighting of the Mie and Rayleigh components, and may be wavelength dependant. An analogous parametrisation was also used by Alexandrov et al. (2012) to retrieve cloud properties using polarised reflectance measurements by airborne instruments.

Mie scattering in water droplets in clouds contributes to the polarisation curve in particular at smaller phase angles. For single scattering, the Mie component can be directly calculated by the scattering matrix elements that depend on the scattering angle, the wavelength, and the properties of the scattering material, in particular its size distribution and its refractory index; see eg, Hansen & Travis (1974) or Bailey (2007).

We calculated the scattering phase-matrices using the latest version of libRadtran (2.0.3)¹, a library of programs to calculate radiative transfer problems in the atmosphere of Earth (Mayer & Kylling 2005; Emde et al. 2016). Its mie utility implements an efficient and well-tested Mie code written by Wiscombe (1980) and allows one to calculate the single-scattering phase matrix of spherical particles for a variety of sizes and size distributions. In addition, the optical properties like the (real and imaginary) index of refraction, and their dependence on wavelength can be specified.

Next, we assemble a grid of Mie models to test if the different scattering properties for populations of particles with different sizes and refractory indices have an observable effect. We aim to retrieve those parameters that fit the observed dependence of polarisation on the phase angle best. This process should allow an estimate of the quality of the retrieval mechanism when compared to those parameters expected and characteristic for Earth.

Particle radii r follow a standard Γ distribution with n(r) ∝ r^α exp−Γr. According to this distribution, effective radii ⟨r_eff⟩ and their characteristic size distribution widths, ⟨v_eff⟩, are determined for the entire population. We systematically vary ⟨r_eff⟩ in the range of 1 and 14 μm (in steps of 1 μm), and probe two typical size distribution widths (⟨v_eff⟩ equal to 0.1 and 0.01).

We use the well-known index of refraction ⟨n_re⟩ of water (having a mean ⟨n_re⟩ = 1.33 at a respective temperature of T = 300 K) and consider its imaginary part and wavelength dependence as used in the code of Wiscombe (1980). The properties of the cloud top layers that cause the cloudbow scattering of Earth should therefore be represented by these reference values. In addition, we construct other sets of models with different refractory indices, varying the real part of refraction from 1.24 to 1.44. We leave the imaginary part equal to zero, and assume no variation of the index of refraction with wavelength for these cases. The different indices of refraction are supposed to bracket a broad range of material and solutions of possible atmospheric constituents in an astrophysical context. For example, liquid methane has ⟨n_re ⟩ = 1.24 at a temperature of ≈100 K, which is characteristic of Titan (Badoz et al. 1992). On the other hand, sulfuric acid has ⟨n_re ⟩ = 1.44 for a mixture of 75% H₂ SO₄ and 25% H₂O, which is characteristic of Venus (Hansen & Hovenier 1974).

In Fig. 4 we show the ratio of the scattering matrix elements − P₁₂∕P₁₁ as examples of Mie calculations for different indices of refraction, effective widths, and three different mean particle radii. As expected, the primary rainbow peak is prominent for all cases, and shifts its peak of maximum polarisation towards smaller phase angles for higher indices of refraction. The secondary rainbow appears at larger phases angles, and is increasingly prominent with increasing mean particle sizes. For the refractive index of water, we compare two size distribution widths: ⟨v_eff⟩ = 0.1 is often considered as a “standard” distribution width of water droplets in Earth clouds. ⟨v_eff⟩ = 0.01 is narrower, and the different orders of the single scattering rainbow peaks appear slightly more pronounced. We note that − P₁₂∕P₁₁ is the fractional polarisation in the (unrealistic) case of single scattering only. Nevertheless, we use it equivalent to P_Mie in Eq. (4), heuristically including multiple scattering effects through dePol and c.

The grid of Mie-models consists of 17 different choices of refractory indices, each calculated for 14 different mean particle radii from 1 to 14 μm, and for four wavelengths (450 nm, 550 nm, 650 nm and 850 nm, representing the bandpasses derived from the observed spectra in B, V, R and I). We use ⟨v_eff⟩ = 0.1 for all values, and in addition ⟨v_eff⟩ = 0.01 and 0.02 for the models with an index of refraction of water.

Final models are constructed according to Eq. (4). For each model, a χ² fit is performed in such a way as to minimise the deviation of Eq. (4) to the data simultaneously in all four wavelength bands. We consider all measurements of the fractional polarisation obtained so far for the same geometrical viewing aspect of Earth (Pacific) from Paper I, and listed in Tables A.1 and A.2. All data points together sample the phase curve from 33° to 138° and contain information about the Rayleigh and the Mie part. Parameters Δα, W, dePol, and c in Eq. (3) are allowed to vary for each band. Those combinations of parameters that minimise χ² are then calculated for each model in the grid. In total, we have 144 independent data points (36 phase angles for each of the four wavelength bins). Four parameters are free for each wavelength, which yields a total of 128 degrees of freedom to calculate reduced χ² values for each model.

Figure 5 displays the results of the four best models according to the χ² -fit procedure, together with the data points of observed fractional polarisation. Data points at low phase angles from this work are indicated with larger symbols compared to others from Paper I.

All four best-fit models appear quite similar, except at very low phase angles below 30°, which are not accessible by our Earthshine measurements. While the width of the distribution is not constrained by the data, a refractive index of (or around) water and characteristic droplet sizes around 6–7 μm are retrieved. It is reassuring that models with a very narrow range of indices of refraction around water (1.33 ± 0.01) are favouredand well constrained given the sensitivity of the rainbow scattering angle on the refractive index.

In addition, the overall shape of the curves represents the measured values quite well for all four bandpasses. The peak value of 15% (in B) and 10% (in R) around the primary rainbow angle is reasonably well retrieved, and the shape of the fitted curves around the rainbow angle is within the uncertainties.

In Table 1 we summarise the fit parameters for the best-fitting model, which are characterised by the refractive index of water and an effective droplet radius of r_eff = 7 μm. For this model we estimate the variation of the fit parameters Δα, W, dePol and c that minimise the χ² deviations, considering the error bars of the data, which are mainly due to the inaccurate determination of the lunar depolarisation factor and are listed in Table A.2. The upper (and lower) bound are derived from assuming the upper (and lower) polarisation values of the observed data set.

It is interesting to note that the parameters for Δα, W, dePol, and c typically do not change by more than 10% among the best-fitting models, which may indicate the robustness of the χ² minimisation techniques applied. The parameter c – indicating the Mie contribution to the polarisation – is only slightly dependant on wavelength, and varies between 12% (in I) and 15% (in V).

The goodness of fit quality for all models is displayed in Fig. 6. The contour levels indicate the same regions of reduced χ² in the grid of models that vary the index of refraction (x-axis) and the mean effective particle radius (y-axis). The lowest contour levels of one-σ are indicated by dark blue, and “banana-shaped” region indicates the best-fitting models according to this procedure. The best-fitting models have a narrow region for the indices of refraction around 1.32–1.33 and effective radii from 4 to 8 μm.

We conclude that a simple, heuristic, two-component model for the polarisation phase curve according to Eq. (4) describes the global polarisation of Earth rather well. It allows robust retrieval of the refractive index and droplet size of liquid water clouds. The most likely material of the scattering agent consistent with our measurements is liquid water with droplet sizes around 6–7 μm.

Fig. 4 Examples for calculations of the Mie contribution due to single scattering. The different panels refer to different indices of refraction ⟨nre⟩. Each panel plots the ratio of the single scattering phase element − P12∕P11 as a function of phase angle, for four wavelengths 450 nm, 550 nm, 650 nm, and 850 nm (B, V, R, and I). Different line styles indicate different mean particle radii distributions for ⟨reff⟩ = 4, 8, and 12 μm.

Fig. 5 Fractional polarisation of Earth as a function of phase angle (Sun–Earth–Moon) with models in B, V, R, and I bandpasses according to Eq. (4). The error bars indicate the effects of uncertain lunar albedo on polarisation. The phase curves of four representative models with minimum χ2 are overplotted with different line styles, together with values of their characteristic sizes and refractory indices.

Fig. 6 Contour map of reduced χ2 to fit the fractional polarisation for sets of models using Eq. (4). The regions of minimum χ2 ≤ 1 are colouredin dark blue.

4.2 Comparison with radiative transfer models

More realistic models of a planetary atmosphere with opaque clouds and reflecting surfaces must include multiple scattering and radiative transfer to correctly calculate the wavelength-dependant Stokes vectors. A full three-dimensional geometry, inhomogeneous surfaces, molecular absorption, and complex micro-physical cloud properties all pose special challenges for solving the problem for all wavelengths and in all directions, which is required to calculate the phase curve. Several simplifyingassumptions are usually made to simulate polarisation spectra of planets in general, and of Earth in particular.

Models by Stam (2008) allow the linear polarisation spectra of Earth to be simulated for all phase angles and resolve spectral lines. This latter author approximates Earth’s inhomogeneous surface albedo (e.g. that of ocean and clouds) by weighted sums of flux reflected on horizontally homogeneous surfaces. Surfaces are covered by an anisotropic Rayleigh scattering atmosphere which itself contains Mie scattering elements that mimic simple water clouds (of a fixed optical depth of ten). This approach can be used to approximate an arbitrary mixture of different components by building linear combinations of weighted sums of polarised fluxes from each constituent coming from scattering and reflection from horizontally homogeneous surfaces. This method allows the ratio of clouded to non-clouded surface regions found in Earthshine observations to be constrained in particular, as demonstrated in Sterzik et al. (2012) for phase angles near quadrature. In principle, Stam’s models also contain signatures of rainbow (cloudbow) that can be compared with our data. While the models show important qualitative similarities to our data, they were not intended to, and do not quantitatively fit the phase curve of Earth for any realistic linear combination of clear ocean surfaces and homogeneous water clouds.

The explanation of cloudbows requires a more realistic description of water clouds, which adds considerable complexity to these models. Clouds on Earth are typically a mixture of high-altitude, moderate-optical-depth ice clouds, and low-altitude, large-optical-depth water clouds, each with different and characteristic scattering properties. For example, models of the Earthshine polarisation spectra from Sterzik et al. (2012) require the presence of a thin layer of high-altitude ice clouds which effectively increases the polarisation in the red parts of the spectra, as demonstrated in Emde et al. (2017) for phase angles around quadrature. On the other hand, ice clouds above liquid water clouds are expected to dampen the effect of a rainbow in the polarisation phase curve, but the feature may still be detectable according to the simulations of Karalidi et al. (2012). In addition, clouds on Earth do not appear as homogeneous entities, but are highly structured and geometrically constrained on many scales in horizontal and vertical dimensions. Trees & Stam (2019) account for patchy clouds above ocean surfaces for which the effect of a sunglint (specular reflected light on the ocean surface) boosts the polarisation spectra at all wavelengths. Their models include water clouds of an optical thickness of 5, with droplets having a standard size distribution (⟨r_eff⟩ = 10 μm and ⟨v_eff⟩ = 0.1). However, the resulting phase curves overestimate the maximum polarisation at the rainbow angle. Their models predict a maximum polarisation of more than 20%, which is significantly larger than the percentage that we observe. While their maximum polarisation at the cloudbow peak is relatively independent of the fractional cloud coverage, their polarisation phase curves for angles ≿ 60° are more sensitive on the fractional cloud coverage, as expected for its impact on the Rayleigh part at higher phase angles.

As demonstrated with the heuristic model, our observations of cloudbow are sensitive on the index of refraction and the mean effective droplet radii. In addition, multiple scattering in clouds, and thus the average optical thickness of water cloud decks, is expected to have a significant impact on the quantitative agreements of the observed and modelled phase curves. In the following, we therefore build a grid of models following the approach of Stam (2008) using fractional weighting and linear combinations of different homogeneous surface types. Simulations are performed with the spherically symmetric version of the radiative transfer model MYSTIC adapted to Earthshine (Emde et al. 2017; Mayer 2009; Emde et al. 2010). MYSTIC is a versatile Monte Carlo code to solve radiative transfer problems in the atmosphere of Earth, and is also available as part of the libRadtran package. The output consists of a polarisation spectrum with all four Stokes parameters (I, Q, U, V).

We simulate three generic constituents embedded in or below a mid-latitude summer standard atmosphere. As ground surface we assume a reflecting open ocean with a realistic bi-directional polarisation distribution function (Mishchenko & Travis 1997) including foam caps and shadowing effects by waves (excited by a constant wind speed of 10 m s⁻¹). As a second component we simulate liquid water clouds located in an atmosphere at a height between 3 and 4 km. We calculate models with different values of optical thickness τ_wc of 5, 7, 10, 15 and 20 in order to analyse the effect of different optical depths of water clouds. Water droplet parameters are assumed as above: spherical droplets follow a Γ distribution with an effective radius ⟨r_eff⟩ = 6 μm and a distribution width of ⟨v_eff⟩ = 0.1. For comparison we also probe effective radii of ⟨r_eff⟩ = 5, 7 and 12 μm. As a final component we simulate ice clouds located at a height of between 10 and 11 km with an optical thickness of one. For the optical properties of the ice cloud, the parameterisations for a general habit mixture from Baum et al. (2014) were used, with ⟨r_eff⟩ = 30 μm. Polarisationspectra for each model and for each phase angle between 1 and 179° were then calculated in steps of 1°. We used 10⁶ photons for each Monte Carlo run, typically lasting 30–60 min of computing time on a standard PC for a single phase angle, with larger computing times for larger cloud optical thicknesses. Spectral importance sampling (Emde et al. 2011) and the REPTRAN molecular absorption parametrisation (Gasteiger et al. 2014) yield a medium spectral resolution of about 1000 for wavelengths between 4000 and 10 000 Å. Efficient Monte Carlo simulations including realistic scattering phase matrices of clouds require sophisticated variance reduction methods, which are also included in MYSTIC (Buras & Mayer 2011).

We build a grid of models by systematic linear combinations of different fractions f_x for each of the three generic constituents, namely: water clouds (wc), ice clouds (ic), and open ocean (oo), and denote these fractions with f_wc, f_ic, and f_oo, respectively. We vary the fractions f_x in steps of 0.1 such that the boundary condition f_wc + f_ic + f_oo = 1 is fulfilled. Stokes parameters (I, U, Q) for the final model are combined as follows: $$I(Q,U)=f_{\text{wc}}\cdot I{(Q,U)}_{\text{wc}}+f_{\text{ic}}\cdot I{(Q,U)}_{\text{ic}}+f_{\text{oo}}\cdot I{(Q,U)}_{\text{oo}}.$$(5)

The degree of polarisation for the final model is calculated as usual by $P=\sqrt{(Q^2+U^2)}/I$. The model grid thus encompasses a total of 66 combinations of the three different components, with each scenery for all 179 phase angles. Each model spectrum is then compared to the observed polarisation at a specific phase angle using only data points of the current observations for phase angles between 30° and 70°. For the simulated spectra, we derive the degrees of polarisation in the same four passbands as for the observations to avoid biases. We interpolate the phase angle observed from neighbouring grid points if necessary, and calculate the difference of the measured from the simulated polarisation for each bandpass and phase angle. In this way, χ² -values are obtained for each model with 84 degrees of freedom (21 distinct phase points × 4 bands).

In Table 2 we list the six best models that minimise the χ² values for all models computed in all four passbands simultaneously. All goodmodels require ⟨r_eff⟩ = 6−7 μm. Smaller or larger ⟨r_eff⟩ yields a larger χ². This result corroborates and confirms the finding of the best-fitting droplet size parameter already obtained in Sect. 4.1 with the heuristic model.

However, the simulations put additional constraints on the optical depth parameter of the water clouds τ_wc and the relative fractions of the generic constituents, f_wc, f_ic, and f_oo. Most importantly, all models with τ_wc < 10 yield significantly poorer fits. For example, the best-fitting model with τ_wc= 7 already increases χ² = 2.7, i.e. a quality significantly below the models with larger τ_wc. Even lower τ_wc = 5 yields a χ² of only ≈5.

The results of the best-fitting model are displayed in Fig. 7 for all passbands, together with the data. The best model has a fractional contribution of 30% water clouds, 60% ice clouds, and 10% ocean, and fits all four bandpasses with a reduced χ² = 2.12. We note that the next five best-fitting models (all with χ² ≤ 2.3) are all constrained to fractions of about 30% water clouds, while the relative fractions of ice clouds and the open ocean surface fractions are much less certain. It is interesting to note that the water cloud fraction and its optical depth can be constrained relatively independently. One might expect that, for example, a higher cloud fraction and lower cloud optical thickness would yield a similar degree of polarisation, but it appears that for phase angles around cloudbow, the radiance is only mildly influenced by Rayleigh, ice clouds, and water surface, i.e. I_wc≫ I_Ray + I_ic + I_oo and Q_wc ≫ Q_Ray + Q_ic + Q_oo. In this case, the total polarisation P for the cloudbow is determined by Q_wc∕I_wc only, and does not depend on the individual fractions f_i. Therefore, for the phase-angle range considered, P only depends on the optical thickness of the cloudy part, allowing its determination.

5 Discussion

How do the parameters obtained compare to values that are known or expected for Earth? Only Earth allows us to determine globally averaged values from manifold high-quality resources, such as globally observed and model datasets available for weather and climate modelling and forecasting. On the other hand, our underlying Earthshine polarisation observations are affected by relatively poor phase sampling and considerable measurement errors. It is therefore not obvious that parameters retrieved match real and/or known values a priori.

The quantitative determination of the correct refraction index of water for Earth might seem trivial. However, the obtained value with a relative accuracy <1% matches the reference value remarkably well, given the very simplified (but robust) heuristic approach used in Sect. 4.1. This shows that the phase sampling around the cloudbow peak, and in particular the location and reliable observation of the polarisation minimum, is sufficient and allows us to sensitively constrain refraction caused by Mie scattering.

The quantitative determination of correct mean droplet sizes and the water cloud optical depths is less obvious. While their values are known qualitatively, they – unlike the index of refraction – vary spatially and temporarily over Earth. Moreover, these parameters are not simple observable quantities for which validated, globally averaged values are available.

In order to compare our values, we derive mean droplet sizes and cloud optical depths based on ECMWF meteorological data, available as ERA5² hourly data for 37 pressure-altitude levels on a spatial grid of 0.25° × 0.25°. The ECMWF meteorological data offer no direct information on the effective radius of cloud droplets and ice particles, but only contain data fields for liquid and ice water content. The water cloud effective radius r_eff can therefore only be expressed indirectly. It is usually expressed as a function of liquid water content L and droplet number concentration N_d using the relation (Martin et al. 1994; Wood 2000): $$r_{\text{eff}}={\left(\frac{3E_{\text{d}}L}{4\pi {\rho }_{\text{w}}{N}_{\text{d}}k}\right)}^{1/3},$$(6)

where ρ_w is the density of liquid water and k is a constant that relates the mean and effective droplet sizes. For water surfaces, k has a value of 0.77. N_d depends on the aerosol number concentration N_a below the cloud base, which can be calculated as a function of the aerosol mass concentration q_a (Boucher & Lohmann 2017; Lowenthal et al. 2004), which itself depends on the surface wind speed W (Erickson et al. 1986; Genthon 1992), and whose zonal and meridional components are ERA5 products. The effective radius enhancement factor E_d (Wood 2000) is dependant on the ratio of drizzle water content to cloud water content and its value is slightly greater than one on average. In addition, the fractional cloud cover f_i is available for each vertical layer. Altogether, the optical depth τ_i and the cloud effective radius r_eff,i can be calculated per layer i and can be vertically averaged. The vertically averaged cloud cover f_wc can then bereconstructed taking into account only those layers that have optical depths τ above a certain threshold. For each observing epoch, we calculated ⟨r_eff⟩, τ_wc, and f_wc. The values for τ_wc and f_wc can be calculated taking into account only those layers that have optical depths above a certain threshold. We use two thresholds, τ_i > 0.1 and τ_i > 2, the latter being indicative of cloud cover without the (minor) contributions of translucent and thin clouds. Spatially averaged values were calculated for an area roughly encompassing the Pacific ocean. Results are summarised in Table 3.

Table 3 shows that the values derived from ERA5 data are quite narrowly distributed around their mean values. The optical depths of the liquid clouds are mostly in the range of 14–17, with mean effective droplet radii of ⟨r_eff⟩ = 6.28 ± 0.11 μm. Obviously, cloud fractions are different for different threshold values. Considering thicker clouds only in each vertical layer, $f_{\text{wc}}^{>2}=0.34\pm 0.02$. These values compare favorably to those retrieved in Sects. 4.1 and 4.2.

On the other hand, globally averaged values of water droplet sizes on top of cloud decks have not been directly observed so far. Experiments carried out by airborne or satellite instruments are sensitive to more regional effects, and typically retrieve somewhat larger sizes of about 10 μm (Mayer et al. 2004; Bréon & Doutriaux-Boucher 2005; Polonik et al. 2020, see e.g.). Particularly over sea surfaces, effective radii up to 11 μm have been derived from ISCCP data (Han et al. 1994), but smaller sizes of 5 μm have also recently been been reported (Alexandrov et al. 2018). Research is ongoing as to how the droplet sizes change at cloud tops due to mixing and evaporation. The impact of changes in droplet size on the global albedo and energy balance of Earth remain to be investigated.

We note that the water cloud fraction and the optical depths of water clouds also appear consistent with the findings of Stubenrauch et al. (2013), who provide an assessment of global cloud properties derived from various different Earth observation satellites.

Fig. 7 Fractional polarisation obtained from the best-fitting radiative transfer model simulated with MYSTIC. Surfaces are assembled as homogeneous mixtures of water clouds, ice clouds, and clear ocean, as explained in the text. The best model yields ⟨reff ⟩ = 6 μm, τwc = 15, fwc = 0.3, fic = 0.6 and foo = 0.1.

6 Conclusion

We observed polarisation spectra of Earthshine at phase angles sampling the characteristics of rainbow and cloudbow scattering features caused by water droplets at the top of liquid clouds. The cloudbow feature seen above the Pacific ocean is very prominent, reaching a maximum measured polarisation of about 14% at a phase angle of 37° in the B band.

These observations can be explained by Mie scattering from water clouds in Earth’s atmosphere. A simple heuristic model with single-scattering Mie cross-sections in combination with a Rayleigh scattering atmosphere explains the degree of polarisation measured in all optical passbands for phases between 30° and 140°. Without a priori knowledge of the physical and chemical scattering components, it sensitively constrains the refractive index and the mean droplet sizes. The derived values are consistent with the properties of clouds on Earth. In particular, the cloud droplet effective radius determined to be around 6–7 μm is of specific relevance; it is the first direct measurement of water droplet sizes ⟨r_eff⟩ at cloud tops above large portions averaged over the Pacific ocean.

In principle, this retrieval method may therefore also be applicable to other planets beyond the Solar System, provided their phase curves can be measured with similar accuracy. However, in these cases, the effects of higher measurement errors and/or a poorer phase coverage must be considered. These latter could for example be modelled using Markov chain Monte Carlo methods.

We increased the fidelity of the heuristic method by detailed vector radiative transfer calculations including multiple scattering in a spherically symmetric approximation of a reflecting ocean surface beneath a standard atmosphere with water and ice clouds. The best models allow us to consistently fit the four bandpasses of the polarimetric phase curve over the visible wavelength regime with good quality. We can corroborate the effective radius of water droplets obtained from the simple model. Moreover, these simulations allow us to infer cloud fractions and their optical depths, given the good sampling of angles across the Mie- and Rayleigh-scattering-dominated regimes in the phase curve. The values retrieved, in particular for water clouds, are consistent with those from Stubenrauch et al. (2013) and those derived from ECMWF-ERA5 models. We find a higher uncertainty in the fractions of ice clouds and ocean. This might decrease when taking into account even larger phase angles which are more sensitive to the sun glint.

The models discussed are simply an approximation of the real situation of Earth during our observations. We are aware that the patchiness of the cloud spatial distribution and the distribution among different cloud phases have to be taken into account to better approximate the observing scenery at any given observation epoch. Therefore, we will pursue three-dimensional Monte Carlo calculations to account for the inhomogeneous surface geometry and the variations of actual cloud cover in the near future.

Acknowledgements

Based on observations collected at the European Southern Observatory under ESO programme P104.C-0048. Observations were carried out in “delegated visitor mode” and we thank the Paranal Science Operations team for impeccable support. Archival data from ESO programs P87.C-0040 and P90.C-0096 were also used.

Appendix A Additional tables

References

All Tables

All Figures

	Fig. 1 Spectra of the fractional polarisation of Earthshine obtained with FORS2 in the 2019-10 to 2020-01 observing campaigns. Individual observations are identified by colours specified in the legend, where each spectrum is denoted by its observation ID and its actual phase angle α in Table A.1.
In the text

	Fig. 2 Fractional polarisation of Earthshine PES as a function of the phase angle (Sun–Earth–Moon). The bandpasses correspond approximately to traditional B, V, R, and I bands (see text). Lines connect observation cycles that belong to sequences of observations during consecutive nights, while different line styles differentiate observations separated by months (see Table A.1).
In the text

	Fig. 3 Fractional polarisation of Earth PE after correction for lunar depolarisation effects as a function of the phase angle (Sun–Earth–Moon). Colours indicate the values for different bandpasses B, V, R, and I. The lines fit a low-order polynomial to the data in each band, and help to guide the eye to follow the general trends of the phase curves between 35° and 70°.
In the text

Fig. 4 Examples for calculations of the Mie contribution due to single scattering. The different panels refer to different indices of refraction ⟨nre⟩. Each panel plots the ratio of the single scattering phase element − P12∕P11 as a function of phase angle, for four wavelengths 450 nm, 550 nm, 650 nm, and 850 nm (B, V, R, and I). Different line styles indicate different mean particle radii distributions for ⟨reff⟩ = 4, 8, and 12 μm.

In the text

Fig. 5 Fractional polarisation of Earth as a function of phase angle (Sun–Earth–Moon) with models in B, V, R, and I bandpasses according to Eq. (4). The error bars indicate the effects of uncertain lunar albedo on polarisation. The phase curves of four representative models with minimum χ2 are overplotted with different line styles, together with values of their characteristic sizes and refractory indices.

In the text

	Fig. 6 Contour map of reduced χ2 to fit the fractional polarisation for sets of models using Eq. (4). The regions of minimum χ2 ≤ 1 are colouredin dark blue.
In the text

	Fig. 7 Fractional polarisation obtained from the best-fitting radiative transfer model simulated with MYSTIC. Surfaces are assembled as homogeneous mixtures of water clouds, ice clouds, and clear ocean, as explained in the text. The best model yields ⟨reff ⟩ = 6 μm, τwc = 15, fwc = 0.3, fic = 0.6 and foo = 0.1.
In the text