© ESO, 2011

1. Introduction

The atmospheric properties of transiting planets can be studied through eclipse spectroscopy. The molecular composition can be constrained by “transmission” spectroscopy during the primary transit or by emission spectroscopy at the secondary eclipse. This has been achieved from both space and the ground for hot Jupiters (e.g. Tinetti et al. 2007, transmission, Spitzer), hot Neptunes (e.g. Stevenson et al. 2010, emission, Spitzer), and more recently for GJ1214b, a 7 M_⊕ planet (e.g. Bean et al. 2010, transmission, VLT). Transits also provide constraints on the thermal structure of the atmosphere (e.g. Knutson et al. 2008, emission, Spitzer), as well as on circulation (e.g. Snellen et al. 2010, transmission,VLT). The thermal phase curves of transiting planets (variation in the apparent infrared emission of the planet with the orbital phase) have been observed in a single band (e.g. Knutson et al. 2009a, Spitzer) or at different wavelengths (e.g. Knutson et al. 2009b, Spitzer), providing constraints on the day/night brightness temperatures and dynamics.

Fig. 1 The orbit and phases as seen by the observer. The inclination is set to 60°. The phases shown here are numbered from 1 to 12 and correspond to the phases presented in Figs. 3 to 5.

However, applying these methods to characterize (or simply detect) the atmosphere of terrestrial exoplanets with the forthcoming generation of telescopes requires the discovery of transiting systems very close to the Sun (roughly within 10 pc). Short-period planets (<50 days) in the 5−20 Earth mass range have been found by the HARPS radial velocity survey around 30 ± 10% of the G and K stars (Mayor et al. 2009b). For the time being, the closest transiting planet with a mass lower than 20 M_⊕ is GJ1214b, found by the MEarth survey at a distance of 13 pc (Charbonneau et al. 2009). Statistically, there should be one transiting planet from this population for distances between 8 and 18 pc (depending upon orbital periods), and possibly as close as 5 pc if we assume a similar population around M dwarfs. Less than a handful of terrestrial targets will thus be suitable for atmosphere characterization through transit spectroscopy with JWST (Belu et al. 2011). Within this population of short-period, low-mass planets, there are  ~10 times more nontransiting planets then transiting ones. A technique able to characterize nontransiting planets and their atmospheres would thus provide enough targets to explore the diversity of terrestrial exoplanet atmospheres.

Ambitious projects aimed at directly detecting of terrestrial planets, such as Darwin (Cockell et al. 2009), TPF-I (Lawson et al. 2007), TPF-C (Levine et al. 2009), and New Worlds (Cash et al. 2009), could measure the planet spectrum (in the visible, near infrared, or mid-infrared, depending on the technique), and its variability. Variations in the broadband mid-IR flux with the orbital phase, have been suggested in the framework of Darwin/TPF-I as a way to identify planets with dense atmospheres (Selsis 2004), a necessary criterion for surface habitability, or to constrain obliquity (Gaidos & Williams 2004). Variations of the reflected light (Pallé et al. 2008) or the infrared emission (Gómez-Leal et al., in prep.) could reveal cloud patterns and the rotation rate of the planet. Until these large observatories have been developed, we will have to study terrestrial exoplanets by observing the combined star+planet light. Although the absolute planetary flux cannot be extracted from such observations, the modulation of the flux at the orbital period of the planet can be attributed to the planet (assuming no pollution of the signal by intrinsic stellar variability or systematic instrumental effects). The phase curve of a hot Jupiter has already been detected in infrared bands in a nontransiting configuration (Cowan et al. 2007; Crossfield et al. 2010) and not detecting the phase curve of the  ~7 M_⊕ planet GJ876d has been used to discuss the presence of a dense atmosphere attenuating the day/night difference in brightness temperature (Seager & Deming 2009). The phase- and wavelength-dependent emission of hot Jupiters have been simulated using either 1D models assuming a local radiative-convective and chemical equilibrium (Barman et al. 2005) or circulation models using a Newtonian scheme for the temperature calculation and either equilibrium chemistry (Burrows et al. 2010) or kinetic calculations (Fortney et al. 2006). The phase-dependent visible flux and its dependence on inclination has been studied for nontransiting by (Kane & Gelino 2011). Cowan & Agol (2011) studied the thermal phase curve of hot Jupiters, including eccentric ones, with a semi-analytic model.

In the present article we address the characterization of the atmosphere of a hot terrestrial planet that does not transit its host star, using spatially-unresolved spectro-photometry of the system. We study how the apparent emission of the planet varies with wavelength and orbital phase, for different atmospheric pressures, in the case of a tidally-locked planet. In another study (Maurin et al. 2011), we discuss the retrieval of the radius, albedo, and inclination of airless planets. In this article we focus on the influence of the atmosphere, as simulated with a 3D climate model.

2. Model

Modeling the orbital phase variations of the thermal emission of a planet observed by a distant observer requires us to model the 3D structure of the atmosphere, including the temperature, pressure, chemical composition, and cloud / aerosol content. Due to the computation time of 3D atmospheric simulations, the possible compositional diversity of planetary atmospheres, and the variety of planetary systems, this preliminary study does not aim to explore the effect of all the parameters controlling the observables (star type, orbital elements, planet rotation, planet size and gravity, atmospheric elemental composition, nature of the surface). At this stage, we chose to consider a simple planet corresponding to a system that we can expect to detect (for instance by radial velocity) in the vicinity of the Sun: a short-period massive terrestrial planet (sometimes called a hot super-Earth) around an M dwarf. The planet we model may be representative of known nearby low-mass exoplanets, which we list in Table 1. Our example planet is on a circular, synchronized orbit, with a null obliquity, a state consistent with strong tidal interactions (see Sect. 4.2). We consider a fixed atmospheric composition and only vary the atmospheric pressure, in order to study how the atmosphere modifies the phase- and wavelength-dependent properties of the apparent infrared emission. The characteristics of the planetary system considered are summarized here:

• The star is an M 3 dwarf (0.31 M_⊙) with a luminosity of 0.0135 L_⊙ and the same spectral distribution as AD Leo (we used the Virtual Planet Laboratory AD Leo data from Segura et al. 2005). For flux and photon-noise calculations, the distance of the star is set to 10 pc.

• The orbit is circular with a semi-major axis of 0.0535 AU and period of 8 days. As seen in Fig. 1, the inclination of the orbit on the line of sight is set to 60°, which the median value for randomly oriented systems.

• The planet is a 1.8 R_⊕ and 9.5 M_⊕ terrestrial planet, consistent with an Earth-like bulk composition (Valencia et al. 2006; Sotin et al. 2007). Its surface gravity is 3 g. We set the surface albedo to 0.2. The planet rotation period is assumed to be equal to its orbital period (8 days).

• The atmosphere is made of carbon dioxide (CO₂) as the unique constituent. We consider three atmospheric pressures: 0.1, 1, and 10 bars, as well as a case with no atmosphere.

In the case with no atmosphere, the surface temperature is in local equilibrium with the stellar irradiation and given by ${\mathit{T}}_{\mathrm{eq}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\mathrm{\left[}\mathit{S}\mathrm{\right(}\mathrm{1}\mathrm{-}\mathit{A}\mathrm{\left)}\cos \mathrm{\right(}\mathit{\theta }\mathrm{)}\mathit{/}\mathit{\sigma }{\mathrm{]}}^{\frac{\mathrm{1}}{\mathrm{4}}}$(1)where S is the stellar flux at 0.05 AU, A the surface albedo, σ the Stefan-Boltzmann constant and θ the zenith angle. The emissivity of the surface is set to 1. On the dark hemisphere, the temperature is set to 37 K. This value is obtained by scaling the intrinsic luminosity of the Earth (~30 TW, due mainly to radioactive decay) to the mass and radius of the planet. Using a night-side temperature of 0 K would not, however, change the disk-integrated properties in an observable way. The substellar temperature is 553 K, the equilibrium temperature is 390 K (stellar irradiation uniformly distributed over the whole planet), and the equilibrium temperature calculated only for the starlit hemisphere is 465 K (stellar irradiation uniformly distributed over the starlit hemisphere only). We do not consider unsynchronized or eccentric cases, and so the irradiation at a given location on the planet always remains the same. Therefore we do not need to include the effects of surface thermal inertia.

For the three cases with an atmosphere, we used the LMD 3D atmospheric general circulation model (GCM). This model, originally developed for the Earth (Sadourny 1975; Hourdin et al. 2006), has been adapted to most terrestrial planet atmospheres in the solar system like Mars (Forget et al. 1999), Titan (Hourdin et al. 1995), Venus (Lebonnois et al. 2010) and even Triton (Vangvichith et al. 2010). Because GCMs are based on physical equations that apply across wide ranges of pressure, temperature, and atmospheric composition, they can simulate accurately the thermal structure and winds of a broad diversity of atmospheres, provided that basic data such as spectroscopic properties are available in the regime of interest. In this study, we used a new, generic version of the LMD GCM that has been developed specifically for exoplanet and paleoclimate studies. It has already been used to study the early Mars climate (Forget et al. 2010) and possible climates on GJ 581d (Wordsworth et al. 2011). In practice, this version uses the LMDZ 3D dynamical core (Hourdin et al. 2006). It is based on a finite-difference formulation of the classical primitive equations of meteorology. The primitive equations are a simplified version of the Navier-Stokes general equations of hydrodynamics based on three main approximations: the atmosphere is assumed 1) to be a perfect gas, 2) to remain in hydrostatic equilibrium vertically, and 3) to have much smaller vertical dimensions than the radius of the planet (thin-layer approximation). In the LMDZ GCM, the potential enstrophy and total angular momentum are numerically conserved for barotropic flows (Sadourny 1975). Scale-selective hyperdiffusion was used in the horizontal plane for stability, and linear damping was applied in the topmost levels to eliminate the artificial reflection of vertically propagating waves. The planetary boundary layer is parameterized using implicit timestepping and the method of Mellor & Yamada (Mellor & Yamada 1982; Galperin et al. 1988) is used to calculate turbulent mixing. In addition to this parametrization, a convective adjustment scheme is used to prevent subadiabatic vertical temperatures gradients. The temperature of the surface is computed from the radiative, sensible and latent heat fluxes at the surface using an 18-level model of thermal diffusion in the soil and assuming a homogeneous thermal inertia of 250 J m^-2 s^−1/2 K^-1 and a surface albedo of 0.2.

GCMs have already been applied to terrestrial exoplanets, especially to study their habitability when locked in synchronous rotation (Joshi et al. 1997; Merlis & Schneider 2010; Pierrehumbert 2011; Heng et al. 2011). The main improvement of the model we use is the radiative transfer scheme, which is similar to the one described in Wordsworth et al. (2010b). We computed high-resolution spectra over a range of temperatures and CO₂ pressures using a 14  ×  9 temperature-pressure grid with values T = 100−750 K, p = 10^-3−10⁵ mbar. The correlated-k method is then employed to produce a smaller database of coefficients suitable for fast calculation in a GCM. The model uses 36 bands between 0.3 and 5 μm for the incoming stellar radiation and 38 bands above 2.5 μm for the planetary infrared emission. The overlap is due to the fact that no source function is considered for the stellar light transfer, only absorption and scattering. Sixteen points were used for the g-space integration, where g is the cumulated distribution function of the absorption data for each band. CO₂ collision-induced absorption (CIA) was included using a parameterization based on the most recent theoretical and experimental studies (Wordsworth et al. 2010a; Gruszka & Borysow 1997; Baranov et al. 2004). A two-stream scheme (Toon et al. 1989) was used to account for the radiative effects of Rayleigh scattering, which was included by the method described in Hansen & Travis (1974). Condensation and CO₂ clouds were not included in these simulations, which is consistent with the range of temperature and pressure found in the 1 and 10 bars cases but not in the 0.1 bar case (see Sect. 4.1).

