# Mapping Earth Analogs from Photometric Variability:

Spin-Orbit Tomography for Planets in Inclined Orbits

###### Abstract

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 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 |
---|---|

(Coordinates) | |

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 |

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:

(1) | |||||

(2) | |||||

(3) | |||||

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

(4) |

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

(5) |

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

(6) |

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

(7) | |||||

(8) |

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:

(9) |

With this weight function, Equation (6) reduces to

(10) |

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

(11) | |||||

where is the phase angle, defined as

(12) |

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.

### 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

(14) |

where denotes the solid angle of th pixel.

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

(15) |

where is noise of the -th data point.

After these discretizations, we obtain a linear discrete equation

(16) | |||||

(17) | |||||

(18) |

As a pixelization scheme for the planetary surface, Hierarchical Equal Area isoLatitude Pixelization (HealPix) is adopted in this paper^{1}^{1}1http://healpix.jpl.nasa.gov/index.shtml.
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

(19) |

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:

(20) |

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

(21) | |||||

(22) |

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

(23) |

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:

(24) | |||||

(25) |

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 ^{2}^{2}2In 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.^{3}^{3}3U.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 Product^{4}^{4}4http://ladsweb.nascom.nasa.gov/. 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)^{5}^{5}5Available at http://modis.gsfc.nasa.gov/.
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.^{6}^{6}6http://modis-snow-ice.gsfc.nasa.gov/
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 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

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

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

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 |

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

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.

#### 5.1.2 Tidally Locked Earth-Twin

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 project^{7}^{7}7http://www.scotese.com/(Scotese, 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.

## 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.

## References

- 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),

(A1) | |||||

(A2) | |||||

(A3) | |||||

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

(A4) | |||||

and

(A5) | |||||

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

(A6) |

We note that and does not depend on .

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

(A7) | |||||

(A8) | |||||

(A9) |

## Appendix B B. Reconstructable Area

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 ,

(B1) |

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

(B2) | |||||

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

(B3) |

The CPDF of the tidally locked planet is expressed as

(B4) | |||||

Both & are shown in Figure 13.

## 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

(C1) |

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

(C2) |

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

(C3) |

one can obtain

(C4) |

The model which maximizes a posteriori probability is given by minimizing

(C5) |

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

(C6) |

## Appendix D D. The L-curve Criterion

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,

(E1) |

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

(E2) | |||||

(E3) | |||||

(E4) | |||||

(E5) |

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

(E6) |