Mapping Earth Analogs

Mapping Earth Analogs from Photometric Variability:
Spin-Orbit Tomography for Planets in Inclined Orbits

Yuka Fujii11affiliation: Department of Physics, The University of Tokyo, Tokyo 113-0033, Japan and Hajime Kawahara22affiliation: Department of Physics, Tokyo Metropolitan University, Hachioji, Tokyo 192-0397, Japan

Aiming at obtaining detailed information of surface environment of Earth analogs, Kawahara & Fujii (2011) proposed an inversion technique of annual scattered light curves named the spin-orbit tomography (SOT), which enables us to sketch a two-dimensional albedo map from annual variation of the disk-integrated scattered light, and demonstrated the method with a planet in a face-on orbit. We extend it to be applicable to general geometric configurations, including low-obliquity planets like the Earth in inclined orbits. We simulate light curves of the Earth in an inclined orbit in three photometric bands (0.40.5m, 0.60.7m, and 0.80.9m) and show that the distribution of clouds, snow, and continents is retrieved with the aid of the SOT. We also demonstrate the SOT by applying to an upright Earth, a tidally locked Earth, and Earth analogs with ancient continental configurations. The inversion is model-independent in the sense that we do not assume specific albedo models when mapping the surface, and hence applicable in principle to any kind of inhomogeneity. This method can potentially serve as a unique tool to investigate the exohabitats/exoclimes of Earth analogs.

Subject headings:
astrobiology – Earth – scattering – techniques: photometric

1. Introduction

One of the ultimate goals in exoplanet study is the discovery of extraterrestrial life. Recently, the abundance of Earth-sized exoplanets has been implied from the radial velocity observations (e.g., Howard et al., 2010), the transit observations (e.g., Borucki et al., 2011), and gravitational planetary micro-lensing (e.g., Cassan et al., 2012) and some in so-called habitable zones (HZs) have been already detected. In near future, direct detection of planetary light will play a crucial role in exploring the surface environment and meteorology of planets in habitable zones (HZ planets), and in detecting possible biosignatures there including O, O, HO absorption lines (e.g., Owen, 1980; Angel et al., 1986; Leger et al., 1993; Selsis et al., 2002; Kaltenegger et al., 2010, and references therein) and the vegetation’s red edge(e.g., Arnold et al., 2002; Woolf et al., 2002; Seager et al., 2005; Montañés-Rodríguez et al., 2006; Hamdani et al., 2006; Kiang et al., 2007a, b).

As a test bed for future observations of HZ planets, Ford et al. (2001) simulated the diurnal variation of the scattered light of the Earth. They found 20% (in realistically cloudy cases) or 150% (in cloud-free cases) variations originated from the inhomogeneity of surface components and cloud cover. Their results encourage the possibility to reconstruct the landscape of exoplanets from scattered light curves. Later on, Oakley & Cash (2009) suggested the inversion using the ocean glint at crescent phase and the albedo difference between ocean and land. Their methods work reasonably well when the cloud cover is lower than the actual Earth. Cowan et al. (2009, 2011) performed a principal component analysis (PCA) on the multi-band photometry of the Earth obtained with EPOXI mission. They found that a dominant eigencolor of the Earth corresponds to the cloud-free land component, which allows one to draw a one-dimensional map of continental distribution along the longitude. Fujii et al. (2010, 2011) decomposed of the total scattered light of the Earth into several albedo models —ocean, snow, soil, vegetation, and clouds— and found that the extracted time variation of them roughly recovers the actual surface inhomogeneities including the distribution of vegetation as well as continents/ocean.

Recently, Kawahara & Fujii (2011, hereafter Paper I) proposed a new inversion method named “spinorbit tomography” or SOT, to sketch a two-dimensional map of the exoplanet’s surface using both spin rotation and orbital revolution of the planet. This method directly maps the reflectivity over the surface with the Tikhonov regularization, which balances between the observational noise and spatial resolution of the surface. This was an improvement over Kawahara & Fujii (2010), where the inversion was regularized by the bounded variables least square method and thereby specific albedo models were required. We stress that the SOT retrieves the spatially-resolved image of exoplanets from the disk-integrated light curve without directly resolving of the planet. In this sense, the SOT can be a precursor to the projects which attempt to resolve planetary image directly with very long baseline, such as Hypertelescope (150km array size with 150 apertures of 3 m; e.g. Labeyrie, 1999). In a hypothetical case where a mock Earth with planetary obliquity 90 is in a face-on orbit, they found that the albedo map reconstructed from a single photometric band traces the distribution of clouds. They also showed that cloud-free surface features could be recovered from the difference of two photometric bands since the cloud reflection spectra are roughly constant over optical wavelengths (e.g. Figure 8 of Fujii et al., 2011). The planetary obliquity was also reasonably estimated within the framework of the SOT.

This paper aims to further develop the SOT by considering various configurations, including low-obliquity planets like the Earth in inclined orbits. In the case where the planetary orbit deviates from face-on, the phase angle of the planet is changed and the anisotropy of scattering, primarily by clouds and the atmosphere, can have non-negligible effects. In this paper, we take account of such effects to make the methodology applicable to more general cases.

The organization of this paper is as follows. In Section 2, we define the geometric configuration and describe the methodology of the SOT in detail. In Section 3, we prepare mock light curves to be applied to the demonstration of the SOT. The simulation scheme, adopted data sets, and examples of simulated light curves are found there. Section 4 shows the results of two-dimensional mapping and obliquity measurements via the SOT, and estimates instrumental requirements for future missions. Case studies with different configurations are demonstrated in Section 5 with examples of a zero-obliquity planet in more inclined orbit and a tidally-locked planet (Section 5.1) as well as Earth models with ancient continental distribution (Section 5.2). Finally, Section 6 summarizes the results of this paper and discusses the implications of the SOT.

2. Spin-Orbit Tomography

In this section we describe the methodology of the SOT, which was originally proposed in Paper I. While we described essentials of the SOT in Paper I assuming the face-on orbit, in this paper we extend our method to be applicable for an arbitrary direction of the line of sight. This generalization adds the orbital inclination to the set of geometrical parameters.