Fig. 2 Equatorial temperature (top) and horizontal wind maps (bottom) at the surface (left) and 40 km (right). Temperatures are given for 0.1, 1 and 10 bar; wind maps are shown for the 10 bar case only. The displacement of the hot spot, compared to the substellar point (longitude = 0°, latitude = 0°) can be seen on the surface plots.

The modeled atmospheres were made of a single atmospheric constituent (CO₂). The GCM could include more species and real atmospheres are not thought to consist of only one constituent, but this has several advantages for a preliminary study. A key point of this study is to show that the variation spectrum (peak amplitude of the phase curves as a function of wavelength), which can be measured in relative photometry for a nontransiting planet, is strongly correlated with the emission spectrum and exhibits molecular signatures. We therefore preferred to study a single component atmosphere as a first step, so that the complex relation between the composition, the spectral- and spatial distribution of the emission, and the amplitude spectrum of the phase curves could be understood more easily. This pure CO₂ case can be seen as an illustration of the principles of this characterization method. Future, more realistic studies will include atmospheric mixtures and, in particular, water vapor and clouds.

For this study the GCM resolution is 32 × 24 × 16 grid (longitude  ×  latitude  ×  altitude). This is found to be high enough to capture the most important features of the atmospheric dynamics, given that it resolves the Rossby deformation radius ${\mathit{L}}_{\mathrm{R}}\mathrm{~}\sqrt{\mathit{RT}}\mathit{/}\mathrm{2}\mathrm{\Omega }\mathrm{~}\mathrm{2.5} {\mathit{R}}_{\mathrm{\oplus }}$, where R is the specific gas constant, Ω is the planetary rotation rate, and T is the atmospheric temperature. We did not test the effect of our model’s finite difference numerical scheme on the results. However, this issue has been studied by Heng et al. (2011), who found quantitative agreement between spectral and finite difference dynamical cores for simulations of a hypothetical tidally locked Earth.

The model starts with an isothermal surface and atmosphere at 300 K. It computes the evolution of the atmosphere during 100 orbits, after which the state of the atmosphere is close enough to a steady state. Using the last orbit, we produce phase curves for 30 of the 38 infrared bands (given in Table 2), between 2.5 and 20 μm. To do that, we use the top-of-the-atmosphere longitude-latitude maps of outgoing fluxes computed by the GCM with a timestep of 1 h (192 points on the orbit). At a given time and wavelength, we sum all these outgoing fluxes, weighted by the area of their corresponding cell and by the cosine of the angle between the normal to the surface and the direction towards the observer. We assume an isotropic distribution of specific intensities at the top of the atmosphere (see Sect. 4.9 for a discussion on the validity of this approximation). For the case with no-atmosphere, the surface emission is calculated using the same longitude-latitude cells and the blackbody emission is integrated over the same spectral bands. The procedure for deriving the phase curve is then exactly the same.

3. Results

The results of our model are presented for initial pressures of 0, 0.1, 1 and 10 bar. In Fig. 2, we present some characteristics of the atmospheric circulation and associated heat transport we have obtained with the model. Wind maps (given in the 10 bar case) show that surface winds converge toward the substellar area, feeding the convective plume generated by the stellar heating. Coriolis forces due to the planet’s rotation produce an east-west asymmetry in the wind field. In the upper atmosphere, super-rotation occurs with an equatorial jet slightly faster than the planet rotation itself. In the uppermost layer, wind speeds become close to the speed of sound, and we may be approaching the limit of our model validity. However, this has negligible consequences on the dynamics at higher pressures that produce the thermal emission. The super-rotation regime starts for pressures above about 1 bar and is not found in the 0.1 bar case.

In the surface wind map for the 10 bar case, the convective plume does not start at the exact substellar location (longitude = 0°, latitude = 0°) but is slightly shifted eastward. This shift of the hot spot is due to heat transport by the atmosphere. It has been predicted for hot Jupiters (Showman & Guillot 2002) and was observed with Spitzer (Knutson et al. 2007) as a lag between the maximum of the phase curve and the secondary eclipse. This delay is also found in our model depending on the wavelength and the probed altitude. Figure 2 also shows the equatorial temperature as a function of longitude. On the surface, a 5 and a 10° displacement of the hot spot can be seen for the 1 and 10 bar atmospheres, respectively, corresponding to a 2.5 and a 5 h lag. However, at 40 km and for the 1 and 10 bar cases, the longitudinal thermal profile is much more complex and no main maximum can be identified. The existence of a clear maximum and its associated lag thus depends strongly upon the wavelength. If observed in a single band, a lag can be regarded as evidence of an atmosphere only if the planet is definitely synchronized. Indeed, surface thermal inertia can produce a similar hot spot displacement if the planet is not synchronized. In particular, the case of eccentric planets that can be either in a pseudo-rotation state or in various spin-orbit resonances (Leconte et al. 2010) must be treated with caution.

Fig. 3 Planetary emission maps as a function of wavelength and phase. Phase numbers refer to Fig. 1.

Fig. 4 Disk-integrated flux as a function of wavelength and phase.

Fig. 5 Brightness temperature maps as a function of wavelength and phase (top) and gas temperature maps as a function of altitude and phase (bottom).

