Simultaneous characterization of the atmospheres, surfaces, and exomoons of nearby rocky exoplanets

1.   Introduction

Multiple space missions (Ricker et al., 2016; Rauer et al., 2014) and ground-based observation facilities (Pepe et al., 2010; Wheatley et al., 2013) are actively searching for potentially habitable exoplanets around nearby stars. The discoveries of Proxima b (Anglada-Escudé et al., 2016) and the TRAPPIST-1 system (Gillon et al., 2017) highlight the rapid progress in this research direction. Retrieving the composition of atmospheres is critical in understanding the habitability and/or biosignature of exoplanets. Because different surface types have different reflection spectra, a method to disentangle the observational effect of surface type from that of atmospheric species is important for interpreting reflection spectral data from exoplanet observations. The Moon is important for the long term stability of Earth's obliquity (Laskar et al., 1993) and climate (Williams et al., 1997; Forgan and Kipping, 2013). The ability to detect exomoons could also be important for our understanding of exoplanetary atmospheres (Rein et al., 2014; Li T et al., 2016). Despite many efforts (Kipping, 2009, Agol et al., 2015, Heller et al., 2016; Heller, 2017, and references therein), no exomoon has been detected up to now.

Inspired by hyperspectral data retrieval methods for Martian mineralogy analysis (the Full Constraint Least Squares algorithm (Heinz and Chang, 2001; Schmidt et al., 2014), referred to as the FCLS algorithm in the following), we developed a method to achieve the goals of obtaining atmospheric composition, retrieving disk-averaged fractions of surface types, and detecting exomoons simultaneously for nearby rocky exoplanet-exomoon systems.

2.   Methods

Following previous works (Rein et al., 2014; Li T et al., 2016), the combined reflection spectra of exoplanet-exomoon systems (REF) are calculated by weighting the reflection spectra of exoplanets and exomoons according to the radius ratio between exoplanets and exomoons (r_ratio=r_moon/r_planet):

where λ is the wavelength, T_planet and T_moon are the atmospheric transmissivity of the planet and moon respectively. For the purpose of this work, we assume airless exomoons (T_moon=1) because the timescale for exomoons at Earth-equivalent orbits to lose their volatiles should be small (Lammer et al., 2014). A_i and B_j are the fractions of surface types i/j on exoplanets and exomoons respectively (A_i ≥0, $\mathop \sum {A_i} = 1$ , B_j ≥0, $\mathop \sum {B_j} = 1$ ). REF_planet and REF_moon are the integrated surface reflectivity of exoplanet and exomoon respectively. REF_{surface, i} (Figure 1) are the reflectivity of 6 typical surface types for exoplanets (the USGS Digital Spectral Library) (Clark et al., 2007) and 5 for exomoons (APOLLO samples 15058 and 70017 in Lunar Soil Characterization Consortium) (LSCC, 2018). It should be noted that equations (1)–(3) apply to the enlightened areas of exoplanets and exomoons observable from Earth. Thus an underlying assumption is that the data analyzed should correspond to the same star-planet-moon configuration.

Figure  1.  Reflectivity of surface types on exoplanets (a) and exomoons (b) used in this work.

DownLoad: Full-Size Img PowerPoint

For all targets in this work, T_planet are calculated from the combination of 3 total pressure, 3 water vapor profiles (corresponding to cold, moderate, and warm climate states), 3 levels of CH₄ concentrations, and 17 CO₂ and 14 O₂ concentration levels (Table 1). Forty reflection spectra (REF) are produced for each T_planet, among which 20 for exoplanet-exomoon systems and 20 for single exoplanet systems, by randomly choosing A_i, B_j, and r_ratio (between 0 to 0.6). Although there are many other possible types of atmospheres and surfaces, the principle of our method should be demonstrated adequately by the chosen combinations of atmosphere and surface type.

Table  1.  Mixing ratios of atmospheric species

	Targets	Database