2.1. Geometry

Figure 1.— Schematic configurations of the planetary system and surface. Panel (a) shows the spherical coordinate fixed on the planetary surface: the complementary latitude, , and the longitude, . The spin rotation axis is indicated by . We denote the spin motion by , the angle between (, ) and the direction of equinox. We define three unit vectors on the planetary surface (the normal vector, , the vector to the observer, , and that to the host star, ) which specify the arguments of the BRDF : the incident zenith angle, , the zenith angle of scattering, , and the relative azimuthal angle between the incident and scattered light, . Panel (b) illustrates the geometrical configuration in the orbital frame. The orbit lies in the plane, corresponding to the directions of the equinox (-axis) and the solstice (-axis) of the planet. The -axis is defined by the normal of the orbital plane. The obliquity and the inclination are defined as the angle from the -axis to the spin axis and to the line of sight, respectively. The orbital phase measured from the superior conjunction (SC) is denoted by . The orbital phase of equinox is indicated by . Panel (c) displays the geometrical configuration in the observer frame where the longitude of the ascending node is indicated.

Figure 1 illustrates the geometrical configurations. We define the spherical coordinate system, (colatitude) and (longitude), fixed on the planetary surface as shown in Figure 1 (a). The spin motion is specified by , which is defined as the angle between the direction of the equinox and the origin of . Denoting the angular velocity of the spin rotation from observation by , we can write . Note that only changes the prime meridian and have no physical consequence. We will discuss the measurement of the spin rotation by the autocorrelation in Section 4.1. The orbital phase is described by , taking the superior conjunction as its origin. We adopt two parameters to specify the planetary spin axis: the obliquity, , and the orbital phase of the vernal equinox, . The inclination of the orbit is indicated by .

Throughout this paper, we assume that the planetary orbit is already determined, i.e., inclination , the time revolution of , and the phase of superior conjunction are known. Geometrical parameterizations are summarized in Table 1.

Symbol Explanation
Colatitude fixed on the planetary surface
Longitude fixed on the planetary surface
Orbital phase measured from superior conjunction
Phase of spin rotation
(Configuration of the System)
Planetary obliquity
Orbital inclination
Orbital phase of vernal equinox
Table 1Coordinate System and Geometrical Configuration

We introduce three unit vectors to describe the direction of the incident ray and the observer: the vector from the planetary surface to the host star, , the vector from the surface to the observer, , and the vector from a planetary center to the surface, as shown in Figure 1 (a). For mathematical simplicity, we can take an auxiliary coordinates where the -axis is the direction of equinox and the -axis is the orbital axis (Figure 1 (b)), and the components of the three unit vectors are:


where is the clockwise rotation matrix around -axis and .

2.2. Modeling of Scattered Light

Scattering property of each patch of the planetary surface is generally expressed by the bidirectional reflectance distribution function (BRDF), , which is a function of the incident zenith angle , the scattering zenith angle , and the relative azimuthal angle between the incident and scattered light (if the surface have no special direction) as shown in Figure 1 (a). The scattered intensity at each planetary pixel is expressed as


where is the BRDF at wavelength , is the incident flux at , is a planetary radius, and . Note that , , and are functions of all of the geometrical parameters listed in Table 1.

For the reflection model, we assume that the scattering is isotropic (Lambertian). Then the BRDF becomes independent on , , or and


where is the albedo at the planetary surface at wavelength . Under the above assumption, Equation (4) reduces to


where the weight functions for the illuminated and visible area are defined by and , respectively (see also Cowan et al., 2009).

Using the three unit vectors (Equations (1)(3)), we obtain


The explicit forms of Equations (7) and (8) are given in Appendix A.

We define the weight function for the illuminated and the visible area as follows:


With this weight function, Equation (6) reduces to


The integral of the weight function over the planetary surface provides the phase function of a uniform Lambert sphere (e.g., Russell, 1906),


where is the phase angle, defined as


Figure 2 displays the behavior of the weight function in the case of , , and . The weight function covers the region of ( or ) and as the planet rotates around the host star and its spin axis, where is the boundary of described in Appendix A.

Figure 2.— Examples of geometrical dependence of the weight function . Annual variations are indicated by different thin lines. Solid, large-dashed, middle-dashed and small-dashed lines indicate and , respectively. Inner and outer contours indicate and . Thick lines indicate the boundary of the visible area, in Equation (A4). Panel (a) indicates an extreme case: , corresponding to the geometry considered in Paper I. Panel (b) corresponds to the fiducial case of the Earth-twin we consider in this paper: . Since the weight function solely depends on , in other words, longitude minus the phase of the spin rotation, we adopt to the axes in panels (a) and (b) instead of showing diurnal variations explicitly. Panel (c) indicates the upright case or the tidally locked case: (). In the tidally locked case, the axis becomes and does not depend on .

2.3. Solving Inverse Problem

In order to linearize Equation (10), we discretize the surface to small pixels. The center and the area of the th pixel are denoted by and , respectively. Equation (10) is then expressed as

We assume that the pixel size is small enough to satisfy


where denotes the solid angle of th pixel.

We consider bins of observational data. We define the -th data point () as


where is noise of the -th data point.

After these discretizations, we obtain a linear discrete equation


As a pixelization scheme for the planetary surface, Hierarchical Equal Area isoLatitude Pixelization (HealPix) is adopted in this paper111 Unless otherwise noted, the resolution parameter for HealPix is set to 8 corresponding pixel numbers of which is = 768 (Górski et al., 2005).

Note that when the orbit is not circular one needs to caution about the incident flux because it changes in time. However, the procedure of analysis after translating the observed planetary intensity into and the observing time into remains essentially same.

2.3.1 Tikhonov Regularization

For the time being, we assume that all geometrical parameters (i.e. not only but also and ) are known. To solve the linear inverse problem (Equation (16)), we add a regularization term to the misfit function; the regularized solution is given by minimizing


where is variance of the observational error and is the penalty function, which measures the regularity of the solution. The regularization parameter controls effective degree of freedom of the model.

In this paper, we adopt the Tikhonov regularization, also known as the dumped least-square method (e.g., Menke, 1989; Tarantola, 2005; Hansen, 2010). The penalty function of the Tikhonov regularization is:


where is the prior of the model, each component of which is fixed to the average of apparent albedo (e.g., Qiu et al., 2003) times . The solution which minimizes Equation (19) with Equation (20) is given by


where we define and . The orthogonal matrices and are given by the singular value decomposition of and is the th eigenvalue of , that is, the -th component of the diagonal matrix (e.g., Menke, 1989; Hansen, 2010). The is the Kronecker delta.

If data have enough information to reconstruct the model, has little influence on , while if data have almost no information, converges with . From the Bayesian viewpoint, the Tikhonov regularization is equivalent to the parameter estimation with a Gaussian-type prior with the average and the covariance matrix . We briefly summarize the statistics of the Tikhonov regularization in Appendix C.

The probability that each pixel illuminates depends on the obliquity, the direction of the observer, and the latitude of the pixel. During a year of the planetary system, some pixels are frequently sampled and others are not. In order to quantify the relative influence of data on estimation of , we introduce the integrated sensitivity vector (Zhdanov, 2002) defined as


The estimated value of is sensitive to data if is large and vice versa. Clearly, should be equal to if is 0.

The regularization parameter is chosen based on “the L-curve criterion”, which is the maximum curvature point of the model norm versus residuals plot (Hansen, 2010). The details of the L-curve criterion is described in Appendix D.

2.4. Obliquity and Equinox

So far we have described the methodology to reconstruct the albedo maps with known geometric parameters. In reality, the weight function and thus the design matrix are functions of and as well:


While direct imaging will be able to estimate and the superior conjunction, two parameters which specifies spin axis (obliquity, , and orbital phase of equinox, ) are difficult to estimate by any other known techniques. In face-on cases, Paper I showed that these two parameters can be simultaneously estimated via the SOT itself by taking them as free parameters and minimizing 222In the case of face-on orbit discussed in Paper I, was measured from the vernal equinox instead of the point of superior conjunction because the latter cannot be defined, and accordingly the measurement of , the phase of the start of observation relative to the equinox, was to be estimated rather than . . In the same manner, we search for the best-fit and by the Nelder-Mead method with varying . At each trial value of and , the best-fit model vectors are derived by Equation (21). The set of the best-fit , and for different allows us to draw the L-curves, and hence we can finally determine the appropriate by the L-curve criterion.

3. Simulation of scattered light curves of an Earth-twin

In this section, we describe our Earth model to be input to the SOT as well as the resultant light curves.

So far, various authors have carried out simulations of scattered light curves of an Earth-twin (Ford et al., 2001; Tinetti et al., 2006a, b; Montañés-Rodríguez et al., 2006; Pallé et al., 2008; Oakley & Cash, 2009; Arnold et al., 2009; Doughty & Wolf, 2010; Fujii et al., 2010, 2011). In this paper, we follow the simulation scheme described by Fujii et al. (2011); we (1) pixelize the surface of the planet into , (2) assign the parameters of the scattering properties to each patch according to the data obtained with the MODIS (MODerate resolution Imaging Spectroradiometer) onboard the Earth Observing Satellites Terra/Aqua, (3) calculate the radiative transfer in each patch by rstar6b (Nakajima & Tanaka, 1988; Nakajima, 1983), and (4) integrate the emergent light from each patch over the illuminated and visible portion.

Since we do not have a specific instrumental design in mind, we simply assume three m width bands: 0.40.5m (Blue), 0.60.7m (Orange), and 0.80.9m (NIR). Since it is fairly time consuming to sum up the radiance over planetary surface and to integrate fine spectrum over photometric bands, we create a look-up table with relevant parameters in the same way as Fujii et al. (2011) and estimate the radiance at a given pixel by linearly interpolating it.

In Section 3.1, we describe the adopted data sets and our assumptions in more detail. Examples of simulated light curves are shown in Section 3.2.

3.1. Assumptions and Input Data

The rstar6b calculates the radiative transfer through the atmospheric layers given the properties of atmosphere, ground albedo, the direction of the host star, and that of the observer. The vertical profiles of pressure, temperature, and molecule composition in the atmosphere are assumed to be uniform for simplicity; we adopt the US standard model.333U.S. Standard Atmosphere, 1976, U.S. Government Printing Office, Washington, DC, 1976 The thermal profile affects little on the resultant spectra at the wavelengths where the thermal radiation is not dominant.

Since cloud cover significantly changes the reflection spectra, we insert a cloud layer according to the actual cloud data of each pixel. In our model, cloud cover is specified by two parameters: cloud cover fraction and cloud optical thickness . The daily global map for these two parameters is taken from the Terra/MODIS Atmosphere Level 3 Product444 Specifically, we adopt the data of 2008. Other parameters for cloud optical properties are fixed to the value adopted by Fujii et al. (2011). In particular, we assume that clouds are composed of pure liquid water and that the cloud layer is at 3.56.5 km. The altitude of the cloud layer affects the continuum level little (e.g. Figure 8 of Fujii et al., 2011).

The boundary condition of the radiative transfer is determined by the BRDF of the surface. The MODIS project adopts the Ross-Li model to specify the BRDF of land surface (e.g., Lucht et al., 2000) and offers the three coefficients in the Ross-Li model (coefficients of isotropic, geometric, and volume terms). We, however, use the isotropic term only in our simulation to match the specification of rstar6b. We obtain the value of the coefficient of the isotropic term at each pixel from “snow-free gap-filled MODIS BRDF Model Parameters”. These data sets are a spatially and temporally averaged product derived from the 0.05 resolution BRDF/albedo data (v004 MCD43C1)555Available at The original data sets are merged into the adopted resolution by averaging.

In order to incorporate the effect of seasonal variation of snow cover, we adopt the monthly snow cover map offered by MODIS project.666 We replace surface albedo by the reflection spectrum of fine snow for pixels with snow cover fraction larger than 50%.

In short, our model in this paper incorporates the effect of time variation of the cloud cover fraction, cloud optical thickness, surface BRDF, and snow cover.

3.2. Light Curves