How efficiently the heat deposited by insolation is redistributed to the night side and to the poles can be evaluated by comparing analytical radiative and dynamical timescales. The mean radiative timescale ${\mathit{\tau }}_{\mathit{R}}\mathrm{=}\frac{{\mathit{c}}_{\mathit{P}}\mathit{P}}{\mathit{\sigma g}{\mathit{T}}^{\mathrm{3}}}$ where c_P is the heat capacity of CO₂ (assumed constant here following Wordsworth et al. 2010b) and T is the mean bolometric brightness temperature of the planet (~400 K in all cases). The mean dynamical timescale where R is the planetary radius and is the density-weighted average of the wind speed. The ratio τ_R/τ_D is  ~0.1 for P = 0.1 bar,  ~1 for P = 1 bar, and  >5 for P = 10 bar. Heat redistribution is thus efficient for pressures above 1 bar.

Fig. 6 Planetary signal as a function of time. The 8-day orbits is divided in 7 exposures. The number of photons is given per exposure, per band, for a 6 m telescope and a 10 pc target. The dotted line represents the amplitude of the stellar photon noise. Wavelengths are indicated in the upper left corner of each plot.

Fig. 7 Same as Fig. 6 but for a 1.5 m telescope.

Fig. 8 Same as Figs. 6 and 7 but the photon flux is given here as the planet/star contrast ratio.

Figures 3 and 4 show the top-of-the-atmosphere outgoing fluxes, as seen by the observer, as a function of wavelength (Y-axis) and phase (X-Axis). The phases represented and their numbering are the same as in Fig. 1. Figure 3 gives the spatial distribution of thermal emission (in arbitrary units proportional to W m^-2μ m^-1). The aim of this figure is to show the spatial distribution of the emission on the apparent planetary disk. Figure 4 gives the corresponding disk-integrated fluxes calculated at 10 pc. The bond albedo of the planet is only slightly affected by the atmospheric pressure: its value is 0.20, 0.22, 0.21, and 0.205 for 0, 0.1, 1, and 10 bars, respectively. Therefore, the global thermal emission of the planet in Watts does not change between the different cases, but its spectral and spatial distribution strongly depends on the pressure due to the effects of the atmosphere on the radiative and heat transport. With increasing pressure, the planet tends to radiate predominantly in the spectral atmospheric windows (where CO₂ absorbs less), while the emission in the dark hemisphere becomes more important due to the more efficient heat transport through atmospheric circulation.

At 0.1 bar, the absorption features of CO₂, in particular the 2.7, 4.3, and 15 μm bands, already affect the spectral distribution, but the thermal emission on the night side remains negligible. In atmospheric windows, for instance between 5 and 8 μm, the flux variations do not differ much from the airless case, but there is slightly higher maximum emission, as surface temperatures are increased by greenhouse warming. At higher pressures, weaker CO₂ bands (for instance between 9 and 11 μm) and collision-induced absorption (dimer absorption at 7−8 μm) also shape the emission spectrum. Night-side cooling is noticeable at 1 bar and becomes a significant contribution to the total emission at 10 bars. Figure 5 presents flux maps similar to those in Fig. 3 but converted into brightness temperatures T_B (=temperature of a blackbody with the same radiance at the considered wavelength). At the bottom of Fig. 5, we have mapped the gas temperature at different altitudes on the observed planetary disk. With increasing pressure, the surface and the lower atmosphere become hotter owing to greenhouse warming and more uniform thanks to efficient heat transport (see Table 3) but the atmosphere also becomes more opaque at a given wavelength. Therefore, the observed emission comes from higher altitudes. With the exception of some nightside regions of the 0.1 case, and some high-altitude substellar inversion in the 1 and 10 bar cases, temperature monotically decreases with increasing altitude. This allows us to roughly determine the altitude at which the photons of a given wavelength are emitted (the τ ~ 1 transition region between the optically thick and thin atmospheric layers) without computing the contribution function, by comparing the maps of T_B(λ) and T(z). For instance, we can see that the 4.26 μm band (framed in Fig. 5) is emitted around 10 km, 35 km, and 60 km, in the 0.1, 1, and 10 bars simulations, respectively . Therefore, the increase in temperature in the lower atmosphere owing to enhanced CO₂ level results in higher fluxes only in atmospheric windows. Windows open at low pressure are closed at high pressures by the CIA continuum and weak features. Also, the day-night temperature contrast is more pronounced at the surface and in the lower atmosphere than at higher altitudes. As a consequence, wavelengths probing the lower atmosphere exhibit large phase-dependent variations, while wavelengths falling in opaque regions of the planet spectrum and probing higher levels exhibit small variations. This can be seen for instance at 8.7 μm, which probes the first few kms of the atmosphere in the 0.1 and 1 bar cases, exhibiting high fluxes and large variations. At 10 bar, the lower atmosphere becomes optically thick due to CO₂-CO₂ dimer absorption and the flux emerges above 10 km with smaller values and variations.

3.1. Observability of the phase curves

The variation in the apparent emission of the planet produces a modulation of the observed star+planet combined flux. Assuming that the planet is detected by radial velocity, its orbital period is known. This allows the observer to extract the flux variations at this specific period. We do not address here the issues of stellar variability and instrumental stability that certainly represent the main obstacles in extracting the planet signature. As a first step, we simply compare the planet emission with the stellar-photon noise to check that the observation can at least be considered. To calculate the stellar photon noise, we need to assume a collecting area, an exposure time, a spectral resolution, and a target distance. In this study, the spectral resolution is imposed by the bandwidths of our radiative transfer scheme (band limits are given in Table 2). The corresponding resolution varies from  ~25 at 2.5 μm to  ~10 at 20 μm. The target distance is set to 10 pc. We assume that one complete orbit of the planet is observed with an exposure time set to one seventh of the orbital period (27.4 h). An odd number of exposures allows us to center one exposure on the maximum phase.

