GLOBAL MAPPING OF EARTH-LIKE EXOPLANETS FROM SCATTERED LIGHT CURVES

Recent progress in observational techniques has revealed various physical properties of exoplanets beyond orbital parameters and planetary mass. Detections of atmospheric components have been reported for several systems using spectroscopy at the planetary transit and secondary eclipse (e.g., Charbonneau et al. 2002; Vidal-Madjar et al. 2003, 2004; Tinetti et al. 2007; Swain et al. 2008, 2009). Interior compositions can be inferred from planetary mass and radius (e.g., Léger et al. 2009; Charbonneau et al. 2009). Constructions of thermal maps of the planetary atmosphere have been proposed by Williams et al. (2006) and Cowan & Agol (2008). A longitudinal thermal map of HD 189733b has been constructed by Knutson et al. (2007) based on the method proposed by Cowan & Agol (2008).

Nevertheless, an identification of planetary surface components still remains an ambitious challenge. One of the promising approaches is to use the scattered light of exoplanets through the direct imaging observation (e.g., Seager et al. 2000; Ford et al. 2001; Sudarsky et al. 2005). Ford et al. (2001) focused on the inhomogeneity of the Earth's surface which causes diurnal variation of the scattered light. They computed the scattered light from a model Earth observed at a distance of 10 pc and showed that time variations of the scattered light in different photometric bands highly depend on the geological and biological features on the planetary surface such as ocean, land, and even vegetation. More detailed characterizations (including spectroscopy) of the scattered light of the Earth and its time variations are discussed via both Earth-shine observations (e.g., Woolf et al. 2002; Arnold et al. 2002; Montañés-Rodríguez et al. 2006) and simulations (Tinetti et al. 2006a, 2006b; Montañés-Rodríguez et al. 2006). These studies have suggested a future possibility of investigating the surface of Earth-like exoplanets by scattered light curves. We note that such time variations of scattered light are also applicable to determine the rotation period as shown by Pallé et al. (2008).

A variety of inversion techniques of the planetary surface from the scattered light curves have been proposed. Surprisingly, the first theoretical study of scattered light curves to make albedo maps was performed at the beginning of the twentieth century, although the authors assume asteroids and satellites for the target (Russell 1906). Cowan et al. (2009) performed principal component analysis (PCA) on multi-band photometric data of the Earth observed by the Extrasolar Planet Observation and Deep Impact Extended Investigation (EPOXI) mission and extracted spectral features which roughly correspond to land and ocean. They also checked the time variation of these components and translated it to the longitudinal distribution of these components based on the formulation by Cowan & Agol (2008). Oakley & Cash (2009) paid attention to the gap of reflectivity between ocean and land and reproduced a longitudinal map of land. Fujii et al. (2010, hereafter F10) have developed a methodology to estimate the areas of ocean, soil, vegetation, and snow from multi-band photometry and showed that the area of these components can be recovered from mock observations of a cloudless Earth.

Since these authors focused on diurnal variations in mapping the surface, the resultant maps have only longitudinal resolution with little information on the latitudinal distribution. One of the goals of the present paper is to develop a method to map inhomogeneous surfaces of exoplanets with both longitudinal and latitudinal resolutions using both diurnal and annual variations of scattered light.

We also consider the determination of the planetary obliquity from time variation of planetary light. The obliquity is an important property of Earth-like planets with its strong implications for climate, habitability (e.g., Williams & Kasting 1997; Williams & Pollard 2003), and planetary formation. N-body simulations of the final stage of terrestrial planet formation indicated that the distribution of the obliquity ζ is isotropic (Agnor et al. 1999; Chambers 2001; Kokubo & Ida 2007). Gaidos & Williams (2004) modeled the infrared light curves of exoplanets and showed how the obliquity affects an annual variation. Oakley & Cash (2009) also pointed out a possibility of determining a planet's obliquity. In this paper, we demonstrate that the obliquity is estimated simultaneously with the global map of the planet by analyzing the scattered light curves over its orbital period.

The rest of the paper is organized as follows. We first review the estimation method of the weighted area from multi-band photometry proposed by F10 and describe our methodology to reconstruct the planetary surface and the obliquity measurement in Section 2. Assuming a future satellite mission for the direct imaging of Earth-like planets, we apply our methodology to mock observations based on real data of the scattering properties of the Earth in Section 3. Finally, we summarize our results in Section 4.

Let us briefly summarize the method of reconstructing a planetary surface from scattered light curves by F10. The scattered light from the planetary surface is characterized by the bidirectional reflectance distribution function (BRDF), f(ϑ₀, ϑ₁, φ; λ)[str^-1], where ϑ₀, ϑ₁, and φ are the incident zenith angle, the scattering zenith angle, and the angle between the incident and scattered light, respectively (Figure 1(a)). The BRDF is also a function of the wavelength λ. Using the BRDF at (ϕ, θ) fixed on the sphere (Figure 1(b)), the total scattered intensity at a given phase is provided by

$$\begin{eqnarray} &&I(\lambda) = F_\ast (\lambda) R_p^2 \int _{S_\mathrm{VI}} f(\phi,\theta ; \lambda) \cos {\vartheta _0 (\phi, \theta) } \cos {\vartheta _1 (\phi, \theta)} d \Omega,\nonumber\\ \end{eqnarray} \tag{ 1 }$$

where F_*(λ) is the incident flux at wavelength λ, R_p is a planetary radius, S_VI is the surface area facing the observer and illuminated by the host star (we call it the visible and illuminated area), and dΩ = sin θdθdϕ. The position on the surface (ϕ, θ) fully specifies ϑ₀, ϑ₁, and φ, and thus we write f(ϕ, θ; λ) ≡ f(ϑ₀(ϕ, θ), ϑ₁(ϕ, θ), φ(ϕ, θ); λ).

Zoom In Zoom Out Reset image size

In order to make a fitting model for the multi-band photometric data, F10 classified the planetary surface into four types: ocean, snow, soil, and vegetation. Denoting the specific BRDF of each type k by f_k(ϑ₀,  ϑ₁,  φ), one can write the local BRDF at position (ϕ, θ) as a summation of the specific BRDFs of the kth-type surface weighted by local cover fractions m^(k)(ϕ, θ):

$$\begin{eqnarray} f(\phi,\theta ; \lambda) &=& \sum _{k=1}^{N_C} f_k (\phi,\theta ; \lambda) m^{(k)} (\phi, \theta), \end{eqnarray} \tag{ 2 }$$

where N_C is the number of surface types and N_C = 4 in this case. Assuming that the scattering is isotropic (Lambertian), they approximate the BRDF by wavelength dependence of albedo a_k(λ),

$$\begin{eqnarray} &&f_k (\phi,\theta ; \lambda) = \frac{a_k (\lambda)}{\pi }. \end{eqnarray} \tag{ 3 }$$