Water Vapor at Planet Surface	1.58×10–4, 9.25×10–4, 7.75×10–3	Identical to the left
CH4	10–5, 10–4, 10–3	Identical to the left
CO2	Evenly distributed between 10–4 and 1.6×10–1 in log10 space	1×10–1, 2×10–1, 3×10–1 and C1× ${{10}^{{-{C}}_{{2}}}}$, where C1=0, 1, 3, 7 and C2=2, 3, 4
O2	10–4 and values evenly distributed between 10–3 and 2.5×10–1 in log10 space	Identical to the above
Note: All cases include 3 atmospheric pressures (0.3, 1, and 3 bars).

 | Show Table

DownLoad: CSV

From the surface reflection spectra, the photon spectra of the hypothetical targets (S_target) are calculated at 40 bands uniformly spaced in frequency between 1 and 3 μm (R~20–60) assuming 1) a telescope with an effective area of 30 m² (similar to JWST); and 2) photon noise as the only significant noise. The latter assumption allows us to add random noise (e) with magnitude of $\sqrt {{{S}_{{\rm{target}}}}}$ to each photon spectrum. We emphasize that this work does not aim at accurately estimating observation time or providing constraints on instrumental noise for future observations. Instead we focus on studying the feasibility of the retrieval method. Thus our simplification on photon noise can be justified.

For Sun-like stars, we use a blackbody photon spectrum at 5800 K and scale it to Earth orbital distance (blue curve in Figure 2). We place all targets at 10-pc distance. As most potentially habitable exoplanets discovered already (Proxima b, TRAPPIST-1 planets, etc.) and to be discovered (TESS planets for example) are orbiting around M dwarfs, we also carry out retrieval analysis for hypothetical exoplanets-exomoons around M dwarfs by using the observed photon spectrum of GJ667C (France et al., 2016), scaled to the orbital distance of GJ667Cc (red in Figure 2). For this group of targets a distance of 6.8 pc is assumed. For both stars, we carry out analysis corresponding to observation time of 1 hour, 1 day, and 10 days. For the sake of simplicity, we assume in this study that the target exoplanets have radii of 1 Earth radius around Sun-like stars and the same size as GJ667Cc around M dwarfs. In reality, the actual sizes of exoplanets could be estimated from broadband photometry before reflection spectra are obtained, which would help to remove uncertainty in data retrieval. We emphasize here that the goal of this work is to retrieve planet properties from rocky exoplanets and exomoons. Thus our discussions do not include possible Earth-size exomoons around Jupiter-mass exoplanets.

Figure  2.  Incident spectra used in this work. The blue curve represents a blackbody spectrum from a Sun-like star with effective temperature of 5800 K scaled to Earth-Sun distance. The red curve is based on observed spectrum of M dwarf GJ667C scaled to the distance of exoplanet GJ667Cc.

DownLoad: Full-Size Img PowerPoint

The FCLS algorithm (Heinz and Chang, 2001) solves the following equations:

where S_obs is the noise-added photon spectrum; M is the photon spectra matrix database combining the transmission spectra of planetary atmospheres (T_planet), the reflection spectra of different surface types (REF_{surface, i}), the target-to-Earth distance, telescope size, and observation time. X is a 1-D array of surface type fractions:

here r_planet is between 0.8 and 1.5r_⊕ (the range of habitable planets) with an interval of 0.05r_⊕ and the term is added to check how uncertainty in planet size might influence the retrieval results. For each S_obs, the FCLS algorithm searches for X that minimizes the cost functions $C = \left\| {{{ S}_{{\rm{obs}}}} - { X} \cdot { M}} \right\|$ while satisfying positivity and sum-to-one constraint.