Fig. 9 Planetary spectrum vs. variation spectrum. Circles, squares, and diamonds correspond the planetary spectrum at 3 different phases indicated at the top of the graph (phases 6, 3, and 12 in Fig. 1). The dashed line shows the amplitude of the 7-point phase curves as a function of wavelength (i.e. variation spectrum).

Figures 6 and 7 show the multiband phase curves for telescope diameters of 6 and 1.5 m, respectively. These values correspond to the size of JWST¹ and EChO, a space telescope proposed to ESA for Cosmic Vision 2. With a 6 m ideal collector, the signal to photon-noise ratio exceeds 5 for λ > 3.5   μm and 10 for λ > 3.5   μm, in the scenario giving the largest variations, which is the case with no atmosphere. With 1.5 m, SNRs of 5 and 10 are obtained at wavelengths above 4.5 and 8.5 μm, respectively. Figure 8 gives the same phase curves but the flux is given relative to the star. We can see that a photometric precision on the order of 10^-5 is required to resolve the planetary signal, which is the aimed precision for EChO in the infrared and the precision reached by Kepler in the visible.

The phase curves have been computed using the thermal emission given by the GCM and neglecting the stellar light scattered by the planet. As shown in Fig. 10, the contribution of the reflected light is only significant at wavelengths below 3.5 μm. At these short wavelengths, the planet/star ratio is lower than 10^-6, making the observation of the planetary signal out of reach even for the next generation of infrared detectors. At higher wavelengths, where the extraction of the planetary emission is conceivable, we can thus safely neglect the reflected light.

3.2. The variation spectrum

Observation of the unresolved star+planet light does not give access to the absolute planetary emission (unless a secondary eclipse is observed). What can be measured is the variation in the planet emission. Valuable information about the nature of the atmosphere can be obtained from the variation spectrum: the peak amplitude of the phase-dependent variations as a function of the wavelength. Figure 9 shows the variation spectrum compared to the planetary spectrum at three different phases. For the 0 and 0.1 bar cases, the mid-IR variation spectrum is extremely close to the emission spectrum at maximum phase. This is because the night side emission is extremely weak. In the 1 bar case, the variation spectrum and the maximum phase spectrum are similar between 3 and 9 μm. Above that, there is no atmospheric window probing the lowermost part of the atmosphere and the variations become very small. At 10 bars, 2 windows remain, at 3.5 and 6 μm, plus a marginal one at 9 μm. At other wavelengths, the planet spectrum is basically phase-independent. In all cases with an atmosphere, the variation spectrum is shaped by the radiative properties of CO₂, and the presence of CO₂ (or another IR absorber) can be inferred from the variation spectrum.

4. Discussions

4.1. CO₂ condensation

Our simulations do not include CO₂ condensation. In the runs at 1 and 10 bars, the pressure found in a steady state never exceeds the vapor pressure of CO₂. In the 0.1 bar case, however, significant supersaturation is found at the surface and in the first 200 m above in some regions of the night side. Figure 11 shows the areas and altitudes where the pressure exceeds the CO₂ vapor pressure. Although the formation of surface CO₂ ice and CO₂-ice clouds in the regions shown in Fig. 11 would have a negligible effect on the disk-integrated thermal emission, these regions could eventually trap the whole atmosphere into a surface deposit of CO₂ ice. A steady state would imply a release of CO₂ balancing the CO₂ loss into this cold trap, which may not be realistic in terms of timescales and which represents an ad hoc situation. Atmospheres modeled without CO₂ condensation require a minimum CO₂ pressure of about 0.5 bars to remain fully unsaturated. The 0.1 bar case presented in this study is thus not self-consistent. Simulations for 0.5 bars produce phase curves very similar to those obtained with 1 bar. It is possible that a pure CO₂ atmosphere could only exist for pressures above 0.5 bars (for the stellar, orbital and planetary parameters considered here). This would mean that once the planet is synchronized, its atmosphere must have already been outgassed, otherwise it would not be able to build up because the released gas would end up as CO₂ ice on the night side. Analysis of this interesting issue is beyond the scope of this study. Nonetheless, we decided to keep the 0.1 bar case, for two reasons:

• Our paper focuses on the information contained in the spectral phase curves that can be measured for nontransiting planets, and comparison of the four chosen cases (0, 0.1, 1, and 10 bars) illustrates the effects of increasing pressure on the spectral and spatial distribution of the thermal emission.

• When assuming a pure CO₂ atmosphere and neglecting CO₂ condensation, we find that the mean surface pressure must be higher than  ~0.5 bar to avoid supersaturation in the dark hemisphere. However, this minimum pressure is expected to be weaker if condensation is actually included with all its associated effects, in particular the release of latent heat and the radiative effect of CO₂-clouds. As the CO₂-ice particles form on the night side and blanket the surface thermal emission, both these effects would tend to warm the region where condensation occurs. Therefore, 0.5 bar is probably an overestimate of the minimum stable CO₂ pressure.

This link between CO₂ condensation and the stability of the atmosphere against the night-side cold trap is an interesting subject for future study.

Fig. 10 Contribution of the reflected light to the phase curve with no atmosphere. The 7-point phase curves are computed with emission only (black circles) and with emission+reflection (diamonds). Band central wavelengths are indicated in the upper left corner of each plot.

Fig. 11 Condensation of CO2 in the 0.1 bar case. Color contours indicate the surface temperature, with the same color scale as in previous temperature maps. Contour lines indicate the region where CO2 is supersaturated at the surface, at 4 m, 10 m and 230 m. Levels above 230 m are not saturated. The 0°−0° lon-lat coordinate corresponds to the substellar point.

4.2. Nonsynchronized planets