Under the above assumption, Equation (1) reduces to

$$\begin{eqnarray} I(\lambda) &=& F_\ast (\lambda) R_p^2 \sum _{k=1}^{N_C} \frac{a_k (\lambda)}{\pi } \int _{S_\mathrm{VI}} m^{(k)} (\phi,\theta)\nonumber\\ &&\times \cos {\vartheta _0 (\phi, \theta)} \cos {\vartheta _1 (\phi, \theta)} d \Omega. \end{eqnarray} \tag{ 4 }$$

Assuming that they know the spectra for each component, F10 developed a reconstruction method for the weighted area of the kth component,

$$\begin{eqnarray} A_k &\equiv & R_p^2 \int _{S_\mathrm{VI}} W(\phi,\theta) m^{(k)} (\phi,\theta) d \Omega \nonumber \\ W(\phi, \theta) &\equiv & \frac{ \cos {\vartheta _0 (\phi,\theta)} \cos {\vartheta _1 (\phi,\theta) }}{ \int _{S_\mathrm{VI}} \cos {\vartheta _0 (\phi,\theta)} \cos {\vartheta _1 (\phi,\theta) } d \Omega }, \end{eqnarray} \tag{ 5 }$$

where W(ϕ, θ) is the weight function. The above definition reduces Equation (4) to a set of linear discrete equations:

$$\begin{eqnarray} C \frac{I(\lambda _b)}{F_\ast (\lambda _b)} &=& \sum _{k=1}^{N_C} a_k (\lambda _b) A_k, \nonumber \\ C &\equiv & \pi \int _{S_\mathrm{VI}} \cos {\vartheta _0 (\phi,\theta)} \cos {\vartheta _1 (\phi,\theta) } d \Omega, \end{eqnarray} \tag{ 6 }$$

where C is a factor that depends on the phase angle, i.e., the planetocentric angle between the host star and the observer, and λ_b is the center wavelength of the bth band. Using a cloudless model Earth, they showed that the weighted area can be approximately recovered by solving Equation (6) with an inequality condition A_k>0, and the summation of A_k(t) is approximately constant over time, ∑A_k(t) ≈ const. ≈ R²_p. It means that one can roughly estimate R_p as well by their method. Therefore, we use the fractional area A_k/R²_p normalized by R²_p in what follows. F10 have considered the spin rotation of a planet only and recovered the weighted area during one day of the planet and the weighted area can be converted to longitudinal maps by the inversion method proposed by Cowan & Agol (2008) and Cowan et al. (2009).

2.1. An Inverse Problem of a Two-dimensional Map

In the present paper, we consider the orbital motion to recover a two-dimensional world map of the planetary surface. We assume a model planet on a face-on circular orbit, which is the most promising case for planet mapping. We also assume that the orbital period is P_orb = 365 day and the spin rotational period is P_spin = 24 hr. We denote the phases of the orbital motion and spin rotation by Θ and Φ (Figure 1(c)), respectively. Our mock observation is assumed to be performed over one year from Θ = Θ_S to Θ = Θ_S + 2π. There are two unknown parameters in this geometry: the obliquity ζ and the initial orbital longitude Θ_S. We discuss the measurement of ζ and Θ_S in Section 2.3.

For our conventional geometry, the weight function is given by

$$\begin{eqnarray} &&W(\phi,\theta ;\Theta ;\Phi ;\zeta) = \frac{3}{2} \cos {\vartheta _0} \cos {\vartheta _1} = \frac{3}{2} [ \cos {\theta } \cos {\Theta } \sin {\zeta } \nonumber \\ &&\quad - \sin {\theta } \sin {(\phi +\Phi)} \sin {\Theta }- \sin {\theta } \cos {(\phi +\Phi)} \cos {\Theta } \cos {\zeta }]\nonumber\\ &&\quad\times \left[ \sin {\theta } \cos {(\phi +\Phi)} \sin {\zeta } +\cos {\theta } \cos {\zeta } \right]. \end{eqnarray} \tag{ 7 }$$

We provide the explicit form of the visible and illuminated area, S_VI in Appendix A. Note that the weight function is normalized to unity by integrating over the sphere:

$$\begin{eqnarray} &&\int W(\phi,\theta ;\Theta ;\Phi ;\zeta) d \Omega = 1. \end{eqnarray} \tag{ 8 }$$

Figure 2 shows the behavior of the weight function W(ϕ, θ; Θ; Φ; ζ) for cases of ζ = 0°, 45°, and 90°. The weight function covers the region of 0° < θ < 90° + ζ and 0° < ϕ < 360° as the planet rotates around the host star and its spin axis. As a result, the scattered light curves contain information about the land distribution up to 50(1 + sin ζ)% of the total surface area in principle.

Zoom In Zoom Out Reset image size

We consider how to deduce the local cover fraction m^(k)(ϕ, θ) if a set of the weighted area throughout the orbital and spin rotational periods is given. In order to solve this problem, we apply one of the techniques of tomography, the linear inverse problem, to the data. We assume that the planetary surface consists entirely of land and ocean and do not consider the other components nor further decompose them into clouds, vegetation, or soil in what follows. Therefore, we use the land cover fraction m(ϕ, θ) for the surface classification (k = land). The surface at (ϕ, θ) covered by land only (ocean only) is expressed by m(ϕ, θ) = 1 (m(ϕ, θ) = 0). We also denote the weighted area of land by A(Φ(t_i), Θ(t_i); ζ). In Section 3, we also consider the components of vegetation and soil using mock photometric data.

To begin with, let us discretize the observing time in N epochs and pixelize the planetary surface in M pixels. The weighted area recovered at the ith epoch t_i is written as

$$\begin{eqnarray} A(\Phi (t_i),\Theta (t_i);\zeta) &=& \sum _{j \,|\,(s_j \, \cap \, S_\mathrm{VI}) \ne \emptyset } \, \int _{s_j} W(\phi,\theta ;\Theta (t_i);\Phi (t_i);\zeta)\nonumber \\ &&\times\, m(\phi,\theta) d \Omega\nonumber \\ &=& \sum _{j \,|\,(s_j \, \cap \, S_\mathrm{VI}) \ne \emptyset } \, \langle m \rangle _{ij} \int _{s_j} W(\phi,\theta ;\Theta (t_i);\nonumber\\ &&\Phi (t_i);\zeta) d \Omega, \end{eqnarray} \tag{ 9 }$$

where s_j indicates the jth pixel on the surface. The weighted cover fraction 〈m〉_ij is defined as

$$\begin{eqnarray} &&\langle m \rangle _{ij} \equiv \frac{\int _{s_j}\! W(\phi,\theta ;\Theta (t_i);\Phi (t_i);\zeta) m(\phi,\theta) d \Omega }{\int _{s_j}\! W(\phi,\theta ;\Theta (t_i);\Phi (t_i);\zeta) d \Omega }. \end{eqnarray} \tag{ 10 }$$