Figure 3.— Yearly light curves of an Earth-twin in the case of , , . Signal-to-noise ratio (S/N) is imposed so that S/N after folding the light curves for 6 days becomes 20). Top: light curve of Blue (0.40.5m) band and phase curves of a Lambert sphere with albedo 0.33. Middle: light curve of Orange (0.60.7m) band and phase curves of a Lambert sphere with albedo 0.25. Bottom: light curve of NIR (0.80.9m) band and phase curves of a Lambert sphere with albedo 0.28. The inserted panels enlarge the diurnal variations at and show theoretical light curves by solid lines for reference. In the analysis in Section 4, data within the region bracketed by dot-dashed arrows are used.

Figure 4.— Top: yearly variations of the difference between the reflectivity in NIR (0.80.9m) and that in Blue (0.40.5m), both of which are shown in Figure 3. The geometric parameter set is , , . Bottom: same as the top panel but the difference between the reflectivity in NIR (0.8-0.9m) and that in Orange (0.6-0.7m).

Figure 3 displays the resultant light curves of the Earth in three photometric bands with signal-to-noise ratio (S/N) = in the case of . The inserted panels show the diurnal variations and theoretical (noiseless) light curves in black dotted lines. The variation is due to the spin rotation and reflects the longitudinal inhomogeneity. The amplitude of variation (%) is consistent with former works (Ford et al., 2001; Oakley & Cash, 2009).

The yearly variation is due to the change in the illuminated and visible area (Section 2), i.e., waxing and waning of the planet. In Figure 3, yearly variations of a hypothetical Lambert sphere are also plotted by dashed lines for reference. Simulated light curves substantially exceed the Lambert phase curves at phase , corresponding to the maximum phase angle. The excess at this phase is likely due to the enhanced forward scattering by clouds and atmosphere (e.g., Robinson et al., 2010). This deviation from isotropic scattering possibly causes a systematic bias on the reconstructed map. Unless the change of phase angle is negligible, one must consider this anisotropic effect. In this paper, we avoid this problem in a simple way: we do not use the data at crescent phase since the fractional deviation from Lambert sphere is quite significant at this phase (see Section 4.2 below for further discussion).

Figure 4 shows the differential light curves between two photometric bands (NIRBlue, NIROrange) for later use.

4. Global mapping of an Earth-twin

In this section, we apply the SOT to the simulated light curves. In Section 4.1, we begin with pre-processing of data to be relevant to the SOT by measuring of the spin rotation period. A two-dimensional mapping is discussed in Section 4.2. The obliquity measurement is studied in Section 4.3. The desired aperture size for reasonable mapping via the SOT is estimated in Section 4.4.

4.1. Pre-processing

Figure 5.— Autocorrelation coefficients computed from the first 14 day (red solid), 30 day (green long-dashed), and 60 day (blue short-dashed) mock observations in NIR in the case of (Figure 3).

First of all, we consider translation of the observation time to phase angles . Orbital phase is known once the planetary orbit is determined, which is one of the assumptions in this paper. On the other hand, spin phase needs to be known by measuring spin period . Therefore, we perform auto-correlation analysis on our mock data as suggested by Pallé et al. (2008).

Figure 5 displays the auto-correlation coefficients from 0.8 to 0.9 m data obtained in the first 7 day, 30 day, 60 day mock observations. We find that the spin rotation period of the Earth, 24 hr, is safely measured from any of single-band observations for 1 month, which enable us to assign . Note that the autocorrelation coefficient at (hr) is not approaching 0 as the number of days increases because of the phase variation. When , the phase angle does not change and the autocorrelation is conversing to 0.

Furthermore, we stack the light curves for every 6 days to increase the signal-to-noise ratio per data point. After all we are left with data sets of for each band (=1830) with S/N=20.

4.2. Two-dimensional Mapping

Figure 6.— (a-1) Two-dimensional mapping from one-year observation in Blue band with S/N=20 in the case of without consideration of bias caused by anisotropic scattering, i.e., data at all phases are used in the analysis. (a-2) Same as the top panel but data at are not used in the analysis since the effect of forward scattering of clouds is not negligible. (b) Annual mean of cloud optical thickness. (c) Cloud cover fraction in 2008 March. (d) Two-dimensional mapping from the difference between two photometric bands with S/N=20. (e) Two-dimensional mapping from the difference between two photometric bands with S/N=100. Plots of integrated sensitivity as a function of latitude are attached for respective cases. Note that does not depend on longitude because spin rate is much faster than orbital motion.

We now adopt the SOT to mock light curves of the Earth. First, we naively used one-year light curves without considering the anisotropic effect. The resultant map displayed in panel (a-1) in Figure 6 exhibits a distinct intense region at the North pole and hardly traces any real features of the Earth. This is due to the fact that the reflectivity at large phase angle ( see also Equation (12])) is much larger than expected from Lambertian phase curves (see Figure 3). In reality, the increase in reflectivity at originates from the forward scattering of clouds and atmosphere as mentioned above. This result implies that we cannot ignore the bias from the anisotropy of scattering in the case of inclined orbit. In contrast, the anisotropy mostly leads to the offset for the face-on case () discussed in Paper I since the phase angle is constant.

In order to reduce this bias, we exclude the data with large phase angle . In this paper we use the data that satisfies the condition where the deviation from the phase curve of a Lambert sphere is significant (see also Figure 2 of Robinson et al., 2010). This condition corresponds to (the bracketed region in Figure 3) based on Equation (12).

Panel (a-2) in Figure 6 shows the albedo map reconstructed from the Blue light curve at . The reconstructed map exhibits the highly reflecting zone at high latitude in the north hemisphere. The high-albedo region at can be interpreted as a pattern of clouds and snow, which have high reflectivity in the Blue band. We can compare the reconstructed map with an annual mean distribution of cloud optical thickness (panel (b)) and snow cover fraction (panel (c)). The intensive region in the south hemisphere is not seen in the resultant map. This is a natural consequence because the region near the south pole is insensitive to the data for this specific geometry (). The integrated sensitivity, , to data is plotted as a function of latitude at the right of each map. The interpretation of the reconstructed albedo map should be accompanied with the sensitivity plot.