Before introducing the results, it is important to point out the difference between our method and those in previous exomoon detection efforts, most of which propose to analyze transit photometry or transit timing/duration variations (Heller et al., 2016; Heller, 2017, and references therein). Several works have been published on exomoon detection through analyzing phase-dependent thermal emission variability of exoplanet-exomoon system during a full exomoon orbit. Moskovitz et al. (2009) suggests that only Mars-size exomoons could be detected around Earth-size exoplanets based on surface temperature and thermal emission calculations. Robinson (2011) suggests that exomoons could outshine their host exoplanets at some wavelengths in thermal emission due to substantial atmospheric absorption of emissions from exoplanet surfaces. Gómez-Leal et al. (2012) suggests that because the lowest few km of the Earth are unseen in thermal emission, diurnal variations of exoplanets could be masked in thermal emission and thus the phase variations of exoplanet-exomoon systems could be dominated by exomoon signals. Forgan (2017) found that distinguishing the lunar component of the thermal emission phase curve of most exoplanet–exomoon systems will require a photometric precision of 10^–5. In this work we focus on the reflection spectra of exoplanet-exomoon systems instead of the thermal emission. By default the data useful for our data retrieval method requires order-of-magnitude reduction of stellar light, which is likely to be achieved through the combination of a space telescope and a coronagraph, the technology of which is still under development. Thus our method may not be useful in the next decade; it is not our intension to propose an exomoon detection method to compete with those already proposed. In this work we aim only at pointing out an interesting aspect of spectra analysis for data to be collected from exoplanet-exomoon systems in the future.

3.   Results and Discussions

3.1   Exoplanets/Exomoons Around Sun-like Stars