Here, let us define the pixel-averaged cover fraction

$$\begin{eqnarray} &&\overline{m}_j \equiv \frac{\int _{s_j}\! m(\phi,\theta) d \Omega }{\int _{s_j}\! d \Omega }. \end{eqnarray} \tag{ 11 }$$

Under the assumption that the weighted cover fraction approximately equals the pixel-averaged cover fraction

$$\begin{eqnarray} &&\langle m \rangle _{ij} \approx \overline{m}_j, \end{eqnarray} \tag{ 12 }$$

we obtain

$$\begin{eqnarray} A(\Phi (t_i),\Theta (t_i);\zeta) &=& \sum _j \overline{m}_j \int _{s_j} W(\phi,\theta ;\Theta (t_i);\Phi (t_i);\zeta) d \Omega.\nonumber\\ \end{eqnarray} \tag{ 13 }$$

We introduce the data vector of weighted area ${\boldsymbol{a}}$, the design matrix G, and the model vector of cover fraction ${\boldsymbol{m}}$ as

$$\begin{eqnarray} {\boldsymbol{a}}&\equiv & \lbrace a_i = A(\Phi (t_i),\Theta (t_i);\zeta)/\sigma _i \, | \, i=1,2,...,N\rbrace, \nonumber \\ G &\,{\equiv}\, & \bigg\lbrace G_{ij} \,{=}\,\! \int _{s_j}\!\! W(\phi,\theta ;\Theta (t_i);\Phi (t_i);\zeta) d \Omega /\sigma _i \, |\, i\,{=}\,1,2,...,N,\nonumber\\ &&\qquad\ \ \ j=1,2,...,M\bigg\rbrace, \nonumber \\ {\boldsymbol{m}}&\equiv & \lbrace \overline{m}_j \,| \, j=1,2,...,M\rbrace, \end{eqnarray} \tag{ 14 }$$

where σ_i is the error of the weighted area. Then we can rewrite Equation (13) as

$$\begin{eqnarray} {\boldsymbol{a}}&=& G \, {\boldsymbol{m}}. \end{eqnarray} \tag{ 15 }$$

The above equation can be solved for ${\boldsymbol{m}}$ by minimizing χ²:

$$\begin{eqnarray} &&\chi ^2 \equiv |{\boldsymbol{a}}_\mathrm{data}- G \, {\boldsymbol{m}}|^2, \end{eqnarray} \tag{ 16 }$$

under inequality constraints:

$$\begin{eqnarray} &&0 \le m_{j} \le 1, \end{eqnarray} \tag{ 17 }$$

where ${\boldsymbol{a}}_\mathrm{data}$ is the observed weighted area with noise. We denote ${\boldsymbol{m}}$ that minimizes Equation (16) by ${\boldsymbol{m}}_{\mathrm{est}}$.

2.2. Solving the Inverse Problem

We are now in a position to solve Equation (16) with condition (17) using an idealized data set. Since one can choose an arbitrary division for the surface modeling, we use the following pixelization:

$$\begin{eqnarray} s_{j} &=& \Bigg\lbrace (\phi,\theta) | 1 - (1+\sin {\zeta }) \frac{j_\theta }{M_\theta } < \cos {\theta } \le 1 - (1+\sin {\zeta })\nonumber\\ &&\times \frac{j_\theta - 1}{M_\theta }, \frac{2 \pi }{M_\phi } (j_\phi -1) \le \phi < \frac{2 \pi }{M_\phi } j_\phi \Bigg \rbrace, \nonumber \\ j &=& (j_\phi - 1) M_\theta + j_\theta \mbox{ for } j_\phi = 1,2,...,M_\phi \nonumber\\ && {\rm and}\quad j_\theta = 1,2,...,M_\theta. \end{eqnarray} \tag{ 18 }$$

The above pixelization has a constant solid angle for a given ζ and the number of pixels does not change for different ζ's. We adopt M_ϕ = 9 and M_θ = 23 (see the right panel of Figure 3 for the ζ = 90° case) and assume the numbers of epochs for diurnal variations N_Φ = 23 and for annual variations N_Θ = 23. In Appendix B, we discuss the dependence of model degeneracies on the pixelization.

Zoom In Zoom Out Reset image size

We create the synthetic data of the weighted area using 1° × 1° fixed land/water masks of the ISLSCP II (input map, left panel of Figure 3). The pixel-averaged land cover fraction of this data (reference map) is shown in the right panel of Figure 3.

We fix the obliquity ζ = 90° and Θ(t₀) = Θ_S = 0° and compute the weighted area A_input(Φ, Θ; ζ = 90°) from the input map (not from the reference map). We compute the data A_data by adding a Gaussian noise N_g to A_input(Φ, Θ; ζ = 90°):

$$\begin{eqnarray} &&A_\mathrm{data} (\Phi, \Theta) = A_\mathrm{input} (\Phi, \Theta ; \zeta =90^\circ) + N_g. \end{eqnarray} \tag{ 19 }$$

We denote the standard deviation of the Gaussian random variable N_g by σ ≡ 〈N²_g〉^1/2.

The linear inverse problem (Equation (15)), or the equivalent least-square problem (Equation (16)), under the inequality constraints (Equation (17)) can be solved with the Bounded Variable Least-Squares Solver (BVLS) developed by Lawson & Hanson (1974, 1995).³ The BVLS is a generalization of the non-negative least square described in Lawson & Hanson (1974, 1995) and uses the QR decomposition (see also Menke 1989). The BVLS works for both overdetermined and underdetermined (ill-condition) problems. Even for ill-condition problems, the BVLS provides one of the (non-unique) solutions (Lawson & Hanson 1995).

By solving the least-square problem (Equation (16)) with the BVLS algorithm, we obtain the recovered map ${\boldsymbol{m}}_{\mathrm{est}}$. Figure 4 displays the two-dimensional images of the weighted area (left), the recovered maps (middle), and the prediction errors, $[A_\mathrm{data} (\Phi (t_i), \Theta (t_i)) - G {\boldsymbol{m}}_{\mathrm{est}}]/\sigma$ (right) for ζ_input = 90°. In this figure, we regard ζ and Θ_S as fixed parameters and use the input values ζ = ζ_input = 90° and Θ_S = Θ_S,input = 0° for Equation (14). We consider four cases of the noise, σ = 0, 0.01〈A_input〉, 0.1〈A_input〉, and 0.3〈A_input〉, where 〈A_input〉 is the average of A_input(ϕ, θ; ζ = 90°) and is approximately equal to 0.3 for the case of Earth. As discussed in Appendix C, we confirmed that the estimated map is the unique solution under the inequality constraint. Comparing the recovered maps with the reference map (right panel in Figure 3), we conclude that one can approximately recover the land distribution from the two-dimensional map of the weighted area.