Surface inhomogeneity is expected to be uncovered by mapping the difference between two photometric bands since the reflection spectra of clouds are roughly constant against wavelengths from visible to near infrared bands (see Figure 3 in Paper I). In this case, we substitute and for and in Equations (16)(21) where subscripts and indicate different bands. We again employ the criterion for usable data because the contribution of clouds is not completely compensated by taking difference of two bands. The reconstructed maps from the difference between two photometric bands with S/N=20 (Figure 4) are displayed in panel (d) in Figure 6. The map for NIRBlue (left panel) can be interpreted as a proxy of the land distribution, since most of the land components have higher albedo than ocean at the near-infrared range (Figure 3 of Paper I). Indeed one can see the bright regions corresponding to the Africa continent and other continents (South America, North America, and Eurasia). On the other hand, the map from NIROrange (middle panel) is sensitive to the vegetated area as well because of the vegetation’s red edge. The bright regions in the NIROrange trace the vegetated area of Africa, South America, and Eurasia well.

We emphasize the fact that the SOT maps the albedo contrast over planetary surface and the map retrieved by the SOT is independent of any assumptions of surface components. The interpretation of the albedo contrast is in our hands and depends on our knowledge about reflection spectra of materials. Hence the degeneracy among components with similar reflection spectra is inevitable. For instance, distinguishing clouds from snow is difficult from the bands we considered because they both can have high and flat reflection spectra in the optical range.

Panel (e) displays the same results as panel (d) except that imposed S/N per data point is 100. The maps with S/N data have higher spatial resolution than those with S/N data and the continental distribution are clearly recovered.

In short summary, the SOT can also be applied to Earth-twins in inclined orbits with observations at the phases larger than . These phases are also favorable from the view point of observation because the intensity of planet is larger at smaller phase angle.

4.3. Measurements of Obliquity and Equinox

Figure 7.— Estimated value of obliquity, (top row), and simultaneously estimated value of equinox phase, (bottom row), as a function of input value of obliquity, . Different colors show the result with varying input value of equinox phase, . Three inclinations, (left column), (center column), and (right column), are studied. Signal-to-noise ratio is 20.

In this section, we test the obliquity measurement introduced in Section 2.4. Unlike Section 4.2, we regard not only but also and as fitting parameters (but the orbital inclination remains to be fixed). We create mock light curves of Earth-twins with varying , , and perform the SOT on them respectively. The same data sets for cloud cover, surface reflectivity, and snow cover are adopted as input in the simulations though the planetary climates may be admittedly inconsistent with obliquity when .

Figure 7 shows the measurement of and from the NIRBlue mapping as a function of input value of . As can be seen, the estimated values are well correlated to the input values. For the highly oblique planets (), the estimation works remarkably well, consistent with the face-on cases discussed in Paper I. Depending on the geometry of observation, there are systematic biases in some cases. Identifying the exact cause of the systematics is complex and we postpone the problem of precise measurement of obliquity until a future paper. Nevertheless, the fact that the SOT can roughly measure the obliquity and equinox is remarkable since the obliquity measurement is significantly difficult with other method.

We also tried the same estimation with single-band light mapping, but the result is much less reliable. Not only latitudinal but also longitudinal inhomogeneity just like the land distribution of the Earth appears to be essential for measuring the planetary spin axis.

4.4. Prospects with Future Instruments

Symbol Quantity Value Unit
Exposure time 24/303600 s
Folded days 6 days
Sharpness 0.08
End-to-end efficiency 0.5
Dark rate 0.001 counts s
Read noise 2 read
QE Quantum efficiency 0.91
Zodiacal light in magnitude 23 mag arcsec
Zodiacal light 2
Table 2Observation Parameters for Noise Estimation
Figure 8.— Signal-to-noise ratio as a function of aperture size on the assumption of parameters listed in Table 2.

Since the development for direct imaging observation is currently ongoing and not yet put into practice, we estimate the required aperture size of space telescopes for the SOT simply as follows. We consider readout noise, dark noise, and noise from exozodiacal light in addition to the photon noise of planetary light itself, and calculate S/N as a function of aperture size with Equations (21)(24) of Kawahara & Fujii (2010). The leakage of the light from the host star might be a dominant noise source for some instruments but depends on the situation. Considering that the occulter concept can safely neglect the leakage of the host star, we do not include it in our S/N estimation. The input values for the parameters to specify the instrumental property are summarized in Table 2.

Based on these assumptions, S/N of 20 with an Earth-sized planet at the distance of 10 pc will require 5 m aperture, and that of 100 will be equivalent to 15 m aperture size (Figure 8). This is significantly smaller than the scale proposed for future missions aiming to obtain spatially-resolved image of exoplanets in a direct manner (e.g. hypertelescope; Labeyrie, 1999). This is the outcome of the SOT which does not resolve the image of the planet instantaneously but instead solves the inverse problem of annual light curves.

5. Case Studies

In this section, we apply the SOT to different configurations. We consider different geometries with different combinations of orbital inclination and obliquity (Section 5.1), and different continental distribution based on ancient Earth maps (Section 5.2).

5.1. Different Geometry

5.1.1 Upright Earth-Twin

Figure 9.— Two-dimensional mapping from NIRBlue (left) and NIROrange (right) light curves in the case of . Data at are used in the analysis.

For the face-on case examined in Paper I, the two-dimensional image can be retrieved for non-zero obliquity. For inclined () cases, however, even non-oblique planets (i.e., ) can provide two-dimensional information of the planetary surface in principle because the visible region on the surface changes according to the orbital revolution as shown in Figure 2(c). To demonstrate such cases, we simulate light curves of a mock Earth with zero obliquity with orbital inclination (Fig. 2(c)). In this case, we may ignore the seasonal variation of polar snow covers and surface BRDFs (we adopt the BRDFs of March all year around) but the same input data of clouds are used.

Figure 9 displays the resultant maps from NIRBlue and NIROrange for and . They shows as good correspondence to the real surface inhomogeneity although the resolution is slightly worse than the case discussed in Section 4.

5.1.2 Tidally Locked Earth-Twin