We restricted our study to a tidally-locked planet on a circular orbit. Atmospheric thermal tides can result in an equilibrium rotation period that is different from the synchronized one (Correia et al. 2008). This is for instance at the origin of the current rotation state of Venus (Correia & Laskar 2001). Therefore, assuming a synchronized planet may not be consistent with a dense atmosphere. Correia et al. (2008) show that the ratio ω/n, where ω is the rotation rate and n the mean motion, scales as ${{\mathrm{(}}^{}\mathit{a}{\mathit{M}}_{\mathrm{\ast }}{\mathrm{)}}^{}}^{\mathrm{2.5}}$ for a given planet and atmosphere. This ratio is 1.92 for Venus. A Venus-like planet at 0.05 AU from a 0.31 M_⊙ star would thus have a ratio of as small as 5 × 10^-5 (which means  ~440 days between sunset and sunrise). Our simulated planet is 11 times more massive than Venus, while its 10 bar atmosphere is  ~30 times less massive than that of Venus, which would make the ratio ω/n even smaller. The departure from synchronization is therefore expected to be far too small to affect our results.

4.3. Photon-noise limited observations

In this paper, the only source of noise we consider is the stellar photon shot-noise, which is clearly a preliminary first step towards a realistic assessment of the observability of phase curves. The astrophysical background and zodiacal light should also be included in a realistic noise estimation. For JWST, one should take available estimates of the known sources of noise (thermal, read-out, dark current) into account for the instruments suitable for spectro-photometric observations, in particular MIRI and NIRSpec. Including these sources of noise should typically yield a factor 2 increase in the noise level (Belu et al. 2010). The stability of the instruments over one or more orbital periods is a key point that should be addressed in depth considering the required photometric precision (10^-4 to 10^-5). EChO is a 1.5 m telescope proposed to the program Cosmic Vision 2 of ESA and dedicated to the characterization of exoplanet by combined star+planet light spectro-photometry. The aim for EChO is to achieve photon-noise limited observations and high detector stability over long periods. Our signal-to-noise basic estimates are therefore more relevant to the case of EChO.

4.4. Stellar variability

Stellar variability is key issue that cannot be solved directly by improving the instrumentation. The intrinsic stellar photometric variations and their wavelength-dependence, at a level of  ~10^-5, represent an obvious impediment to the observation of nontransiting terrestrial planet phase curves. Basri et al. (2010) studied the photometric variability of  ~150 000 Kepler stars over 33.5 days with a cadence of 30 min. This work shows that fewer than about 20% of main-sequence stars exhibit a variability lower than  ~10^-3, and this fraction tends to decrease for low-mass dwarfs. What can contaminate the planetary phase curve, though, is the stellar component that varies periodically with the known period of the planetary modulation. This means that a minimum of two orbital periods will have to be observed in most cases in order to filter out anything that is aperiodic or periodic at another period. In-depth characterization of the stellar variability, including its origin, periodicity, and wavelength-dependence, is certainly required prior to any attempt to measure the phase curves of terrestrial planets. On an optimistic note, the optical phase variations of Kepler 10 b with a peak amplitude of only a few ppm has been measured, showing that an accuracy of 10 ppm must be achievable, at least in some cases. Our expertise is primarily in planetary physics, and we encourage experts of stellar variability to use their models and observations to assess the actual observability of the planet signature.

4.5. Orbital distance and stellar types

An 8-day period planet around an M 3 dwarf (the case modeled in the present study) is not the most favorable case for attempting a phase curve observation of a terrestrial exoplanet. Hotter planets, with a shorter orbital period, are frequent and can provide better targets, depending on the stellar type and the wavelength. Known terrestrial exoplanets can be as hot as a few thousand K, such as Corot-7b (Léger et al. 2009) and Kepler-10b (Batalha et al. 2011), and Table 1 shows that planets with terrestrial masses that are hotter than our model are known in the immediate neighborhood of the Sun.

Fig. 12 Amplitude of the phase curve for an airless planet as a function of stellar mass, orbital period and wavelength. The albedo of the planet is 0.2 and the inclination is 60°. Grey circles indicate the case modeled in this paper. Stellar radii and Teff are from Baraffe et al. (1998).

Figure 12 shows how the peak amplitude of the thermal phase curve (given in planet-to-star flux ratio) varies with the stellar type for two fixed orbital periods (2 and 8 days) and three wavelengths (3, 9, and 15 μm). One can see that phase curve at wavelengths around 3 μm can only be measured for very short-period objects around Sun-like stars. The search for the 2.7 and 4.3 μ CO₂ signatures may be restricted to such systems. We chose to model a planet below 600 K because one can argue that dense atmospheres may be rare around hotter planets due to atmospheric erosion. However, the existence and nature of dense atmospheres on terrestrial exoplanets, and their frequency as a function of orbital distance, stellar type, planet mass, must be eventually addressed by observation. It is therefore important to also target very hot planets, keeping in mind that the multiband phase curves are easier to measure in the absence of an atmosphere due to the high amplitude of the phase variations. For airless planets, multiband phase curves can be used to constrain the orbital inclination, the planetary radius and the surface albedo (Maurin et al. 2011).

4.6. Realistic atmosphere composition and their influence on phase curves

As said previously, the atmosphere prototype that we used in this study is most certainly too simple compared to any of the realistic compositions that we can expect. Within the limited sample of known cases, the Venus atmosphere may represent an analog of some hot terrestrial exoplanets atmospheres. On Venus, the presence of clouds and the extremely dense atmosphere would make detecting the phase curve difficult, as all infrared wavelengths probe high altitude levels with little day-night variation. But even in such a case, if the signal-to-noise ratio were good enough, it would be possible to infer the presence of a dense atmosphere from not detecting the phase variation, as attempted by Seager & Deming (2009). The ability to distinguish planets with dense atmospheres from airless bodies is essential for addressing the questions of atmosphere formation and survival at short orbital distances.