Zoom In Zoom Out Reset image size

We also solve the inverse problem in the case of the obliquity of Earth, ζ_input = 2345 (Figure 5). Even for a small obliquity like Earth's, one can recover the land distribution of 50 [1 + sin(ζ = 2345)] ≈ 70 % of the planetary surface. We conclude that our methodology can recover main features of the continents on the planetary surface even if the planet has a small obliquity as is the case with Earth.

Zoom In Zoom Out Reset image size

2.3. Measurement of Obliquity

So far, we have assumed that the obliquity ζ and the initial phase Θ_S in the orbit are known a priori. Next, we try to estimate the obliquity ζ by our method itself. Figure 6 indicates the obliquity dependence of images of the weighted area. We determine the best-fit value of ζ and Θ_S by minimizing χ² distance:

$$\begin{eqnarray} \chi ^2 (\zeta, \Theta _S) &\equiv & \big| {\boldsymbol{a}}- G(\zeta,\Theta _S) {\boldsymbol{m}}_{\mathrm{est}}^{\zeta,\Theta _S}\big|^2 \nonumber \\ G(\zeta,\Theta _S) &\equiv & \bigg\lbrace G_{ij} (\zeta,\Theta _S) = \int _{s_j} W(\phi,\theta ;\Theta (t_i);\Phi (t_i);\zeta) d \Omega /\sigma _i \nonumber \\ &\,& | \, i=1,2,...,N, j=1,2,...,M\bigg\rbrace \nonumber, \\ \Theta (t_0) &=& \Theta _S, \end{eqnarray} \tag{ 20 }$$

where ${\boldsymbol{m}}_{\mathrm{est}}^{\zeta,\Theta _S}$ is determined by minimizing χ² with given ζ and Θ_S using the BVLS fitting. Figure 7 displays χ²(ζ, Θ_S) for ζ_input = 45° and Θ_S,input = 180°. We use the amoeba routine (Press et al. 1992) to search for ζ and Θ_S, which minimizes χ²(ζ, Θ_S). We estimate errors of ζ and Θ_S by the bootstrap resampling because errors are originated not only from additional Gaussian noise but also from pixelization, linearization, and other systematics. After 10³ iterations of the bootstrap, we obtain ζ_est = 463 ± 06 and Θ_S,est = 1818 ±  08 for 10 % additional noise and ζ_est = 503 ± 23 and Θ_S,est = 1830 ± 23 for 30 % additional noise. The estimated values of both obliquity and initial phase in orbit agree remarkably well with the input values.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Now, we apply our methodology to a mock multi-band observation of a cloudless Earth as a more realistic demonstration. We simulate scattered light curves of a cloudless Earth in the same manner as the simulation in F10. We use the solar spectrum as an incident flux on the model Earth. We assume an observer seeing the planet on a face-on and circular orbit at a distance of 5 pc (Case A) or 10 pc (Case B) from it. The scattered light from land is computed using the BRDF of actual land surface assigned by the MODIS data set ("snow-free gap-filled MODIS BRDF Model Parameters"). We adopt the data of April and ignore the seasonal variation of the BRDF at each point on the planetary surface. The scattered light from ocean is calculated with the BRDF model for wavy ocean described in Nakajima & Tanaka (1983). In our model, the snow cover regions around the poles are replaced by ocean. We also include the effect of Rayleigh scattering by the atmosphere with a single scattering approximation between the atmosphere and the underlying surface. Neither the effect of clouds nor the molecular absorption is, however, incorporated.

Our assumed observation parameters are listed in Table 1. As an observing system, we assume a satellite telescope with 1.1 m aperture. This assumption comes from the basic architecture proposed for the occulting ozone observatory (O3) (Kasdin et al. 2010), which is a satellite mission for direct imaging of exoplanets in UV/optical/near-IR bands shading the light from the host star by a 30 m external occulter. The scattered light is computed at the central wavelength of each band and multiplied by its band width. We consider four bands centered at the wavelengths listed in Table 2, which correspond to the bands of the MODIS land BRDF data. We assume 0.1 μm as the bandwidth.

Table 1. Observation Parameters

Symbol	Quantity	Value	Unit
texp	Exposure time	24/23×3600	s
n	Folded days	14	days
D	Diameter of telescope aperture	1.1	m
Ψ	Sharpness	0.0433
α	Pixel scale	0.03125	arcsec pixel−1
h	End-to-end efficiency	0.5
υ	Dark rate	0.001	counts s−1
κ	Read noise	2	$\sqrt{{\rm counts}}$/read
QE	Quantum efficiency	0.91
Ωz	Zodiacal light in magnitude	23	mag arcsec2
Z	Zodiacal light	1
F0	Zero flux	1.4 × 104 (band 1)	cts cm-2 nm s−1
		9.6 × 103(band 2)	cts cm-2 nm s−1
		7.2 × 103(band 3)	cts cm-2 nm s−1
		4.1 × 103(band 4)	cts cm-2 nm s−1

Download table as:  ASCIITypeset image

We divide the orbital phase into 23 (Θ_i = 2πi/23,  i = 0, 1,..., 22) and the spin rotational phase into 23 (Φ_i = 2πi/23,  i = 0, 1,..., 22). We compute the scattered light curves in four bands with the exposure time t_exp = 24/23 hr throughout one year. Then, we stack data for 14 days centered at each orbital phase Θ_i and fold the light curves according to its spin rotation period to obtain diurnal curves at the orbital phase, assuming that the spin rotational period is already known through periodogram analysis (Pallé et al. 2008). Thus, the resultant light curves have 23 × 23 data points for each band.

As for the noise, we consider the photon shot noise, the exzodiacal light, dark noise, and read noise. We ignore the leakage of the light from the host star. Thus, the signal-to-noise ratio (S/N) is expressed as

$$\begin{equation} {\rm S/N} = \frac{S}{\sqrt{S+N_{\rm z}+N_{\rm d}+N_{\rm r}}}, \end{equation} \tag{ 21 }$$

where S, N_z, N_d, and N_r are the scattered light from the planet, and the contribution from zodiacal light, dark noise, and read noise, respectively. Specifically, they are given as

$$\begin{eqnarray} N_{\rm z} &\equiv & hZF_0 \alpha ^2 10^{-\frac{\Omega _z}{2.5}} \frac{1}{\Psi }\,\pi \left(\frac{D}{2} \right)^2 t_{\rm exp} n \Delta \lambda, \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} N_{\rm d} &\equiv & \frac{\upsilon t_{\rm exp} n}{\Psi }, \end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray} N_{\rm r} &\equiv & \frac{\kappa ^2 n}{\Psi }. \end{eqnarray} \tag{ 24 }$$