Figure 10.— Left: albedo mapping from NIR band light curve of a hypothetical tidally locked atmosphere-less Earth. Number of sampling per orbit (with constant intervals) are 200 and all data are used in the analysis (). S/N is imposed. The white dotted lines are boundaries of illuminated and visible area. Right: albedo mapping from NIRBlue light curve of a hypothetical tidally-locked planet with Earth-twin surface and cloud cover. Number of sampling per orbit (with constant intervals) are 400 and data at are used in the analysis (). S/N is imposed. Assuming S/N does not improve the maps significantly.

For future projects of direct imaging of HZ planets, those around late-type stars are potentially favorable targets since late-type stars are most abundant in the universe. Since HZs around the late-type stars are so close to the star, planets there are likely tidally locked (Kasting et al., 1993). While the illuminated area does not change for the tidally-locked planet, the visible area does change unless a planet is in a face-on orbit (Figure 2(c)), which makes the SOT applicable to such cases as well, though the reconstructable area is inevitably limited to the illuminated half of the surface (see Appendix B).

Here we consider two tidally locked cases (i.e., ): one is an atmosphere-less Earth to see the results for the idealized case. Another is an Earth-twin with cloud cover same as the previous section. In both cases, again, we ignore seasonal variation of snow cover and surface albedos. Although those inputs may be unrealistic for actual tidally locked Earth analogs, we avoid getting deeply involved since the aim of this section is to see the basic applicability of the SOT to the tidally locked situation.

We assume that the sampling frequency is 200/orbit and 400/orbit, respectively, taking into account the fact that the orbit of tidally locked planets is much shorter than Earth’s (for instance, the orbital period of an Earth-mass planet around a star with is 2400 hr). Accordingly, we reduce the number of pixels to ( in order to avoid underdetermined problems.

The left column of Figure 10 displays the recovered maps from the atmosphere-less light curve at the NIR band with S/N in the case of . For demonstrative purpose, we show the maps with , and . This inversion is in a sense straightforward since our simulation assumes that the surface is Lambertian. As expected, the continental distribution is recovered fairly well.

The right column shows the results of mapping from the NIRBlue light curves of a tidally locked Earth-twin with cloud cover at . Though the resolution is now highly deteriorated, there is still inhomogeneity which certainly traces from the real features. In this case, increasing S/N to 100 does not improve the image quality significantly.

5.2. Different Land Distribution

In this section, we consider the sensitivity of reconstruction to the continental distribution. One of the ways to do this is to use the continental distribution of ancient Earth (see Sanromá & Pallé, 2012). We adopt the Earth map in the latest Jurassic period (150 Ma) and that in the early Cambrian period (540 Ma) taken from the paleomap project777, 2001) shown by the left column of Figure 11. We simply assume that everywhere on the continents has the same soil-like spectrum with reflectivity 0.1 in Blue band and 0.3 in NIR band. Since the primary features of cloud pattern are controlled by the global circulation which depends largely on the spin rotation rate and planetary radius, one may expect that global cloud pattern is not completely different in these periods. Therefore we simply overlaid clouds of the present Earth on the ancient continental configuration.

The middle column of Figure 11 displays two-dimensional maps retrieved by the same procedure from the NIRBlue light curves with S/N. The resultant map of the Jurassic Earth shows the primary continental mass concentrated on the center and a smaller feature at the northeast, which trace the input map with lower resolution. Likewise, the map of Cambrian Earth also recovers the input feature. The difference among the resultant maps of Earths at different epochs (Figures 6 and 11) are evident. These results imply that the SOT is applicable to various continental configurations.

Figure 11.— Left: input data of land distribution based on paleomap of the latest Jurassic period. Right: reconstructed map from NIRBlue light curves with signal-to-noise ratio 20 based on the latest Jurassic (upper panels) and early Cambrian (lower panels) land distribution. Geometric parameters are . Data at are used in the analysis. The integrated sensitivity is plotted for each case in the right column.

6. Summary and Discussion

In this paper, we extended the SOT, a method to recover two-dimensional information of exoplanetary surface from annual light curves proposed in Paper I, to be applicable to various geometric cases. In Paper I, we considered a planet in a face-on orbit, for which global maps with highly oblique planets are available. In this paper, we extend the applicable range of this method to Earth-twins in an arbitrary inclined orbit, which, in some geometric cases, enables the sketch of almost global surface map of low-obliquity planet like the Earth. Generally, the effect of anisotropic scattering by cloud can cause bias, but we have shown that omitting the data at crescent phases allows us to reasonably recover the surface inhomogeneity of Earth-twins. In the same way as Paper I, an albedo mapping of an Earth-twin from single-band photometry primarily contains the information of the distribution of cloud and snow. Information of bare surface under the clouds can be estimated by mapping the difference of light curves of two different photometric bands, thanks to the fact that cloud reflection spectra are roughly independent of wavelengths. On the simple assumptions for observational noise described in Section 4.4, a space telescope with aperture size 15 m will allow as to obtain two-dimensional map of an Earth-twin at a distance of 10 pc with satisfactory spatial resolution from continuous 1/2 year observation.

We tested the reliability of the simultaneous obliquity measurement with mock light curves of Earth-twins with varying directions of the spin axis and the observer. We confirmed that the global inhomogeneity of surface albedo allows us to reasonably estimate the planetary obliquity.

In the sense that primary features in scattered light of the Earth come from clouds because of their global cover (typically more than a half of the entire surface is covered by clouds with optical thickness larger than 5), clouds can be regarded as a substantial noise source if we aim to know the surface components (Ford et al., 2001; Pallé et al., 2008; Oakley & Cash, 2009; Cowan et al., 2009; Fujii et al., 2010, 2011). It is even annoying since cloud pattern is variable in the timescale of hour and there is also a seasonal variation. Indeed, the short-time variation of cloud cover causes 10% fluctuation on diurnal light curves, which is comparable to the observational noise we typically imposed ( 12%). Nevertheless, it is known that the annual mean of cloud cover has a characteristic pattern as shown in the bottom panel of Figure 6, which is related to the global atmospheric circulation of the Earth. Reconstructed albedo maps are likely to infer such averaged features of cloud cover since the SOT utilizes long-term light curves. Therefore we are no longer bothered by temporarily localized cloud cover.