Figure 3 shows the photon spectra of 4 example target systems around Sun-like stars before photon noise is added, and the corresponding retrieval results with observation time of 1 hour. The properties of the targets and their corresponding retrieval are listed in Table 2. Observations of Figure 3 and Table 2 show that: 1) the retrieved photon spectra agree well with the target spectra (Figure 2); 2) the atmospheric pressure and the concentrations of water vapor and CH₄ in the retrieval results are identical to those in the targets; but 3) the other parameters (concentrations of CO₂ and O₂, A_i, B_j, r_planet, and r_ratio in the retrieval results deviate from those in the target atmospheres (Table 2).

Figure  3.  Photon spectra of 4 target systems and their retrieval results. Observation time is 1 hour in all cases. The properties of the target systems and retrieval results are listed in Table 2.

DownLoad: Full-Size Img PowerPoint

Table  2.  Properties of some targets and their corresponding retrieval results

Panels		Exoplanets												Exomoons
		Atmosphere Properties					Fractions of Surface Types						Radius	Fractions of Surface Types					rratio
		Pressure	f(H2O)	f(CH4)	f(CO2)	f(O2)	Augite	Basalt	Ocean	Olivine	Sand	Snow		B. P.	G. P.	Pla.	Oli.	Ilm.
a	Target	1 bar	9.25E-04	1.00E-04	1.58E-01	1.00E-01	0.15	0.09	0.28	0.11	0.12	0.26	1.00	0.09	034	0.18	0.15	0.23	0.26
	Retrieval	1 bar	9.25E-04	1.00E-04	2.00E-01	1.00E-01	0.07	0.00	0.66	0.04	0.11	0.13	1.37	0.58	0.00	0.00	0.42	0.00	0.20
b	Target	1 bar	1.58E-04	1.00E-05	1.00E-02	1.00E-04	0.25	0.10	0.22	0.29	0.01	0.13	1.00	0.25	0.21	0.31	0.08	0.16	0.23
	Retrieval	1 bar	1.58E-04	1.00E-05	1.00E-02	0.00E+00	0.19	0.05	0.46	0.21	0.00	0.09	1.18	0.48	0.22	0.14	0.16	0.00	0.20
c	Target	3 bars	1.58E-04	1.00E-03	1.00E-03	1.00E-02	0.17	0.26	0.19	0.18	0.15	0.05	1.00	0.34	0.03	0.14	0.21	0.28	0.39
	Retrieval	3 bars	1.58E-04	1.00E-03	1.00E-03	1.00E-02	0.10	0.00	0.28	0.14	0.48	0.00	0.93	0.32	0.00	0.21	0.24	0.23	0.40
d	Target	3 bars	7.75E-03	1.00E-03	1.58E-01	1.00E-04	0.21	0.13	0.27	0.25	0.08	0.07	1.00	N/A					0.00
	Retrieval	3 bars	7.75E-03	1.00E-03	2.00E-01	1.00E-02	0.16	0.43	0.28	0.11	0.00	0.02	1.20						0.00
Notes: B. P.=Brown Pyroxene, G. P.=Green Pyroxene, Pla.=Plagioclase; Oli.=Olivine; Ilm.=Ilmenite; f(gas) means the mixing ratio of corresponding gas; N/A means there is no satellite in the planetary system; 9.25E-04 means 9.25×10-4.

 | Show Table

DownLoad: CSV

The extremely good retrieval results observed in item 2 in the above is actually a consequence of including the target atmospheric composition in the database matrix M. A proof of this explanation is that when the target CO₂ or O₂ atmospheric concentration is included in database matrix M, the retrieval results of these species can also become extremely good. On one hand, this practice suggests that expanding the database will help improve the quality of retrieval. On the other hand, because we cannot foretell the composition of exoplanets and it is impractical to include infinitely large numbers of atmospheric compositions in our database, a statistical approach is necessary to evaluate the effectiveness of the retrieval method.

Figure 4 shows the statistical analysis of the retrieval results for CO₂ (blue in Figure 4a), O₂ (yellow in Figure 4a), exomoon (Figure 4b and c), and r_planet (Figure 4d). Δ is the difference between the retrieved values and those. For CO₂, it is clear that most of the retrieval are bound within ±0.5, an acceptable result for exoplanet observations. For exoplanet-exomoon target systems, the deviations of r_moon for >80% of the retrieval are smaller than 10% ( Figure 4b). The retrieval results for single exoplanet target systems have extremely low probability of false positive—very small percentage of them are retrieved as single planets (Figure 4c). Thus detection of exomoons similar to, or even smaller than, our Moon through analysis of reflection spectra of exoplanet-exomoon systems is promising. Most of the retrieved exoplanets have radii between 0.9 and 1.3 (Figure 4d) when the observation time is 1 hour. The effect of exoplanet size retrieval can be substantially improved with 1-day and 10-day observation times. Although the retrieval results for CO₂ are good, they are not good for O₂. In particular there is a high probability for ΔO₂=–1, which represents a strong trend for false negative: target systems containing O₂ in their atmospheres are interpreted as O₂-free. This situation is improved only slightly when the observation time is increased from 1 hour to 10 days.

Figure  4.  Retrieval results of hypothetical targets systems for 1-hour observations. Panel a is the deviation distributions of retrieved atmospheric O₂ (yellow) and CO₂ (blue), $\Delta \rm{Gas/Gas} = \frac{{\rm{Gas}\left({{\rm{retrieved}}} \right) - \rm{Gas}\left({{\rm{target}}} \right)}}{{\rm{Gas}\left({{\rm{target}}} \right)}}$ , which shows that our method is efficient for retrieving atmospheric CO₂. Panel b is the deviations of retrieved exomoon radii for planet-moon target systems, $\Delta {r_{{\rm{moon}}}} = {r_{{\rm{moon}}}}\left({{\rm{retrieved}}} \right) - {r_{{\rm{moon}}}}\left({{\rm{target}}} \right)$ , which shows that our method is effective at retrieving the size of existing exomoons. Panel c is the distribution of retrieved exomoon radii for single planet target systems, which shows that the probability of false positive is low. Panel d shows the distribution of retrieved planet radii (target radii are set to unity). The red lines in panels b & d represent the situation where the retrieved systems are identical to the targets: $\Delta {r_{{\rm{moon}}}} = 0$ and ${r_{{\rm{planet}}}} = 1$ .

DownLoad: Full-Size Img PowerPoint

Figure 5 shows the correlations between CO₂ and O₂ on targets (vertical axis) and in the retrieval results (horizontal axis). Large circles represent a higher percentage in all retrieval results. For demonstration purpose all retrieval results with less than 10^–4 retrieved O₂ are added to the count for O₂ =10^–4. It is clear that targets with CO₂ concentrations > 10 ^–4 can be retrieved rather accurately even with 1 hour observation time, consistent with the distribution of blue bars in Figure 4a. For targets with O₂ concentrations >10 ^–2, most of the retrieval results are good—close to the dashed line. If a retrieval shows that an exoplanetary atmosphere might contain 10^–2 O₂, the probability for the target atmosphere to have O₂ concentration less than 4×10^–3 is <10%. On the other hand, most of the retrieved atmospheres are with zero O ₂ if the atmospheric O₂ concentrations on targets are lower than a few times 10^–3—the wide spread of yellow bars inFigure 4a is mainly a result of hypothetically low target O₂ concentrations. We note that extending observation time from 1 hour to 10 days does not improve the situation of O₂ retrieval. Kaltenegger et al. (2007) suggests that atmospheres with O₂ concentration lower than 10^–3 cannot be observed. Because the best wavelengths to detect O₂ are 0.69 and 0.76 micron, which are outside of the wavelength range analyzed in this work, our conclusion regarding O₂ detection is more conservative than that in Kaltenegger et al. (2007).

Figure  5.  Distributions of CO₂ (left panel, 1 hour observation time) and O₂ (right panel, 10 day observation time). For demonstration purpose all retrieval results with O₂ less than 10^–4 are shown with O₂ = 10^–4. The dash lines represent exact retrieval of the gas concentration. The color of the circle symbols represent target gas concentration and the areas are proportional to the percentage of retrieval results.

DownLoad: Full-Size Img PowerPoint

The deviation distributions of surface type retrieval for 1 hour observation time are shown in Figure 6. It is clear that the retrieval effect of surface types are strongly affected by the characteristics of their reflectivity spectra (Figure 1). Ocean has the lowest reflectivity and the least spectral feature among all surface types considered in this work (Figure 1). Thus the retrieval result of oceans on exoplanets (Figure 5c) is worse than those of the other surface types. The situation of basalt is better than that of ocean (Figure 5b) but also has a broader distribution than the other surface types. When the observation time is increased from 1 hour to 1 day, the retrieval results of ocean and basalt become much better constrained (Figure 6a and b), which suggest that the data retrieval method can indeed be effective at retrieving bulk surface properties of exoplanets and exomoons.

Figure  6.  Deviation distributions of some retrieved surface types corresponding to Figure 4 (Sun-like stars + 1 hour observation). The definition of Δ is identical to that in Figure 4b. The deviation distributions of olivine on retrieved exoplanets are similar to those of snow (panel e). The deviation distributions of green pyroxene, plagioclase, olivine, and ilmenite on retrieved exomoons are all similar to that of brown pyroxene (panel f).

DownLoad: Full-Size Img PowerPoint

There is a weak negative correlation between the deviations of retrieved ocean fraction and basalt fraction. The reflectivity spectra of these 2 surface types are both flat, with that of basalt containing some features toward the long wavelength end. Thus there is a weak degeneracy between the fractions of these 2 surface types (Figure 6c). A strong positive correlation is present between the deviation of ocean fraction and the deviation of planet radii (Figure 6d)—the reflection spectrum of a retrieved planet with higher ocean fraction and smaller planet size is similar to that of the target spectrum because of the featureless reflectivity of ocean. Despite this, our method can still generate acceptable retrieval results for ocean fractions (Figure 7b)

Figure  7.  Deviation distributions of retrieved basalt and ocean fractions on hypothetical exoplanets around Sun-like stars with 1-day observation time (panels a and b), and the correlations between deviations of ocean fraction, basalt fraction, and planet radii (panels c and d). The retrieval of ocean fraction (panel b) is improved significantly in comparison to results from 1-hour observations (Figure 6c). The weak correlation between ΔBasalt and ΔOcean and strong correlation between $\Delta {r_{\rm planet}}$ and ΔOcean are the results of featureless reflectivity of basalt and ocean. Color from green to red represents retrieved r_planet varying from 0.8 to 1.5.

DownLoad: Full-Size Img PowerPoint

3.2   Exoplanets/Exomoons Around M dwarfs

The retrieval results for CO₂, O₂, exomoons, and r_planet (Figure 8 and 9) are similar to those for Sun-like stars. Because the photon flux used in this analysis for planets around M dwarfs is stronger than that for planets around Sun-like stars (red vs. blue in Figure 2), and the exoplanets used in the M dwarf analysis are larger than those used in the Sun-like star analysis, a careful comparison of Figure 9 and Figure 5 reveals that the distributions of retrieved O₂ or CO₂ in the M dwarf cases focus more closely to the target values (dashed lines) when observation times are identical. Nevertheless target atmospheres with O₂ concentration lower than a few times 10^–3 are still difficult to retrieve (Figure 9).

Figure  8.  Retrieval results using M dwarf photon spectrum (1 hour observation time). See Figure 4 for explanations of items in the panels. Note that r_planet in panel d is scaled to the planet radius used in producing the target photon spectra.

DownLoad: Full-Size Img PowerPoint

Figure  9.  Distributions of retrieval O₂ & CO ₂ using M dwarf photon spectrum. See Figure 5 for explanations of items in the panels.

DownLoad: Full-Size Img PowerPoint

The enhanced photon flux and larger planet size in the M dwarf analysis have more apparent effect on the retrieval of surface types. For Sun-like stars, the retrieved fractions of basalt and ocean have broad distributions when using the 1 hour observation time (Figure 6b and c). For M dwarfs, the distributions of retrieved basalt and ocean are narrower (Figure 10b and c), which is consistent with the improved performance when using longer observation times (Figure 7a, b vs. Figure 6b, c).

Figure  10.  Deviation distributions of retrieved surface type fractions corresponding to Figure 8. See Figure 6 for the meanings of all panels.

DownLoad: Full-Size Img PowerPoint

Figure  11.  Correlations between retrieval deviations of all surface types on exoplanets in this work. The change of color from green to red corresponds to an increase of the retrieved r_planet from 0.8 to 1.5.

DownLoad: Full-Size Img PowerPoint

Figure 7c demonstrates that there is a weak correlation between ΔBasalt and ΔOean because both surface types are featureless and have low reflectivity. Figure 11 shows the correlations in retrieval results between all surface types on exoplanets used in this work, with color from green to red represents retrieved r_planet varying from 0.8 to 1.5. Because the data retrieval method uses normalized r_planet on targets, red means that the retrieved exoplanets are larger than the targets and green means that the retrieved ones are smaller. The most apparent pattern in Figure 11 is that most of the red data points correspond to largest positive ΔOean (~0.6)—there is a strong degeneracy between retrieved planet size and retrieved ocean fraction. Because ocean has low reflectivity, a planet with a small ocean could be recognized as a larger planet with larger ocean fraction. What is interesting is that although basalt also has featureless and low reflectivity spectrum (higher than that of ocean), the red data points focus on negative ΔBasalt. One possibility is that although the overall reflectivity of basalt is low, it is actually greater than the reflectivity of other surface types at various wavelengths (Figure 1). Thus it is more difficult for the retrieval algorithm to misinterpret the fraction of basalt. The same trend (red data points focusing on negative Δ and green ones focusing on positive Δ) appears in the correlation patterns of all other surface types (Figure 11), which could indicate that an overcompensation between ocean fraction and r_planet is occurring in the data retrieval results.

The surface types with the tightest constraints are snow, augite, and olivine. The presence or absence of large fractions of snow and ice on an exoplanet will provide key information on planet climate. The variability of snow/ice could indicate seasonal cycles on exoplanets. Augite and olivine are two most common silicates on Earth, both of which are products of igneous processes. Olivine could be altered into iddingsite readily in the presence of water. Thus the presence or absence of olivine could provide indications of the efficiency of weathering processes on exoplanets. Thus the fact that the data retrieval method can provide reliable constraints on the fractions of snow, augite, and olivine is interesting. We emphasize that this work is only an exploratory test to establish the effectiveness of the data retrieval method. Much future work will be needed to perform careful and more proper selection of surface types.

4.   Conclusions

In this work a data retrieval method based on an analysis method for Martian surface minerals is developed to retrieve information based on near-IR (between 1 and 3 micron) reflection spectra of nearby rocky exoplanets. The analysis shows that the existence of exomoons, the size of exomoons, and the concentrations of some atmospheric species (H₂O, CH₄, and CO₂) on Earth-like exoplanets can be retrieved with high levels of reliability. In addition, the method can provide well-constrained fractions of some basic surface types (snow, augite, and olivine) on the targets because of the characteristic spectral features of atmospheric species and surface types in the analyzed spectral range.

We find that atmospheric O₂ concentrations lower than 10^–3 typically leads to false negative results (interpreted as zero O₂) even with 10-day observation time. This problem is caused by the weakness of O₂ absorption feature, which could be easily mixed up with other species. This conclusion is consistent with the findings in Kaltenegger et al. (2007).

Our analysis is based on observations in 40 bands uniformly spaced in wavenumber (spectral resolution R=20–60). Because data from different bands are not equally valuable for characterization of exoplanets, observation capability tailored specifically for exoplanet characterization will be helpful to further improve retrieval efficiency. Our analysis shows that because planets at the Earth-equivalent distance from an M dwarf receive more NIR radiation than the Earth receives from the Sun, less observation time will be needed to achieve similar retrieval effects for exoplanets orbiting M dwarfs. We also conclude that including more samples of transmission spectra of atmospheric species in the database will improve the quality of retrieval.

One underlying assumption in this work is that the same combination of surface types is observed repeatedly. Although exoplanets are orbiting around their stars and rotating, knowledge on the rotation rates and phases of habitable zone exoplanets and exomoons around M dwarfs can be obtained prior to observation planning and data analysis. Nevertheless multi-orbit observations using future large space telescopes will be required to accumulate adequate photons for data retrieval described in this paper. The orbital periods and phases of exomoons will likely remain unknown to us and our analysis did not consider this issue explicitly. The multi-orbit observations described previously could contain reflections from exomoons during different phases, which could add one more dimension to future data analysis. If observations are carried out with high enough frequency, it might be possible to obtain periodic spectra variations due to phase changes of exomoons, which will be an interesting and independent method to confirm the existence of detected exomoons through photometry.

The main uncertainty that might influence the results of our retrieval method is the cloud effect, which will be important to address in future work. In addition, the non-linear effects of grain size on data retrieval (Schmidt et al., 2014) and instrument noise should be explored. Despite these caveats, our analysis shows that in principle atmospheric composition, distributions of surface types, and the presence of exomoons could be retrieved simultaneously by analyzing NIR reflection spectra (with moderate spectral resolution R=20–60) of nearby Earth-size rocky exoplanets.

Acknowledgement

Wenliang Cui, Duo Cui, Tong Li, Feng Tian are supported by the National Natural Science Foundation of China (1166116101, 4164104). , Xiaomeng Huang is supported by the National Natural Science Foundation of China (41375102). Jinsu Zhang acknowledges supports from the College Student Research and Career-Creation Program of Beijing City (201710003B057) and the Spark Program at Tsinghua University. Frederic Schmidt acknowledge support from “Institute National des Sciences de l’Univers” (INSU), the "Centre National de la Recherche Scientifique" (CNRS) and "Centre National d’Etude Spatiale" (CNES), through the "Programme National de Planétologie", the PICS program and MEX/OMEGA and MEX/PFS programs. We thank Hao Zhang at the Chinese University of Geosciences for help on confirming sources of mineral spectra data. This research utilizes mineral spectra acquired by Crustal Imaging and Characterization Team and Central Region Minerals Team from USGS (led by R. N. Clark), and C. Pieters’s group with the NASA RELAB facility at Brown University.