The parameters and their values we adopt are listed in Table 1 (see, e.g., Savransky et al. 2010). The zero flux F₀ for calibration is calculated by linearly interpolating the spectrum of Vega given by Tables 1 and 2 in Colina et al. (1996) and shown in Table 3 in this paper. In order to quantitatively assess each noise level, we list the amplitudes of these noises per epoch in Table 3 together with the total average of the signal S in the case of ζ = 90°. The averaged S/N per epoch is also listed in Table 3. We note that the averaged signal per epoch is comparable with the noise level.

Table 3. Assumed Signal and Noises in Four Bands

Band b	1	2	3	4
Averaged signal $\bar{S}$ (case A; 101 cts epoch−1)	40.9	40.3	36.9	40.9
Averaged signal $\bar{S}$ (case B; 101 cts epoch−1)	10.2	10.1	9.2	10.2
Exozodi Nz (103 cts epoch−1)	12.4	9.2	7.0	4.0
Dark noise Nd (103 cts epoch−1)	1.2
Read noise Nd (103 cts epoch−1)	1.3
Averaged S/N (case A; 1 epoch−1)a	3.3	3.6	3.7	4.9
Averaged S/N (case B; 1 epoch−1)a	0.83	0.93	0.94	1.3

Note. ^aThe averaged S/Ns are defined by $\bar{S}/\sqrt{\bar{S}+N_z+N_d+N_r}$.

Download table as:  ASCIITypeset image

The obtained mock data I(λ_b) at each epoch (Θ(t_i), Φ(t_i)) is decomposed into the weighted area of four surface types, (ocean, vegetation, soil, and snow) by solving Equation (6) following F10. The weighted area of land is computed by summing up that of soil, vegetation, and snow, that is, A(Θ, Φ) = A_soil(Θ, Φ) + A_vegetation(Θ, Φ) + A_snow(Θ, Φ). While F10 used nonlinear fitting in order to constrain A_k to be positive, we use the BVLS linear solver instead without upper constraints. The resultant A(Θ, Φ) as a function of Θ and Φ is displayed in the left panels of Figure 8.

Zoom In Zoom Out Reset image size

We now apply the reconstruction method described in Section 2 (see Equations (14), (16), and (17)) to the set of obtained A(Θ, Φ). For the errors σ_i in Equation (14), we first compute the variance of the weighted area by 100 times resampling of photon counts according to Poisson statistics, and we call such variance as σ²_i,res. The σ_i,res, however, turned out not to relevant for σ_i in Equation (14) because it often results in zero or very small errors when A_i = 0 due to the boundary condition (A_i ⩾ 0). Therefore, we substitute an average value of errors σ = ∑^N_i=1σ_i,res/N to Equation (14) instead. For the case that all estimated errors of A_k are positive, the results with σ_i = σ_i,res do not change significantly. In the fitting, we assume that the local cover fraction of each surface type remains unchanged and thus ignore the seasonal variations of vegetation, soil, and snow. The center and right rows of Figure 8 display recovered maps and prediction errors for ζ_input = 90° of the cases A and B. The major features of the surface can be seen in the resultant maps. The averaged prediction errors $|A_\mathrm{data} (\Phi (t_i), \Theta (t_i)) - G {\boldsymbol{m}}_{\mathrm{est}}|/N$ are 0.53 〈A_data〉 (case A) and 0.85 〈A_data〉 (case B), respectively.

The land and ocean distributions have been considered so far. Next, we create color composition maps of vegetation and soil on the land and ocean recovered map (Figure 9). Adopting the same method of the land recovery to the weighted areas of soil (A_soil(t_i)) and vegetation (A_vegetation(t_i)), we derive the soil and vegetation distributions and blend them by yellow (soil) and green (vegetation) over Figure 8, in other words, the land (white) and ocean (blue) distributions. Even using the weighted areas of soil and vegetation separately for the recovery, the resulting vegetation and soil maps almost distribute over the land recovered with the summation of the weighted areas (A_land(t_i) = A_soil(t_i) + A_vegetation(t_i) + A_snow(t_i)). In this figure, snow is ignored because the contribution of it is significantly small. The concentrations of vegetation at South America and soil at middle Africa seen in case A indeed correspond to the Amazon forest and the Sahara desert, respectively.

Zoom In Zoom Out Reset image size

Now, we try to measure the obliquity for the mock observation. Figures 10–13 display the light curves of four bands with Θ_S,input = 180°, corresponding to different obliquities, ζ_input = 90°, 60°, 45°, and 30°, respectively. The noises assumed in these figures are computed in case A. Performing the χ² fitting described in Section 2.3, we estimate the best-fit obliquities and the initial orbital longitude. The predicted curves with the best-fit values are shown by solid lines in Figures 10–13. The input and the best-fit values and the reduced χ² are listed in Table 4. We compute the χ² for the land distribution (Equation (20)). Therefore, the degree of freedom is 23 × 23 − 23 × 9 = 322. It is likely that the inequality constraints and the other systematics such as the Lambert assumption induce the relatively high reduced χ² values (∼2). One can see both diurnal and annual variations of light curves in these figures. The larger variations are seen in band 4 (0.8585 μm) due to large albedos of the land component (soil and vegetation; see Figure 7 in F10). There is "pinching" of the light curve at Θ_S = 0° (i ∼ 264) for ζ_input = π/4 because the maximum of the weight function is located at the planet's pole. In this period, the physical position of the maximum stays at the pole for the daily motion. As a result, there are little variations in the period. However, the "pinching" seen at Θ_S ∼ 180°(i ∼ 0) for ζ_input = 90° has a different origin. In this period, the weight function moves around the Southern Hemisphere, that is, ocean mainly.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Estimated Obliquity and the Initial Orbital Longitude by χ² Minimization and the EIC

Case	ζinput	ΘS,input	The χ2 Fit	The χ2 Fit	χ2	χ2/dof	EIC	EIC
			ζest	ΘS,est			ζest	ΘS,est
A	90°	180°	898 ± 12	1838 ± 22	650.3	2.02	890	1845
A	60°	180°	647 ± 48	1753 ± 35	661.2	2.05	653	1765
A	45°	180°	444 ± 33	1798 ± 41	731.1	2.27	445	1775
A	30°	180°	325 ± 30	1886 ± 56	677.5	2.10	307	1905
B	90°	180°	898 ± 35	1908 ± 175	596.1	1.85	881	1875
B	60°	180°	797 ± 81	1718 ± 248	568.5	1.77	485	1556
B	45°	180°	382 ± 216	2318 ± 436	609.5	1.89	267	1326
B	30°	180°	892 ± 171	2103 ± 327	714.4	2.23	79	1386

Download table as:  ASCIITypeset image