The SOT has a potential to reveal the presence and rough distributions of the continents. For life on the Earth, the presence of continents has played an essential role in evolution of life as a source of nutrients by weathering. Nutrients such as phosphorus are indispensable to microbes in ocean such as cyanobacteria, who has been creating oxygen in Earth’s atmosphere. Therefore, discovery of continents on HZ planets may encourage us to search for life there. In addition, planets with continents can be advantageous targets for the detection of direct indicators of photosynthesis, because the signature of chlorophyll of microbes in ocean is more difficult to detect than the red edge of land plants (Knacke, 2003). The distribution of continents, if exist at all, and the planetary obliquity can greatly affect the clime of the planet (e.g., Williams & Kasting, 1997; Abe et al., 2005). Knowing these properties via the SOT will potentially enable the synthetic discussion of habitability from a variety of aspects.

We are thankful to the anonymous referee for many constructive comments. We deeply thank Y. Suto, A. Taruya, E. L. Turner, and C. Pichon for useful discussions and comments. Y.F. thanks the members of the Department of Astrophysical Sciences of Princeton University for their kind hospitality and suggestive discussions during her visit. Y.F. gratefully acknowledges support from the Global Collaborative Research Fund (GCRF) “A World-wide Investigation of Other Worlds” grant and the Global Scholars Program of Princeton University. Y.F. and H.K are supported by JSPS (Japan Society for the Promotion of Science) Fellowship for Research, DC:23-6070 and PD:22-5467, respectively. We are grateful to OpenCLASTR project for providing the rstar6b package.


  • Abe et al. (2005) Abe, Y., Numaguti, A., Komatsu, G., & Kobayashi, Y. 2005, Icarus, 178, 27
  • Angel et al. (1986) Angel, J. R. P., Cheng, A. Y. S., & Woolf, N. J. 1986, Nature, 322, 341
  • Arnold et al. (2009) Arnold, L., Bréon, F., & Brewer, S. 2009, Int. J. of Astrobiol., 8, 81
  • Arnold et al. (2002) Arnold, L., Gillet, S., Lardière, O., Riaud, P., & Schneider, J. 2002, A&A, 392, 231
  • Borucki et al. (2011) Borucki, W. J., et al. 2011, ApJ, 736, 19
  • Cassan et al. (2012) Cassan, A., et al. 2012, Nature, 481, 167
  • Cowan et al. (2009) Cowan, N. B., et al. 2009, ApJ, 700, 915
  • Cowan et al. (2011) Cowan, N. B., et al. 2011, ApJ, 731, 76
  • Doughty & Wolf (2010) Doughty, C. E., & Wolf, A. 2010, Astrobiology, 10, 869
  • Ford et al. (2001) Ford, E. B., Seager, S., & Turner, E. L. 2001, Nature, 412, 885
  • Fujii et al. (2010) Fujii, Y., Kawahara, H., Suto, Y., Taruya, A., Fukuda, S., Nakajima, T., & Turner, E. L. 2010, ApJ, 715, 866
  • Fujii et al. (2011) Fujii, Y., Kawahara, H., Suto, Y., Fukuda, S., Nakajima, T., Livengood, T. A., & Turner, E. L. 2011, ApJ, 738, 184
  • Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., Wandelt, B. D., Hansen, F. K., Reinecke, M., & Bartelmann, M. 2005, ApJ, 622, 759
  • Hamdani et al. (2006) Hamdani, S., et al. 2006, A&A, 460, 617
  • Hansen (2010) Hansen, P. C. 2010, Discrete Inverse Problems: Insight and Algorithms (Philadelphia, PA: SIAM)
  • Howard et al. (2010) Howard, A. W., et al. 2010, Science, 330, 653
  • Kaltenegger et al. (2010) Kaltenegger, L., et al. 2010, Astrobiology, 10, 89
  • Kasting et al. (1993) Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108
  • Kawahara & Fujii (2010) Kawahara, H., & Fujii, Y. 2010, ApJ, 720, 1333
  • Kawahara & Fujii (2011) Kawahara, H., & Fujii, Y. 2011, ApJ, 739, L62 (Paper I)
  • Kiang et al. (2007a) Kiang, N. Y., et al. 2007a, Astrobiology, 7, 252
  • Kiang et al. (2007b) Kiang, N. Y., Siefert, J., Govindjee, & Blankenship, R. E. 2007b, Astrobiology, 7, 222
  • Knacke (2003) Knacke, R. F. 2003, Astrobiology, 3, 531
  • Labeyrie (1999) Labeyrie, A. 1999, in ASP Conf. Ser. 194, Working on the Fringe: Optical and IR Interferometry from Ground and Space, ed. S. Unwin & R. Stachnik (San Francisco, CA: ASP), 350
  • Leger et al. (1993) Leger, A., Pirre, M., & Marceau, F. J. 1993, A&A, 277, 309
  • Lucht et al. (2000) Lucht, W., Schaaf, C. B., & Strahler, A. H. 2000, IEEE Trans. Geosci. Remote Sensing, 38, 977
  • Menke (1989) Menke, W. (ed.) 1989, Geophysical Data Analysis: Discrete Inverse Theory (New York: Academic)
  • Montañés-Rodríguez et al. (2006) Montañés-Rodríguez, P., Pallé, E., Goode, P. R., & Martín-Torres, F. J. 2006, ApJ, 651, 544
  • Nakajima (1983) Nakajima, T. 1983, J. Quant. Spectrosc. Radiat. Transfer, 29, 521
  • Nakajima & Tanaka (1988) Nakajima, T., & Tanaka, M. 1988, J. Quant. Spectrosc. Radiat. Transfer, 40, 51
  • Oakley & Cash (2009) Oakley, P. H. H., & Cash, W. 2009, ApJ, 700, 1428
  • Owen (1980) Owen, T. 1980, Strategies for the Search for Life in the Universe, ed. M. D. Papagiannis (Astrophysics and Space Science Library, Vol. 83; Dordrecht: Reidel), 177
  • Pallé et al. (2008) Pallé, E., Ford, E. B., Seager, S., Montañés-Rodríguez, P., & Vazquez, M. 2008, ApJ, 676, 1319
  • Qiu et al. (2003) Qiu, J., et al. 2003, J. Geophys. Res. (Atmospheres), 108, 4709
  • Robinson et al. (2010) Robinson, T. D., Meadows, V. S., & Crisp, D. 2010, ApJ, 721, L67
  • Russell (1906) Russell, H. N. 1906, ApJ, 24, 1
  • Sanromá & Pallé (2012) Sanromá, E., & Pallé, E. 2012, ApJ, 744, 188
  • Scotese (2001) Scotese, C. R. 2001, Digital Paleogeographic Map Archive on CD-ROM, PALEOMAP Project, Arlington, TX
  • Seager et al. (2005) Seager, S., Turner, E. L., Schafer, J., & Ford, E. B. 2005, Astrobiology, 5, 372
  • Selsis et al. (2002) Selsis, F., Despois, D., & Parisot, J.-P. 2002, A&A, 388, 985
  • Tarantola (2005) Tarantola, A. 2005, Inverse Problem Theory and Methods for Model Parameter Estimation (Philadelphia, PA: SIAM)
  • Tinetti et al. (2006a) Tinetti, G., Meadows, V. S., Crisp, D., Fong, W., Fishbein, E., Turnbull, M., & Bibring, J. 2006a, Astrobiology, 6, 34
  • Tinetti et al. (2006b) Tinetti, G., Meadows, V. S., Crisp, D., Kiang, N. Y., Kahn, B. H., Fishbein, E., Velusamy, T., & Turnbull, M. 2006b, Astrobiology, 6, 881
  • Williams & Kasting (1997) Williams, D. M., & Kasting, J. F. 1997, Icarus, 129, 254
  • Woolf et al. (2002) Woolf, N. J., Smith, P. S., Traub, W. A., & Jucks, K. W. 2002, ApJ, 574, 430
  • Zhdanov (2002) Zhdanov, M. S. 2002, Geophysical Inverse Theory and Regularization Problems, Methods in Geochemistry and Geophysics (New York: Elsevier)