In future simulations, we will explore in detail the phase curve signatures of more complex atmospheres, containing a mixture of absorbers, condensable species, and aerosols.

4.7. Null detection and albedo

For nontransiting planets, the absence of photometric modulation can come from the presence of a dense atmosphere, to a very low inclination (see next section), or to a Bond albedo close to 1 (highly reflective planet). Such a high surface albedo is known only for ices, which can be ruled out for the starlit hemisphere of hot planets. Clouds can possibly produce a high Bond albedo (close to 0.7 for Venus), but cloud cover is also a signature of an atmosphere.

We should, however, remain careful about the properties of planetary surfaces of unknown composition. For instance, the phase variation of the transiting planet Kepler 10b has been observed in the 0.4−0.9 μm band of Kepler (Batalha et al. 2011). This phase variation cannot be attributed to the thermal emission of the planet, which remains too low even for a null albedo. To fit the phase variation of this exoplanet when assuming no atmosphere and Lambertian scattering requires an albedo over 0.5 (possibly close to unity). One explanation could be that the Lambertian approximation cannot be applied at optical wavelength and that the actual scattering phase function peaks towards the star, hence towards the observer at the maximum of the phase curve. For extremely hot planets like Corot 7b and Kepler 10b, the existence of a surface magma ocean underneath a tenuous transient atmosphere is possible in the substellar area (Léger et al. 2011). The optical and reflective properties of such a surface are unknown.

In summary, a Bond albedo close to unity can be responsible for an undetected thermal emission but no known solid, non-icy planetary surface has this property.

4.8. Null detection and inclination

With an inclination close to 0°, the observer always sees half of the day- and night-sides, and the day-night brightness temperature contrast does not produce phase variations. In this configuration, periodic variations in the thermal emission can only be due to seasonal changes, because of either an obliquity or an eccentric orbit. Tidally evolved planets on circular orbits would not exhibit seasonal variations.

Fig. 13 Thermal phase curves at 3 wavelengths for a 60° (left) and a 30° (right) inclination. Photon fluxes are calculated for a 1.5 m telescope. Graphs on the left are the same as in Fig. 7.

The probability of a very low inclination is, however, low (only 6% and 1.5% of randomly oriented systems have an inclination lower than 20° and 10°, respectively) and the sensitivity to inclination is weak for inclinations over 20°. This weak sensitivity to the inclination is due to the small variation of cos(θ) (where θ is the zenith angle) over a wide angular region around the substellar point. As a result, a large fraction of this hot substellar region contributes to the emission at the maximum phase even for inclinations down to about 20°. Figure 13 compares phase curves obtained at 60° and 30° inclinations. In addition, the range of possible inclinations can often be constrained. When the measured projected rotation velocity of the star Vsin(i) is measured, a minimum value for i can be deduced from the possible range of values for V for the type and age of the star, and V can also be measured for active stars by the photometric variability due to the stellar spots. This was done for the planet-hosting star GJ 876 (Correia et al. 2010). This reasoning assumes that the i in Vsin(i) is close to the orbital inclination. This can be very wrong for hot Jupiters, but it might be more reasonable for other systems. In multiple systems the orbital inclinations can be constrained directly when fitting dynamical orbits to the radial velocity data, as done for planets b and c, also around GJ 876 (Correia et al. 2010). In the GJ 581 system, where the hot terrestrial planet e would be a fine target for phase curve observation, inclinations lower than 60° would not be consistent with a stable system, due to the resulting high planetary masses and associated planet-planet interactions (assuming that mutual inclinations between the planets remains small). Astrometric measurements could also be used to determine the inclination, or at least constrain it, although the short-period planets we consider here are not easy targets for astrometry.

The number of configurations in which a low inclination prevents us from inferring any atmospheric properties should thus remain low, and these cases are possible to identify.

4.9. Validity of the Lambertian approximation

In each latitude-longitude cell, the flux F_λ at the top of the atmosphere for each band is an output of the GCM. To compute the flux received by an observer, we need the angular distribution of specific intensities. We assume an isotropic (Lambertian) distribution of intensities (I_λ = F_λ/π). This approximation neglects limb darkening (or limb brightening, which is possible when a temperature inversion exists). To quantify this approximation, we used a detailed 1D line-by-line radiative transfer code. We considered an uniform atmosphere (same P(z) and T(z) everywhere) and computed the specific intensities at the top of the atmosphere with a resolution of 5°. We used these specific intensities to directly compute the fluxes received by a remote observer, and we also integrated them to calculate the top-of-the atmosphere fluxes. From these top-of-the atmosphere fluxes, we calculated the approximated disk-integrated fluxes received by the distant observer, assuming an isotropic distribution of intensities. Comparison of the two disk-integrated fluxes showed that the error due to limb-darkening is below 5% in all our bands.

It is, however, likely that the error could be more significant for a nonuniform atmosphere. Indeed, the hot substellar region is responsible for the largest phase curve modulations. At some wavelengths, this region should be fainter than it is in our modeling when it is near the edge of the apparent planetary disk. At low inclinations (<30°), this could lower the amplitude of the variations. At higher inclinations, the amplitude of the modulation should not be significantly affected, but the variation should be less sharp at wavelengths where the visibility of the hot spot decreases when approaching the edges of the planetary disk. Doing the full line-by-line spectrum synthesis in 3D would be too time-consuming but, for future studies, we are developing a 3D band model at the resolution of the GCM.

4.10. Transiting planets