The χ² maps for the four input obliquities are displayed in Figure 14 and the best-fit obliquities and initial longitudes are listed in Table 4. For case A, the χ² fitting provides the estimated values close to the input ones within 10°. Figure 15 shows the difference of the predicted light curves with the best-fit obliquity ζ_est = 44.4 and that shifted ±33(∼1σ level estimated by the bootstrap resampling), ζ = 477 and 411.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The results for case B are also listed in Table 4. Figure 16 displays the χ² map for case B. The estimated values for ζ_input = 30° and 60° tend to deviate from the input one. As shown in the bottom right panel, the result for ζ_input = 30° is significantly biased. It is likely that the systematics affect the obliquity estimates for this noise level. We also perform the other method called the extended information criterion (EIC) in Appendix D. Although the tendency that the estimated value is biased to ζ = 90° vanishes, the best-fit value does not improve significantly.

Zoom In Zoom Out Reset image size

In short, the planetary obliquity can be estimated well for case A, while, for case B which has larger noises, it is marginal to derive reliable values of the estimated obliquity by our method.

We have developed a method of reconstructing the two-dimensional planetary surface via diurnal and annual variation of the scattered light. Applying the method to the mock photometric data, we have demonstrated that our method works for the mock Earth model, while this model has a lot of simplifying assumptions as follows: (1) cloudlessness, (2) a face-on circular orbit, (3) known reflectance spectra, (4) lack of atmospheric absorption, (5) known rotation rate, (6) a static map, and (7) the absence of a moon. We also found that the planetary obliquity can be estimated by this method. With our method, future satellite missions such as the occulting ozone observatory (Kasdin et al. 2010) might provide "a global map" of Earth-like exoplanets. While only the terrestrial planets have been considered in this paper, our method might be applicable to any planets with an inhomogeneous surface, including Jupiter-like exoplanets.

In this paper, we ignored the effect of clouds and expected that clouds affect the estimation like a statistical noise because of the relatively short time variation of clouds. The PCA performed by Cowan et al. (2009) is a promising approach because they could separate the land and ocean compositions even though they used the EPOXI data that contain the cloud effect. The effect of clouds is discussed elsewhere (Y. Fujii et al. 2010, in preparation). We also assumed a face-on circular orbit in this paper. This assumption might be too severe for practical applications. We will generalize our method in the next paper.

We are deeply grateful to Yasushi Suto, Atsushi Taruya, and Edwin L. Turner for helpful and insightful discussion. We are thankful to Dmitry Savransky, Jeremy N. Kasdin, and David N. Spergel, who enlightened us about the observing instruments designed for future satellite missions. We also thank the referee, Nicolas Cowan, for a lot of constructive comments. H.K. is supported by a Japan Society for Promotion of Science (JSPS) Grant-in-Aid for science fellows. This work is also supported by a Grant-in-Aid for Scientific Research from JSPS and from the Japanese Ministry of Education, Culture, Sports, Science, and Technology (20·10466 and 22·5467), and by the JSPS Core-to-Core Program "International Research Network for Dark Energy."

The visible and illuminated area is expressed as

$$\begin{eqnarray} &&S_\mathrm{VI}= \left\lbrace \begin{array}{@{}l@{}} \lbrace (\phi,\theta) | 0 \le \theta \le {\rm min}[\theta _V (\phi ;\Phi ;\zeta), \theta _{I,+} (\phi ;\Theta ; \Phi ;\zeta)] \rbrace \nonumber \\ \qquad\mbox{for $\quad 0 < \Theta < \pi /2$ or $3 \pi /2 \le \Theta < \pi $ }\!, \nonumber \\ \lbrace (\phi,\theta) | \theta _{I,-} (\phi ;\Theta ; \Phi ;\zeta) \le \theta \le \theta _{V} (\phi ;\Phi ;\zeta) \rbrace \nonumber \\ \qquad\mbox{for $ \quad\pi /2 \le \Theta < 3 \pi /2$}. \nonumber \end{array} \right.\\ \end{eqnarray} \tag{ A1 }$$

The boundary lines of visible area θ_V(ϕ; Φ; ζ) and illuminated area θ_I,+(ϕ; Θ; Φ; ζ) and θ_I,−(ϕ; Θ; Φ; ζ) are obtained by solving the following equation:

$$\begin{eqnarray} && W(\phi,\theta ;\Theta ;\Phi ;\zeta) = (\sin {\theta } \cos {(\phi +\Phi)} \sin {\zeta } +\cos {\theta } \cos {\zeta }) \nonumber \\ &&\quad \times (\cos {\theta } \cos {\Theta } \sin {\zeta } - \sin {\theta } \sin {(\phi +\Phi)} \sin {\Theta }\nonumber\\ &&\quad - \sin {\theta } \cos {(\phi +\Phi)} \cos {\Theta } \cos {\zeta }) = 0. \end{eqnarray} \tag{ A2 }$$

The explicit equations of θ_V(ϕ; Φ; ζ), θ_I,+(ϕ; Θ; Φ; ζ), and θ_I,−(ϕ; Θ; Φ; ζ) are expressed as

$$\begin{eqnarray} &&\theta _V (\phi ;\Phi ; \zeta) = \cos ^{-1}{\left(-\frac{\cos \!{(\phi +\Phi)} \sin {\zeta }}{\sqrt{\cos ^2{\zeta } + \cos ^2{(\phi +\Phi)} \sin ^2{\zeta }}} \right)},\nonumber\\ \end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray} \theta _{I,+} (\phi ;\Theta ; \Phi ;\zeta) &=& \cos ^{-1} (\eta (\phi +\Phi ;\Theta ;\zeta)) \end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray} \theta _{I,-} (\phi ;\Theta ; \Phi ;\zeta) &=& \cos ^{-1} (- \eta (\phi +\Phi ;\Theta ;\zeta)), \end{eqnarray} \tag{ A5 }$$

where

$${\fontsize{6.8}{8.9}\selectfont{\begin{eqnarray*} &&\hspace{-0.2pc}\eta (\phi ;\Theta ;\zeta)\nonumber\\ &&\equiv \big(\frac{\cos \zeta \cos \Theta \cos \phi + \sin \Theta \sin \phi }{\sqrt{\cos ^2 \zeta \cos ^2 \Theta \cos ^2 \phi + 2 \cos \zeta \sin \Theta \cos \Theta \sin \phi \cos \phi +\sin ^2 \zeta \cos ^2 \Theta +\sin ^2 \Theta \sin ^2 \phi }} \big). \end{eqnarray*}}}$$

We note here that the design matrix alone does not fully specify ${\boldsymbol{m}}_{\mathrm{est}}$ without any constraints. We perform singular value decomposition (SVD) of the design matrix:

$$\begin{eqnarray} &&G = U \Lambda V^T, \end{eqnarray} \tag{ B1 }$$

where U and V are N × N and M × M unitary matrices, and Λ is a N × M matrix:

$$\begin{eqnarray} \Lambda &=& \left(\begin{array}{@{}c@{}} \tilde{\Lambda }\\ 0 \end{array} \right), \end{eqnarray} \tag{ B2 }$$