Appendix A A. Weight Function

Here we present the explicit form of the weight function in Equation (9),


The boundary curves of the visible and illuminated area are expressed as




Using the spherical coordinate to the observer on the planetary surface, one can rewrite Equation (A4) to


We note that and does not depend on .

Converting to the spherical surface at (i.e. ) , we obtain the transformation formula to as


Appendix B B. Reconstructable Area

Figure 12.— Contour plot of the cover fraction of the illuminated and visible surface over a planetary year as a function of and . Solid contours show and . Each panel displays different obliquity corresponding to , and .

Figure 12 shows the contour of the area cover fraction of the surface sampled at least once during a planetary year , that is the maximum reconstructable area, for various geometric configurations. We also derive the cumulative probability distribution function . The whole region of the visible area is covered in a planetary year, we only have to consider the cover fraction of the visible area ( does not depend on the obliquity). Considering diurnal rotation, we can express the cover fraction as a function of ,


Using the uniform distribution of the line of sight and Equation (B1), we obtain


The does not depend on the obliquity and one can get above 95% (80%) of the surface information in 43.5% (80%) for the randomly selected exoplanets, in principle.

When a planet is in a synchronous orbit, a half of the planetary surface is always dark. We obtain


The CPDF of the tidally locked planet is expressed as


Both & are shown in Figure 13.

Figure 13.— Cumulative probability distribution function (CPDF) of the cover fraction on the assumption of the random line of sight (), . This function is independent of the obliquity . The probability is 43.5% and is 80%. The dashed line indicates the CPDF for tidally locked planets: . The probability is 43.5% and is 80%.

Appendix C C. Bayesian Interpretation of the Tikhonov Regularization

Following Tarantola (2005), we discuss the Bayesian interpretation of the Tikhonov regularization. The probability of a posterior information in the model space is given by


Here we neglect modelization uncertainties. Introducing a Gaussian model for the a priori information on the model parameters with the covariance matrix ,


and assuming the statistical noise is an independent Gaussian noise with the covariance matrix


one can obtain


The model which maximizes a posteriori probability is given by minimizing


If taking the independent Gaussian noise and , we obtain the Tikhonov regularization ,


Appendix D D. The L-curve Criterion

Figure 14.— dependence of the best-fit model. The left panel shows the L-curve for the mock observations with , and the additional 5% noise. The input data are a land/ocean binary map of the Earth, which has 1 on land and 0 on ocean. A solid curve indicates the parametric plot of and as a function of the regularization parameter . A dotted curve shows the curvature of the L-curve, . The definition of is written in Appendix E. We took three different regularization parameters: (A), , which corresponds to the maximum curvature (B), and (C). The recovered models are displayed in the right panels: the upper, middle, and lower panels correspond to A, B, and C.

The regularization parameter is chosen based on “the L-curve criterion” following Hansen (2010). The L-curve is a parametric plot of the model norm and the residuals as a parametric function of . In the left panel of Figure 14, we demonstrate the L-curve of the mock observations of an atmosphere-less Earth with and 3% noise. As increases, rapidly decreases until B (the L-curve’s corner), and then the slope becomes gentler. The large value of (small ) implies that the perturbation error due to observational noise is large, while large (small ), indicates that the bias due to the regularization (regularization error) is large. The L-curve corner serves as the point of good balance between these errors. The L-curve corner is defined by the maximum curvature in the log-log space, as described in Appendix E(Hansen, 2010). Adopting different values for to Equation (21), we see the dependence of the recovered models in the right panels. As shown in the left panel, the L-curve criterion gives at the maximum curvature point (B). Right panel of Figure 14 shows the recovered maps with (A), the L-curve corner (B), and (C). The recovered map for A exhibits large perturbation noise, while the resolution of the model is degraded for the case (C). The L-curve criterion balances these two opposite effects (B).

Appendix E E. Maximum curvature of the L-curve

Here, we briefly describe the derivation of the curvature of the L-curve following Hansen (2010). We first introduce and . The curvature of the L-curve is defined in the space,


where the prime indicates the derivative by . The derivatives of and are


Using Equations (E2) and (E2), one obtain

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description