Phase curves can also of course be measured for systems in a transit configuration. In this case, the absolute flux and not just the variations can be obtained, as the flux of the star alone is inferred from the secondary eclipse. In addition, the day-side emission spectrum of the planet can be measured by comparing the star + planet emission just after or before the secondary eclipse with the emission during the secondary eclipse. However, as in the nontransiting configuration, the phase curve provides precious information on the atmospheric circulation, with higher accuracy thanks to the ability to measure the absolute planetary flux. In addition, the variation spectrum is not limited by the transit duration. For instance, the transit duration for the system studied in this paper is 1 h 30 min, while the phase curve can be sampled with seven exposures of  ~28 h to cover the 8-day orbit. The phase curve and secondary eclipse obtained on Kepler 10b illustrate this point nicely. Although in this case it is the scattered light that has been observed, Fig. 13 of Batalha et al. (2011) shows that the amplitude of the phase curve is better measured than the depth of the secondary eclipse.

5. Conclusions

By using 3D climate simulations of pure CO₂ atmospheres on a tidally-locked, moderately hot terrestrial planet, we have studied the influence of an atmosphere on the multiwavelength thermal phase curve. We described the effect of increasing pressures on the spectral and spatial distribution of the thermal emission and on the observed phase curves.

We have shown that the variation spectrum (the peak amplitude of the observed modulation due to the planet as a function of wavelength), which can be measured for a nontransiting planet, is shaped by molecular bands, just like the emission spectrum. This is because the day-night temperature contrast changes with altitude, and different wavelengths probe different altitudes. Signatures of atmospheric molecules can thus be extracted from spectro-photometric observations of nontransiting planets.

The existence of an atmospheric window probing the near-surface layer is, however, necessary for obtaining enough phase curve variations. A high atmospheric column density, a mixture of absorbers covering the whole wavelength range of the observations, or clouds, could prevent such windows and restrict the characterization to the sole inference of a dense atmosphere.

The 4.3 μm band of CO₂ appears in the variation spectrum at all pressures and is thus a promising way to detect CO₂. However, measuring the phase curve at this wavelength would require a hotter planet than the one studied here, in order to reach observable planet/star contrasts and SNR.

On an 8-day synchronous planet around an M 3 dwarf, a pure CO₂ atmosphere is stable against a condensation collapse but only for mean surface pressures above a certain threshold, estimated to be in the 0.1−0.5 bar range. In a future study, we intend to refine the actual value of this minimum pressure and its dependence upon parameters such as surface albedo, stellar type, orbital period, and planetary radius.

A photon-noise limited instrument of only 1.5 m with an extremely stable detector able to perform photometry at a level of 10^-5 over several days simultaneously in several spectral bands could in theory detect modulation by a nontransiting planet with a precision allowing the atmosphere to be characterized. The level of characterization that can be achieved depends on the nature and density of the atmosphere.

The wavelength-dependent phase curves computed in this study can be used to assess the observability of nontransiting terrestrial exoplanets by taking intrinsic stellar variability into account. The stellar variability and photometric stability of the instrument are certainly the main challenges when measuring the phase curves of nontransiting and spatially-unresolved terrestrial exoplanets.

Acknowledgments

F.S. acknowledges support from the European Research Council (Starting Grant 209622: E₃ARTHs).

References

All Tables

All Figures

	Fig. 1 The orbit and phases as seen by the observer. The inclination is set to 60°. The phases shown here are numbered from 1 to 12 and correspond to the phases presented in Figs. 3 to 5.
In the text

	Fig. 2 Equatorial temperature (top) and horizontal wind maps (bottom) at the surface (left) and 40 km (right). Temperatures are given for 0.1, 1 and 10 bar; wind maps are shown for the 10 bar case only. The displacement of the hot spot, compared to the substellar point (longitude = 0°, latitude = 0°) can be seen on the surface plots.
In the text

	Fig. 3 Planetary emission maps as a function of wavelength and phase. Phase numbers refer to Fig. 1.
In the text

	Fig. 4 Disk-integrated flux as a function of wavelength and phase.
In the text

	Fig. 5 Brightness temperature maps as a function of wavelength and phase (top) and gas temperature maps as a function of altitude and phase (bottom).
In the text

	Fig. 6 Planetary signal as a function of time. The 8-day orbits is divided in 7 exposures. The number of photons is given per exposure, per band, for a 6 m telescope and a 10 pc target. The dotted line represents the amplitude of the stellar photon noise. Wavelengths are indicated in the upper left corner of each plot.
In the text

	Fig. 7 Same as Fig. 6 but for a 1.5 m telescope.
In the text

	Fig. 8 Same as Figs. 6 and 7 but the photon flux is given here as the planet/star contrast ratio.
In the text

	Fig. 9 Planetary spectrum vs. variation spectrum. Circles, squares, and diamonds correspond the planetary spectrum at 3 different phases indicated at the top of the graph (phases 6, 3, and 12 in Fig. 1). The dashed line shows the amplitude of the 7-point phase curves as a function of wavelength (i.e. variation spectrum).
In the text

	Fig. 10 Contribution of the reflected light to the phase curve with no atmosphere. The 7-point phase curves are computed with emission only (black circles) and with emission+reflection (diamonds). Band central wavelengths are indicated in the upper left corner of each plot.
In the text

	Fig. 11 Condensation of CO2 in the 0.1 bar case. Color contours indicate the surface temperature, with the same color scale as in previous temperature maps. Contour lines indicate the region where CO2 is supersaturated at the surface, at 4 m, 10 m and 230 m. Levels above 230 m are not saturated. The 0°−0° lon-lat coordinate corresponds to the substellar point.
In the text

	Fig. 12 Amplitude of the phase curve for an airless planet as a function of stellar mass, orbital period and wavelength. The albedo of the planet is 0.2 and the inclination is 60°. Grey circles indicate the case modeled in this paper. Stellar radii and Teff are from Baraffe et al. (1998).
In the text

	Fig. 13 Thermal phase curves at 3 wavelengths for a 60° (left) and a 30° (right) inclination. Photon fluxes are calculated for a 1.5 m telescope. Graphs on the left are the same as in Fig. 7.
In the text