$$\begin{eqnarray} \tilde{\Lambda }&\equiv & \mbox{diag}(l_1,..., l_M), \end{eqnarray} \tag{ B3 }$$

with l_j being the singular value of G in descending order. We denote the jth row of V by the vector ${\boldsymbol{v}_j}$.

In general, there are N_null zero components in the singular values, $(l_{M-N_\mathrm{null}+1}=l_{M-N_\mathrm{null}+2}=l_{M}=0)$. Corresponding ${\boldsymbol{v}}$ vectors are called null vectors. Linear combination of the null vectors does not affect the data vector,

$$\begin{eqnarray} &&G \sum _{s=1}^{N_\mathrm{null}} \beta ^{\mathrm{null},s} {\boldsymbol{v}}^{\mathrm{null},s} = {\boldsymbol{0}}, \nonumber \\ &&\quad{\boldsymbol{v}}^{\mathrm{null},s} \equiv {\boldsymbol{v}}_{M-N_\mathrm{null}+s}, \end{eqnarray} \tag{ B4 }$$

where β^null,s is an arbitrary coefficient. Equation (B4) indicates that the predicted model has a freedom to add any linear combination of the null vectors.

The solid curve in the left panel of Figure 17 displays the singular values in descending order (this plot is referred to as spectrum of data kernel) in the case of M_ϕ = N_Φ = N_Θ = 23 and N_θ = 9. Due to the truncation errors, the singular values obtained by numerical algorithms of the SVD do not vanish exactly. However, only three of them (j = M − 2, M − 1, and M) have significantly smaller values than other diagonal elements. We regard the elements smaller than 10⁻⁴ as zero and the corresponding row vectors of V as null vectors.

Zoom In Zoom Out Reset image size

Although the presence of null vectors does not directly lead to the indefiniteness if one adopts the inequality constraints (Equation (17)), a smaller number of null vectors makes it easier to interpret an uncertainty of the recovered map. Therefore, it is useful to know how the number of null vectors depends on pixelization. Then, one can decrease the number of the null vectors by an appropriate choice of pixelization. For simplicity, we fix M_ϕ = N_Φ = N_Θ and estimate the number of the null vectors by changing M_θ and M_ϕ. The left panel in Figure 17 shows the spectra of the singular value in descending order for three cases; M_ϕ = 22 (dotted), M_ϕ = 23 (solid), and M_ϕ = 24 (dashed) with M_θ = 9. The case of M_ϕ = 23 has a steeper spectrum than the others. It means that the choice of M_ϕ = 23 makes it easier to identify null vectors. The right panel of Figure 17 indicates the number of null vectors for different sets of M_ϕ and M_θ, where we define vectors with singular value below 10⁻⁴ as null vectors. There is a general tendency that the odd number bin of M_ϕ has lesser N_null. Although the selection of smaller values of M_θ or M_ϕ generally leads to a smaller number of the null vectors, it degrades the resolution of the map. As a compromise, we decided to adopt M_ϕ = 23 and M_θ = 9 for our fiducial values of the model in this paper. Figure 18 indicates the (unit) null vectors for M_ϕ = 23 and M_θ = 9, ${\boldsymbol{v}}^{\mathrm{null},1} = {\boldsymbol{v}}_{M-2}, {\boldsymbol{v}}^{\mathrm{null},2} = {\boldsymbol{v}}_{{M-1}}$, and ${\boldsymbol{v}}^{\mathrm{null},3} = {\boldsymbol{v}}_{M}$, where ${\boldsymbol{v}}_j$ is the jth row of V. These three null vectors for M_ϕ = 23 and M_θ = 9 affect the large-scale structure of the map but do not change the small-scale distribution like the continents of the present Earth.

Zoom In Zoom Out Reset image size

In this appendix, we examine whether there is the freedom for adding any components of null vectors to the recovered map derived by the BVLS. Because of the boundary condition for ${\boldsymbol{m}}_{\mathrm{est}}$ (Equation (17)), the estimated values of m_est j are classified into three types: (1) m_est j = 1 (on the upper boundary), (2) m_est j = 0 (on the lower boundary), and (3) 0<m_est j < 1 (in between). We focus on types (1) and (2) to consider the uniqueness of our recovered map. We denote the index j of type (1) by j− and that of type (2) by j+. In order not to violate the boundary condition, the linear combination of null vectors should satisfy the following constraints:

$$\begin{eqnarray} \sum _{k=1}^{3}\; \beta ^{\mathrm{null},k} v_{j-}^{\mathrm{null},k} &\le & 0 \mbox{\;\;\;for any $j-$}, \end{eqnarray} \tag{ C1 }$$

$$\begin{eqnarray} \sum _{k=1}^{3}\; \beta ^{\mathrm{null},k} v_{j+}^{\mathrm{null},k} &\ge & 0 \mbox{\;\;\;for any $j+$}. \end{eqnarray} \tag{ C2 }$$

These inequalities can be solved in terms of β^null,1 using the fact that v_j is positive (see Figure 18):

$$\begin{eqnarray} \beta ^{\mathrm{null},1} &\le & -\frac{v_{j-}^{\mathrm{null},2}}{v_{j-}^{\mathrm{null},1}} \beta ^{\mathrm{null},2} -\frac{v_{j-}^{\mathrm{null},3}}{v_{j-}^{\mathrm{null},1}} \beta ^{\mathrm{null},3}, \end{eqnarray} \tag{ C3 }$$

$$\begin{eqnarray} \beta ^{\mathrm{null},1} &\ge & -\frac{v_{j+}^{\mathrm{null},2}}{v_{j+}^{\mathrm{null},1}} \beta ^{\mathrm{null},2} -\frac{v_{j+}^{\mathrm{null},3}}{v_{j+}^{\mathrm{null},1}} \beta ^{\mathrm{null},3}. \end{eqnarray} \tag{ C4 }$$

Then, the following inequality is required so that β^null,1 exists

$$\begin{equation} \frac{v_{j-}^{\mathrm{null},2}}{v_{j-}^{\mathrm{null},1}} \beta ^{\mathrm{null},2} + \frac{v_{j-}^{\mathrm{null},3}}{v_{j-}^{\mathrm{null},1}} \beta ^{\mathrm{null},3} \le \frac{v_{j+}^{\mathrm{null},2}}{v_{j+}^{\mathrm{null},1}} \beta ^{\mathrm{null},2} +\frac{v_{j+}^{\mathrm{null},3}}{v_{j+}^{\mathrm{null},1}} \beta ^{\mathrm{null},3}. \end{equation} \tag{ C5 }$$

This equation can be reduced to the following expression:

$$\begin{equation} {\bf p}_{(j+, j-)} \cdot {\bf b} \ge 0, \end{equation} \tag{ C6 }$$

where p_(j+,j−) and b are two-dimensional vectors:

$$\begin{eqnarray} {\bf p}_{(j+, j-)} &\equiv & \left(\frac{v_{j+}^{\mathrm{null},2}}{v_{j+}^{\mathrm{null},1}} - \frac{v_{j-}^{\mathrm{null},2}}{v_{j-}^{\mathrm{null},1}}, \frac{v_{j+}^{\mathrm{null},3}}{v_{j+}^{\mathrm{null},1}} - \frac{v_{j-}^{\mathrm{null},3}}{v_{j-}^{\mathrm{null},1}} \right) \end{eqnarray} \tag{ C7 }$$

$$\begin{eqnarray} {\bf b } &\equiv & (\beta ^{\mathrm{null},2}, \beta ^{\mathrm{null},3}). \end{eqnarray} \tag{ C8 }$$

Now, let us denote the arguments of p and b by θ_α(j +, j−) and θ_β, respectively, i.e.,

$$\begin{equation} {\bf p }_{(j+, j-)} \equiv | {\bf p}_{(j+, j-)} | \;(\cos \theta _{\alpha }(j+, j-), \sin \theta _{\alpha }(j+, j-)), \end{equation} \tag{ C9 }$$

$$\begin{eqnarray} {\bf b } &\equiv & | {\bf b} | (\cos \theta _{\beta }, \sin \theta _{\beta }). \end{eqnarray} \tag{ C10 }$$

If |b| ≠ 0, the solution of Equation (C6) should satisfy

$$\begin{equation} \theta _{\alpha }(j+, j-) \le \theta _{\beta } \le \theta _{\alpha }(j+, j-). \end{equation} \tag{ C11 }$$

We calculate the above constraints for any combination of (j +, j−) and found that there is no θ_β that satisfies Equation (C11) for all combinations of (j +, j−). Thus, the only solution allowed is β^null,2 = β^null,3 = 0. From equations (C3) and (C4), it follows that β^null,1 = 0. As a result, we confirmed that there is no room for adding extra components of null vectors to the recovered map.

In this appendix, we try to perform the model selection of ζ_est and Θ_S,est using the EIC (e.g., Efron 1983; Konishi & Kitagawa 1996; Ishiguro et al. 1997) based on the Kullback–Leibler divergence. Although the model selection for ζ and Θ_S based on the χ² minimization is very convenient, the result might be biased for large noises compared with signal. The EIC is a bootstrap-based extension of Akaike Information Criterion defined as

$$\begin{eqnarray} &&\mbox{EIC} \equiv - 2 \sum _{i=1}^{N} \log { \mathcal {F}({\boldsymbol{a}}|{\boldsymbol{m}}_{\mathrm{est}}, \zeta _{\mathrm{est}}, \Theta _{S,\mathrm{est}})} + 2 b_B, \end{eqnarray} \tag{ D1 }$$

where $\mathcal {F}({\boldsymbol{a}}|{\boldsymbol{m}}_{\mathrm{est}}, \zeta _{\mathrm{est}}, \Theta _{S,\mathrm{est}})$ indicates the likelihood function. The best model or parameter set is given by minimizing the EIC. The bootstrap estimate of the bias b_B is given by Ishiguro et al. (1997):

$$\begin{eqnarray} b_B &=& \frac{1}{N_\mathrm{EIC}} \sum _{k=1}^{N_\mathrm{EIC}} \big[\log { \mathcal {F}\big({\boldsymbol{a}}^{\ast,k} |{\boldsymbol{m}}_{\mathrm{est}}^{\ast,k}, \zeta _{\mathrm{est}}^{\ast,k}, \Theta _{S,\mathrm{est}}^{\ast,k}\big)}\nonumber\\ && - \log { \mathcal {F}({\boldsymbol{a}}^{\ast,k} |{\boldsymbol{m}}_{\mathrm{est}}, \zeta _{\mathrm{est}}, \Theta _{S,\mathrm{est}})} + \log { \mathcal {F}({\boldsymbol{a}}|{\boldsymbol{m}}_{\mathrm{est}}, \zeta _{\mathrm{est}}, \Theta _{S,\mathrm{est}})}\nonumber\\ && - \log { \mathcal {F}\big({\boldsymbol{a}}|{\boldsymbol{m}}_{\mathrm{est}}^{\ast,k}, \zeta _{\mathrm{est}}^{\ast,k}, \Theta _{S,\mathrm{est}}^{\ast,k}\big)} \big], \end{eqnarray} \tag{ D2 }$$

where ${\boldsymbol{a}}^{\ast,k}, {\boldsymbol{m}}_{\mathrm{est}}^{\ast,k}, \zeta _{\mathrm{est}}^{\ast,k}$, and Θ*^,k_S,est indicate the kth bootstrap resample, ${\boldsymbol{m}}_{\mathrm{est}}$, ζ_est, and Θ_S,est, and N_EIC is the number of trials of the bootstrap. Assuming a Gaussian distribution for ${\boldsymbol{a}}$ with σ = 1, we derive

$$\begin{eqnarray} \mbox{EIC} &=& | {\boldsymbol{a}}- G(\zeta,\Theta _S) {\boldsymbol{m}}_{\mathrm{est}}|^2 + N\log {2 \pi } + b_B \nonumber \\ b_B &=& - \!\frac{1}{2 N_\mathrm{EIC}} \sum _{k=1}^{N_\mathrm{EIC}} \big[ \big| {\boldsymbol{a}}^{\ast,k} - G\big(\zeta ^{\ast,k},\Theta _S^{\ast,k}\big) {\boldsymbol{m}}_{\mathrm{est}}^{\ast,k} \big|^2 \nonumber\\ &&- \big| {\boldsymbol{a}}^{\ast,k} - G(\zeta,\Theta _S) {\boldsymbol{m}}_{\mathrm{est}}\big|^2 + \big| {\boldsymbol{a}}- G(\zeta,\Theta _S) {\boldsymbol{m}}_{\mathrm{est}}\big|^2\nonumber\\ && - \big| {\boldsymbol{a}}- G\big(\zeta ^{\ast,k},\Theta _S^{\ast,k}\big) {\boldsymbol{m}}_{\mathrm{est}}^{\ast,k} \big|^2 \big]. \end{eqnarray} \tag{ D3 }$$

We adopt N_EIC = 500 for bootstrap resampling. Figure 19 displays the dependence of the EIC on ζ_est and Θ_S,est. The minimum EIC estimate listed in Table 4 is derived by a stepwise search at 1° interval.

Zoom In Zoom Out Reset image size

The minimum EIC estimate listed in Table 4 is derived by a stepwise search at 1° interval. For case A, the EIC provides almost the same values as listed in Table 4. Figure 19 displays the EIC map for case B. Although the estimated obliquities for ζ_input = 30°, 45°, and 60° by the EIC do not agree well with the input value (Table 4), the bias of ζ_input = 60° and 30° seems to be improved compared to that from the χ² fitting.