© ESO, 2013

1. Introduction

The number of discovered exoplanets is growing at a very fast pace, reaching more than seven hundred in May 2013¹. Nearly one-third of them transit in front of their host stars, a lucky circumstance that allows us to access their radius through simple geometrical assumptions. Nevertheless, the great majority of these planets lack a complete characterization, which should include the knowledge of their internal structure, atmospheric composition, and evolutionary history (Seager & Deming 2010). This is even more limiting if one takes into account that most extrasolar planets known so far have physical properties that are completely different from those found in our solar system, i.e. we cannot rely on “easy” analogies (Zhou et al. 2012).

Most discoveries of the early years (1995−2007) were strongly biased towards “hot Jupiter”-type planets, i.e., gas giants orbiting at a few stellar radii from their host star. Thanks to dedicated space missions, such as CoRoT (Baglin et al. 2006) and Kepler (Borucki et al. 2010), we are now accessing a broader region of the parameter space, probing smaller planetary radii and masses, and detecting planetary systems hosted by a wider variety of parent stars. Not surprisingly, the observational landscape is getting even more complicated, with wholly new – and sometimes unexpected – classes of planets appearing.

“Super-Earths” (2 ≲ M_p ≲ 10 M_⊕) and “Neptunian” planets (15 ≲ M_p ≲ 50 M_⊕) are now puzzling theoreticians, revealing a complex diversity that cannot be easily explained by models or constrained by observations (Haghighipour 2011). Super-Earths and Neptunes were once hypothesized to be separate and well defined classes of rocky and icy planets, respectively. Now we realize that their measured densities range from ρ_p = 0.27 (Kepler-18d; Cochran et al. 2011) to ~10 g   cm^-3 (Kepler-10b; Batalha et al. 2011), with a significant overlap between super-Earths and Neptunes (Wright et al. 2011). In this mass range, average densities are unable to put a firm constraint on the inner structure of the planet, and theoretical models can lead to several degenerate solutions (e.g., Adams et al. 2008). One of the missing key quantities is the atmospheric scale height H of the planet, defined, for a well-mixed isothermal atmosphere in equilibrium, as $\mathit{H}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{b}}{\mathit{T}}_{\mathrm{eq}}}{\mathit{\mu g}}\mathit{,}$(1)where k_b is the Boltzmann constant, T_eq the equilibrium temperature at its surface, μ the mean molecular weight of the atmosphere, and g the planet surface gravity.

An emblematic case is that of GJ1214b (Charbonneau et al. 2009), for which allowed scenarios include a “mini-Jupiter” with a small solid core and large envelopes of primordial H/He, a Neptune-like planet with an atmosphere made of sublimated ices, a “water world”, or a rocky super-Earth with an outgassed atmosphere (Rogers & Seager 2010). Newer examples of such degeneracy are HAT-P-26b (Hartman et al. 2011) and GJ3470b (Bonfils et al. 2012). Additional information other than ρ_p is required to break the model degeneracy, and an effective way to do it is to search for spectral signatures coming from the planetary atmosphere to probe its composition, or to put constraints on μ (Seager & Deming 2010).

Additional information can come from transmission spectroscopy, a technique based on the observation of exoplanetary transits at different wavelengths, searching for changes of its effective radius R_p as a function of λ owing to absorption and/or scattering processes undergone by stellar light after traveling through the atmospheric limb (Sing et al. 2011). Of course, low-resolution transmission spectroscopy can also be performed at a higher efficiency by gathering photometric light curves through different passbands and investigating the resulting R_p(λ) dependence (e.g., de Mooij et al. 2012).

While transmission spectroscopy is most frequently carried out to search for atomic or molecular absorption features, useful information can be also inferred from the detection of scattering processes affecting the continuum. The most straightforward one to interpret is Rayleigh scattering, due to molecules such as H₂ (Lecavelier des Etangs et al. 2008b) or to haze from condensate particles (Fortney 2005; Lecavelier des Etangs et al. 2008a). Rayleigh scattering has a typical signature: its steep λ^-4 dependence, which translates into a larger R_p on the blue side of the optical spectrum. Under simplified assumptions, we can analytically approximate the scale height H (and hence μ, once a T_eq is assumed) from the steepness of R_p(λ) where Rayleigh scattering is dominant (Lecavelier des Etangs et al. 2008a, among others). On the other hand, detailed atmospheric models are available, but at the price of increased complexity and number of free parameters (Howe & Burrows 2012).

GJ3470b has been defined a “hot Uranus” planet by its discoverers, who first detected the 3.33-day radial velocity (RV) modulation of its M1.5V host star with HARPS and then caught its transits through ground-based photometry (Bonfils et al. 2012). Very recently, a follow-up analysis based on Spitzer, WIYN-3.5m and Magellan data (Demory et al. 2013) has resulted in updated values for the mass and radius of GJ3470b: M_p = 13.9 ± 1.5 M_⊕ and R_p = 4.83 ± 0.21 R_⊕, making it one of the less dense low-mass planets known with ρ_p = 0.72 ± 0.13 g   cm^-3. Showing this peculiarity, having a total mass at the boundary between super-Earths and Neptunes, GJ3470b stands out as an ideal case to test the present planetary structure and evolution theories, which predict a significantly extended envelope of primordial H and He for it (Demory et al. 2013; Rogers & Seager 2010) and thus a large scale height. From the observational point of view, GJ3470b is promising thanks to its relatively large transit depth (R_p/R_⋆)² ≃ 0.006, being hosted by a small star (M_⋆ ≃ 0.54 M_⊙, R_⋆ ≃ 0.57 R_⊙; Demory et al. 2013). Nevertheless the level of chromospheric activity of GJ3470 appears to be low, given its spectral type, and his <0.01 mag long-term photometric variability in the I_c band (Fukui et al. 2013).

In this paper we present new photometric data from LBT, consisting of two simultaneous transit light curves of GJ3470b gathered through passbands at the extremes of the optical spectrum (λ_c = 357.5 and 963.5 nm). Our primary aim was to search for signatures of Rayleigh scattering in the planetary atmosphere. We describe the instrumental setup and observing strategy in Sect. 2, and our data reduction techniques in Sect. 3. Next we illustrate how we extracted the photometric parameters of GJ3470b from our light curves (Sect. 4) including a correction on k = R_p/R_⋆ due to the presence of unocculted spots. Finally, we discuss our results in Sect. 5.

2. Observations

We observed a full transit of GJ3470b at LBT (Mt. Graham, Arizona) in the night between 2013 Feb. 16 and 17, under a clear, photometric sky. The Large Binocular Camera (LBC; Giallongo et al. 2008) consists of two prime-focus, wide-field imagers mounted on the left and right arms of LBT, and optimized for blue and red optical wavelengths, respectively. We exploited their full power to gather two simultaneous photometric series, setting the U_spec filter on the blue channel (a filter with Sloan u response but having an increased efficiency, centered at λ_c = 357.5 nm) and the F972N20 filter in the red channel (an intermediate-band filter centered at λ_c = 963.5 nm, originally designed for cosmological studies)². These choice of passbands have been made to maximize the wavelength span of our observations, while at the same time avoiding most of telluric lines in the red channel.

We set a constant exposure time of 60 s on both channels. We chose to read out only a 2300 × 1800 pixel window from the central chip of each camera (chip#2) in order to decrease the technical overheads and optimize the efficiency of the series, achieving a 60% duty-cycle on average. The resulting 7.6′ × 6.7′ field of view (FOV) was previously tailored to include a suitable set of reference stars, required to perform differential photometry with a signal-to-noise ratio (S/N) as high as possible in both channels. Stellar point spread functions were defocused to ~7′′ FWHM (full width at half maximum; 31 physical pixels) in U_spec and ~13′′ (58 pix) in F972N20 to avoid saturation and to minimize intra-pixel and pixel-to-pixel systematic errors (Nascimbeni et al. 2011a). Autoguiding is not allowed on LBC during a photometric series; instead, the telescope was passively tracking throughout the series, drifting by about 80 pixels from the first to the last frame.

A total of 134 U_spec frames and 146 F972N20 frames were secured from 3:22 to 7:06 UT, covering also 55+55 minutes of off-eclipse photometry before and after the external contacts of the transits.

3. Data reduction

Both sets of frames were bias-corrected and flat-fielded using standard procedures and twilight master flats. Then we extracted the U_spec and F972N20 light curves of GJ3470 using STARSKY.

Fig. 1 Left plot: light curve of GJ3470b gathered with the blue channel of LBC in the Uspec band, plotted with the original sampling cadence (upper panel; the black line corresponds to the best-fit model adopted in Table 1) and residuals from the best-fit model (lower panel). Right plot: same as the left plot, but in the F972N20 band of the red channel of LBC. The two light curves are simultaneous.

STARSKY is a photometric pipeline optimized to perform high-precision differential aperture photometry over defocused images, originally developed for the TASTE project (The Asiago Search for Transit timing variations of Exoplanets; Nascimbeni et al. 2011a,b). Its prominent feature is an empirical algorithm to choose the aperture radius and to weight the flux from the comparison stars in such a way that the off-transit scatter of the target star is minimized (Nascimbeni et al. 2013). In a second step, STARSKY searches for linear correlations between the differential magnitude and a set of standard external parameters, including the (x,y) position of the star, its FWHM, sky background, and time t. In our case, a 22.8-pixel (5.1′′) aperture was selected as the optimal one for the U_spec data set, followed by a decorrelation against FWHM and t. For the F972N20 series, the optimal aperture diameter was 40.0 pixels (9′′), followed by a decorrelation against t. It is remarkable that we did not detect any significant (x,y)-dependent systematic effect in either the blue channel or the red one, suggesting that the “extreme” defocusing strategy we applied is effective against most residual flat-field and intrapixel errors.

The resulting light curves are plotted in Fig. 1 in their original cadence. Their overall root mean square (RMS) scatter is 1.19 mmag (U_spec) and 0.28 mmag (F972N20) over a mean cadence of 95 s, only slightly larger than the theoretical values estimated by standard “white noise” formulae: 1.08 mmag and 0.24 mmag, respectively (Howell 2006). It is worth noting that the F972N20 series is, to our knowledge, the most precise light curve ever obtained from a ground-based facility, breaking by a narrow margin the previous record set by Tregloan-Reed & Southworth (2013) at the New Technology Telescope (NTT). The prospects from these achievements will be discussed further in Sect. 5.

4. Data analysis

4.1. Light curve modeling

The JKTEBOP code³ version 28 (Southworth et al. 2004) was employed to fit a transit model over the LBC light curves. We chose to derive the uncertainties over the best-fit parameters with a Monte Carlo bootstrap algorithm (Southworth et al. 2005), to avoid underestimation due to correlations between parameters and residual systematic noise (also called “red noise”; Pont et al. 2007). A quadratic law for modeling the stellar limb darkening (LD) is adopted (Claret 2004).

We first fitted the F972N20 curve leaving five free parameters: the radius ratio k = R_⋆/R_p, the sum of the planet and star radii scaled by the orbital semi-major axis Σ = (R_⋆ + R_p)/a, the orbital inclination i, the central time T₀, and the linear LD coefficient u₁. Following a common practice to improve the robustness of the fit (Southworth 2008), we did not fit for the quadratic LD coefficient u₂, because it is strongly degenerate with u₁. Even fitting for the weakly correlated parameters c₁ = 2u₁ + u₂ and c₂ = u₁ − 2u₂ (as suggested by Holman et al. 2006) or other linear combinations of u₁ and u₂ leads to unphysical results due to overfit. Instead, u₂ was injected into the Monte Carlo bootstrap as a Gaussian prior, extracting its theoretical value from the tables by Claret et al. (2012), interpolated by adopting the stellar atmospheric parameters published by Demory et al. (2013) and taking their uncertainties on T_eff, log g and [Fe/H] into account. As the F972N20 band is not tabulated by Claret et al. (2012), we adopted the SDSS z′ values basing on the fact that z′ and F972N20 overlap with each other in a spectral region where the specific intensity I_ν(θ) is nearly devoid of spectral features (θ being the angle between the observer and the normal to the stellar surface). Moreover, for our given set of atmospheric parameters, the quadratic term u₂ is nearly insensitive to λ, oscillating between 0.260 and 0.366 over the full optical-NIR range between Sloan u to 2MASS J, with a mean and RMS of 0.30    ±    0.03. After a first best-fit solution is found on real data, we ran a bootstrap analysis on 20   000 resampled light curves having the same statistical properties as the real one. The final best-fit values and the asymmetric errors on each fitted quantity were then estimated from the 50th, 15.87th, and 84.13th percentiles of the bootstrapped distribution (Table 1, second column). We emphasize that the fitted value of u₁ = 0.25    ±    0.04 is perfectly consistent with the theoretical estimate u_1,th = 0.275 for the z′ band by Claret et al. (2012), supporting our previous assumption u₂(z′) ≃ u₂(F972N20). As a further cross-check, we verified that even by setting the input u₂ at the aforementioned extremal values 0.26 and 0.37, the resulting best-fit parameters change by less than 1-σ with respect to our adopted solution.

The U_spec light curve has a much lower S/N compared to the F972N20 one. We then improved the robustness of the fit by fixing the value of i (which is a purely geometrical parameter, not dependent on λ) to what was found on the F972N20 series (88.1° ± 0.3), by injecting a Gaussian distribution with the same mean and standard deviation into the bootstrap analysis. The same is done for T₀. Instead, Σ (and so derived quantities such as a/R_⋆) was left free, since a small dependence of these variables on λ is expected when R_p(λ) itself is not constant. The value of both u₁ and u₂ are fixed to their Claret et al. (2012)u′ theoretical Gaussian priors as above. The resulting best-fit values for the two fitted parameters (Σ and k) are reported in the third column of Table 1. Other quantities of interest, such as the impact parameter b and the full duration T₁₄ (Winn et al. 2009) are extracted for both LBC light curves from the bootstrapped distribution, assuming zero eccentricity (e = 0), as suggested by the RV measurements (Bonfils et al. 2012).

We note that all the best-fit parameters we estimated on both LBC curves are in excellent agreement with each other and with the Demory et al. (2013) values (Table 1) within their errorbars. The only significant exception is k, which is offset respectively by $\frac{\mathit{k}\mathrm{\left(}{\mathit{U}}_{\mathrm{spec}}\mathrm{\right)}\mathrm{-}\mathit{k}\mathrm{(}\mathrm{4.5}\mathit{\mu }\mathrm{m}\mathrm{)}}{\mathrm{\Delta }\mathit{k}\mathrm{\left(}{\mathit{U}}_{\mathrm{spec}}\mathrm{\right)}}\mathrm{\simeq }\mathrm{3.6}\hspace{0.17em}\mathit{\sigma ,}\hspace{0.17em}\hspace{0.17em}\frac{\mathit{k}\mathrm{\left(}\mathit{F}\mathrm{972}\mathrm{\right)}\mathrm{-}\mathit{k}\mathrm{(}\mathrm{4.5}\mathit{\mu }\mathrm{m}\mathrm{)}}{\mathrm{\Delta }\mathit{k}\mathrm{\left(}\mathit{F}\mathrm{972}\mathrm{\right)}}\mathrm{\simeq }\mathrm{-}\mathrm{6.4}\hspace{0.17em}\mathit{\sigma }\mathit{.}$(2)This offset is not driven by a different assumption on the orbital inclination (or, equivalently, on the impact parameter), as our best-fit value is in perfect agreement with found by Demory et al. (2013) and later adopted by Fukui et al. (2013) as a prior in their analysis. As a further check, we reran our analysis, adopting a Gaussian prior and obtaining statistically indistinguishable values of k for both light curves (Δk/σ(k) ≤ 0.3). As suggested by the anonymous referee, we also redid a full analysis of both light curves by injecting Gaussian priors on c₁ and c₂ defined as above and leaving everything else unchanged. Again, the resulting best-fit parameters and their associated errors are fully consistent with those listed in Table 1, meaning that correlation between the LD parameters cannot explain the measured discrepancy between k(U_spec) and k(F972). In the next section we apply a correction to our k values to take the presence of unocculted starspots on the GJ3470 photosphere into account. Then in Sect. 5 we discuss the significance and origin of this chromatic signature.

The extremely high precision on T₀ we achieved on the F972N20 transit (about 9 s) allowed us to compute a greatly improved orbital ephemeris of GJ3470b, by combining our measurement with the two T₀ published by Demory et al. (2013) and fitting a linear ephemeris by ordinary weighted least squares: $\begin{array}{ccc}{\mathit{T}}_{\mathrm{0}}\mathrm{\left(}{\mathrm{BJD}}_{\mathrm{TDB}}\mathrm{\right)}& \mathrm{=}& \mathrm{(}\mathrm{2}\hspace{0.17em}\mathrm{456}\hspace{0.17em}\mathrm{340.72559}\mathrm{\pm }\mathrm{0.00010}\mathrm{)}\\ & & \mathrm{+}\mathrm{(}\mathrm{3.336}\hspace{0.17em}\mathrm{649}\mathrm{\pm }\mathrm{0.000}\hspace{0.17em}\mathrm{002}\mathrm{)}\mathrm{·}\mathit{N}\end{array}$(3)where the integer epoch N is the number of transits elapsed from our LBC observation. The combination of time standard and reference frame here adopted is BJD_TDB (barycentric Julian Day, barycentric dynamical time), as we follow the prescription by Eastman et al. (2010).

4.2. Correction for unocculted starspots

Starspots affect the stellar flux variations in a transit light curve in two ways. Spots that are not occulted by the planet produce a decrease in the out-of-transit flux of the star, which is the reference level to normalize the transit profile. Spots that are occulted during the transit produce a relative increase in the stellar flux because the flux blocked by the planetary disk is lower than in the case of the unperturbed photosphere. Therefore, if the star is spotted, and f_o and f_u are the spot filling factors of the occulted and unocculted photosphere, respectively, the relative radius of the planet k_a derived from the fitting of the transit profile, is related to the true relative radius k by ${\mathit{k}}_{\mathit{a}}\mathrm{\simeq }\mathit{k}\left[\mathrm{1}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{A}}_{\mathit{\lambda }}\mathrm{(}{\mathit{f}}_{\mathrm{o}}\mathrm{-}{\mathit{f}}_{\mathrm{u}}\mathrm{)}\right]\mathrm{,}$(4)where A_λ ≡ 1 − F_s(λ)/F_⋆(λ), with F_⋆(λ) and F_s(λ) the emerging fluxes from the unperturbed photosphere, and the spotted area considered as a photospheric model computed at the mean spot temperature, respectively.

Fig. 2 Upper panel: reconstructed transmission spectrum of GJ3470b. The Uspec and F972N20 data points are extracted from our LBC light curves (blue circles), the g′RcIcJ points from Fukui et al. (2013) (green circles), and the 4.5  μm point from Demory et al. (2013) (red circle). The model plotted with a dark green line is scaled from a Howe & Burrows (2012) model computed for a 10-M⊕ cloud-free, haze-free planet with Teq = 700 K and a pure water (H2O) atmosphere. The open squares represent the model points integrated over the instrumental passbands. Second panel: same as above, but with a cloud-free, haze-free, metal-poor atmosphere dominated by H and He, computed by the same authors (orange line). Abundances are scaled from the solar ones (Z = 0.3 × Z⊙). Third panel: same as above, but adding the contribution of a scattering haze from Tholin particles at 1 μbar. Lower panel: instrumental passbands employed for each data point.

Bonfils et al. (2012) suggest that GJ3470 is not a very active star according to the low projected rotational velocity (≲2 km   s¹) and the H_α (6562.808 Å) line observed in absorption, which are signatures of a mature star older than ~300 Myr. The same authors suggest that GJ3470 is younger than 3 Gyr, since its kinematic properties match those of the young disk population. Indeed, Fukui et al. (2013) performed a photometric monitoring of GJ3470, finding a peak-to-valley variability in the I_c band on the order of 1% in an interval of ~60 days. Assuming that such variability is due to the presence of starspots, these have a filling factor f of ~0.01 (Ballerini et al. 2012). In Eq. (4), the effects of the unocculted and occulted spots tend to compensate to each other, although the filling factor of the occulted spots may generally be higher because starspots tend to appear at low or intermediate latitudes in Sun-like stars. The filling factor of the out-of-transit spots f_u can be assumed to be constant because the visibility of those spots is modulated on timescales much longer than that of the transit, i.e., on the order of the stellar rotation period, which is generally several days. Instead, the filling factor of the occulted spots f_o is in general a function of the position of the center of the planetary disk along the transit chord. A variable spot distribution in the occulted area gives rise to bumps in the observed transit light curve, whereas the lack of obvious bumps in our data (Fig. 1) suggests a more or less continuous background of spots along the transit chord. In the observed case, f_o can be assumed to be almost constant during the transit because otherwise we would have detected the individual spot bumps.

The maximum deviation in the relative radius k due to starspots can be computed in the unphysical hypothesis that all the spots are concentrated in the occulted area (f_o = 0.01 and f_u = 0) or outside the occulted area (f_o = 0 and f_u = 0.01). Using the grid given by Ballerini et al. (2012), where GJ3470 is very close to their case #4, the maximum correction of the relative radius in our passbands of interest is $\begin{array}{c}\\ \mathrm{\Delta }{\mathit{k}}_{\max }\hspace{0.17em}\mathrm{\left(}{\mathit{U}}_{\mathrm{spec}}\mathrm{\right)}& \mathrm{=}& \mathrm{0.00540}\\ \mathrm{\Delta }{\mathit{k}}_{\max }\hspace{0.17em}\mathrm{\left(}{\mathit{g}}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& \mathrm{0.00538}\\ \mathrm{\Delta }{\mathit{k}}_{\max }\hspace{0.17em}\mathrm{\left(}{\mathit{R}}_{\mathrm{c}}\mathrm{\right)}& \mathrm{=}& \mathrm{0.00529}\\ \mathrm{\Delta }{\mathit{k}}_{\max }\hspace{0.17em}\mathrm{\left(}{\mathit{I}}_{\mathrm{c}}\mathrm{\right)}& \mathrm{=}& \mathrm{0.00500}\\ \mathrm{\Delta }{\mathit{k}}_{\max }\hspace{0.17em}\mathrm{\left(}\mathit{F}\mathrm{972}\mathit{N}\mathrm{20}\mathrm{\right)}& \mathrm{=}& \mathrm{0.00500}\\ \mathrm{\Delta }{\mathit{k}}_{\max }\hspace{0.17em}\mathrm{\left(}\mathit{J}\mathrm{\right)}& \mathrm{=}& \mathrm{0.00391}\\ \mathrm{\Delta }{\mathit{k}}_{\max }\hspace{0.17em}\mathrm{\left(}\mathrm{4.5}\mathit{\mu }\mathrm{m}\mathrm{\right)}& \mathrm{=}& \mathrm{0.00355}\mathrm{.}\end{array}$(5)These corrections should be only considered an upper limit to the actual impact of activity on our estimated k(λ). For this reason, in the following discussion and plots, we adopt the uncorrected values of k(λ) and R_p(λ), unless noted otherwise, and discuss separately the limit case when Δk_max is applied.

A Δk_max correction ranging from 0.00355 to 0.00540 could appear quite large since it is an order of magnitude more than our best errors on k (~0.0005 in F972N20 and 4.5  μm). One has to consider, however, that the absolute scaling of this correction (which is related to the “solid” radius of the planet R_p,0 measured at H = 0) is of secondary interest for our study. Instead, the differential effect of the Δk_max correction is crucial, because it changes the slope of the Rayleigh scattering absorption and the shape of the overall transmission spectrum. It is worth noting that this may not hold when gathering transits at different epochs. In such a case, due to starspot evolution or activity cycles, the filling factors f_o and f_u may change, adding an unknown systematic offset between non-simultaneous measurements. However, our main result is based on simultaneous measurements (in U_spec and F972N20), so they are not affected by such an effect.

5. Discussion and conclusions

5.1. The transmission spectrum of GJ3470b

If we examine the transmission spectrum of GJ3470b as reconstructed by our LBC measurements and from those by Demory et al. (2013) and Fukui et al. (2013) (Fig. 2, large circles), it appears evident that no constant value of k could explain the data (χ² = 67.8 for 5 degrees of freedom, corresponding to a reduced value of ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}\mathrm{=}\mathrm{13.5}$). The color dependence of k is then significant, with a steep increase toward the blue side of the spectrum, which is usually associated with scattering processes. How can we interpret this spectrum k(λ)? A detailed atmospheric modeling of GJ3470b is beyond the scope of this paper. Instead, we rescaled the model transmission spectra computed by Howe & Burrows (2012) by taking the grid point of their simulations closer to GJ3470b (T_eq = 700  K, M_p = 10  M_⊕) and investigating which class of atmospheric composition best matches the observed data points and which ones can be excluded.

The most striking result is that the large observed variations of k(λ) are totally incompatible with all model sets by Howe & Burrows (2012) having an atmosphere of high mean molecular weight  μ, such as those made of pure H₂O (Fig. 2, upper panel) or pure methane (CH₄). All those models predict a very flat spectrum owing to their tiny scale height and are rejected by the available data at >10σ. The same conclusion can be drawn by considering the LBC data alone, since k(U_spec) and k(F972N20) differ by 6σ.

On the other hand, we investigated the fitness of cloud-free, haze-free atmospheric models of low μ, which is largely dominated by H and He and with a large scale height H. This is the most probable atmospheric composition expected for GJ3470b following current theories of planetary interiors (Demory et al. 2013; Rogers & Seager 2010). A metal-poor model atmosphere having a composition with Z = 0.3 × Z_⊙ rescaled from a solar abundance (Howe & Burrows 2012; second panel of Fig. 2) provides us with a much better fit on the red/IR passbands (F972N20, J, 4.5  μm). Nevertheless it fails to reproduce the steep rise observed on the blue side of the spectrum, especially for our U_spec measurement. That is, the scale height required to output a Rayleigh scattering from H₂ large enough to match our UV measurement would also give rise to infrared molecular features which are inconsistent with the reddest data points, in particular with the 4.5  μm one, which is in a spectral region particular sensitive to species, such as H₂O and CO. Similar models, but with the addition of clouds, perform even worse, flattening the spectrum at intermediate wavelengths without significantly improving the overall quality of the fit.

An additional source of scattering, other than H₂, seems to explain the data better. Fortney (2005) shows that a wide variety of condensate particles can be very efficient in this regard. Among them Tholin, for instance, is a polymer produced by UV photochemistry of CH₄, said to occur on Jupiter and giant exoplanets (Sudarsky et al. 2003). Indeed, the Howe & Burrows (2012) class of models with a metal-poor H/He-dominated atmosphere and high-altitude hazes from Tholin provide us with a satisfactory fit with ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}\mathrm{=}\mathrm{5.2}$ (black line in Fig. 2, third panel; here the model points integrated over the instrumental passbands are plotted in magenta for clarity). Even after the maximum correction for unocculted spots (Eq. (5)) has been applied on all k(λ) points, the “haze” model with ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}\mathrm{=}\mathrm{7.9}$ still fits the data set better than a flat function or any other class of Howe & Burrows (2012) models. We emphasize that Tholin was chosen here only as a test model to check whether the transmission spectrum of GJ3470b can be explained through an additional source of scattering due to condensate particles. Still, we should keep in mind that a large number of molecules are able to reproduce a similar result by means of the same physical process, and their characteristic spectral fingerprints are so subtle that they cannot be discerned in a low-resolution spectrum. Additional high S/N observations will be required to confirm this scenario, especially in the visible region at 400−900 nm (where the S/N of the Fukui et al. 2013 measurements is too low to constrain the slope of the Rayleigh scattering independently from our data) and in the near IR, for instance in the H and K spectral region, where molecular bands such as H₂O, CO, or CH₄ are expected to dominate the trasmission spectrum, giving rise to strong features (Fig. 2). At the estimated equilibrium temperature of GJ3470b, T_eq ~ 700 K, narrow alkali metal lines, such as NaI and KI, begin to appear but are still too weak to hope in a ground-based detection, even with the current top-class facilities. Nevertheless they become very prominent at 1000 K (Fig. 6 in Howe & Burrows 2012), so searching for them could put a tight and independent upper limit to the atmospheric effective temperature of GJ3470b.

The available data can also be interpreted through a simplified analytical approach, first developed by Lecavelier des Etangs et al. (2008a) under the usual approximation of a well-mixed and isothermal atmosphere in chemical equilibrium. Following this work, we assume a scaling law of index α for the planetary cross section σ = σ₀(λ/λ₀)^α. Under reasonable assumptions, the slope of the planetary radius as a function of wavelength λ can be expressed as a function of the scale height H (Eq. (1)) as dR_p/dlnλ = αH (Lecavelier des Etangs et al. 2008a). The expression for the planetary equilibrium temperature then reduces to $\mathit{\alpha }{\mathit{T}}_{\mathrm{eq}}\mathrm{=}\frac{\mathit{\mu g}}{\mathit{k}}\frac{\mathrm{d}{\mathit{R}}_{\mathrm{p}}}{\mathrm{d}\ln \mathit{\lambda }}\mathrm{·}$(6)If the main physical process involved is Rayleigh scattering and if we neglect atomic/molecular absorption, α =  −4 and the mean molecular weight can be estimated as $\mathit{\mu }\mathrm{=}\mathrm{-}\mathrm{4}\mathit{k}{\mathit{T}}_{\mathrm{eq}}{\left(\mathit{g}\frac{\mathrm{d}{\mathit{R}}_{\mathrm{p}}}{\mathrm{d}\ln \mathit{\lambda }}\right)}^{-1}\mathrm{·}$(7)We assume from Demory et al. (2013) a planetary surface gravity g = 5.75 m   s^-2 and an equilibrium temperature T_eq = (1 − A_b)^0.25(683    ±    27) K. As simple limiting cases, we can adopt two distinct values for the Bond albedo of GJ3470b: A_b = 0.3, which matches the albedo measured for Solar System icy giants (Uranus and Neptune); and A_b = 0, corresponding to a cloud-free, very dark surface as theoretically predicted for some classes of “temperate” Neptunian planets by Spiegel et al. (2010), and inferred from observations of the hot Neptune GJ436b (Cowan & Agol 2011). With these assumptions, the temperature becomes T_eq = 624 ± 25 K (A_b = 0.3) and T_eq = 683 ± 27 K (A_b = 0), respectively.

Fig. 3 Upper panel: same as Fig. 2, but zoomed on the optical region. The red line corresponds to the best weighted linear fit in the (lnλ,k) plane over all data points at λc < 1μm (ours and from Fukui et al. 2013) under the assumption of pure Rayleigh scattering; see text for details. Lower panel: instrumental passbands employed for each data point.

By fitting a straight line in the (k,lnλ) plane by weighted least squares and estimating the uncertainties through an ordinary bootstrap algorithm we find $\begin{array}{ccc}\mathit{k}& \mathrm{=}\mathrm{(}\mathrm{-}\mathrm{0.0076}\mathrm{\pm }\mathrm{0.0013}\mathrm{)}\ln \mathit{\lambda }\mathrm{+}\mathrm{(}\mathrm{0.123}\mathrm{\pm }\mathrm{0.008}\mathrm{)}& \mathrm{(}{\mathit{A}}_{\mathrm{b}}\mathrm{=}\mathrm{0.3}\mathrm{)}\\ \mathit{k}& \mathrm{=}\mathrm{(}\mathrm{-}\mathrm{0.0075}\mathrm{\pm }\mathrm{0.0012}\mathrm{)}\ln \mathit{\lambda }\mathrm{+}\mathrm{(}\mathrm{0.127}\mathrm{\pm }\mathrm{0.008}\mathrm{)}& \mathrm{(}{\mathit{A}}_{\mathrm{b}}\mathrm{=}\mathrm{0}\mathrm{)}\mathit{,}\end{array}$(8)which is equivalent to dR_p/dlnλ ≃  − 3000 ± 500 km for both cases once we adopt R_⋆ = 0.568 R_⊙ as the stellar radius from Demory et al. (2013) and propagate the uncertainties. From Eq. (7), and taking the errors on all variables into account, we derive a mean molecular weight of $\mathit{\mu }\mathrm{=}\mathrm{1.2}{\mathrm{0}}_{-0.17}^{\mathrm{+}\mathrm{0.24}}$ atomic mass units (amu; A_b = 0.3) or $\mathit{\mu }\mathrm{=}\mathrm{1.3}{\mathrm{2}}_{-0.19}^{\mathrm{+}\mathrm{0.27}}$ amu (A_b = 0) for GJ3470b, that is, a mean molecular weight consistent with that of a H/He-dominated atmosphere with subsolar chemical abundances. The error bars are largely dominated by the uncertainty on k. Our result is virtually independent of the correction Δk_max in Eq. (5), which leads to a nearly unchanged μ = 1.3 ± 0.25 amu (A_b = 0.3) when it is applied. We emphasize however that the last approach does not take into account that the k value at 4.5  μm seems to be inconsistent with a cloud-free H/He envelope. On the other hand, the (k,lnλ) fit also demonstrates that the increase in k (or R_p) at blue wavelengths is significant at 5.8σ, regardless of the underlying physical process. The significance is largely driven by our U_spec point, and this emphasizes the importance of the optical UV bands in exoatmospheric studies.

Fig. 4 Upper panel: synthetic 1-R⊕ transit (black line) injected in the residuals of our F972N20 light curve (red points and error bars), showing that a short-period Earth around an M1.5V star would have been detectable from our photometry. Lower panel: same as above, binned over 20-min intervals.

5.2. Pushing the limits of ground-based photometry

As a final note we highlight that the accuracy of our F972N20 light curve is close to the theoretical limit expected for a photometric series gathered with a 8.4-m telescope. It is interesting to investigate the consequences of this unprecedented accuracy on the ground-based detection of very small planets. A way to do it is to extract the F972N20 residuals from our best-fit model and to inject on them a synthetic transit corresponding to a 1-R_⊕ planet hosted by a star having the same radius of GJ3470, i.e. a M1.5V dwarf. Having employed real residuals to perform this simulation, we are confident that the instrumental, atmospheric, and astrophysical noise contributions are preserved. The result is plotted in Fig. 4, with the original 90-s sampling rate (upper panel) and binned over 20-min intervals to increase the S/N (lower panel). The transit signature is clearly detectable with high confidence (9σ), opening up new prospects for a future targeted search for short-period, Earth-sized transiting planets around M dwarfs without necessarily resorting to much more expensive, dedicated space missions.

The short-term photometric accuracy of the data presented here is probably not achievable on a longer term (night-to-night) baseline, owing to atmospheric and instrumental systematic drifts. This limitation, combined with the extremely low density of bright M dwarfs in the sky, prevents us from a “classical” transit search of terrestrial planets with a camera such as LBC. However, when a sample of low-mass planets hosted by red dwarfs and discovered by an RV survey is available, it will be possible to predict the instant of the inferior conjunction within a few hours of uncertainty and to search for planetary transits in that observational window. For a scaled-down version of GJ3470b, the a priori transit probability is around 8%, meaning that at least a dozen RV planets need to be followed up in order to detect one transit. Given the importance of these planets for the exoplanetary studies, thanks especially to their closeness to the habitable zone and to the relatively easy access to their atmospheric signatures, we suggest that a targeted search would be worth implementation.

Note added in proof. Simultaneously with this report, an analysis of optical and infrared transit photometry of GJ 3470b was released that also suggests a hazy atmosphere (Crossfield et al. 2013).

Acknowledgments

This work was partially supported by PRIN INAF 2008 “Environmental effects in the formation and evolution of extrasolar planetary system”. V.N. and G.P. acknowledge partial support by the Università di Padova through the “progetto di Ateneo #CPDA103591”. V.N. acknowledges partial support from INAF-OAPd through the grant “Analysis of HARPS-N data in the framework of GAPS project” (#19/2013). Some tasks of our data analysis have been carried out with the VARTOOLS (Hartman et al. 2008) and Astrometry.net codes (Lang et al. 2010). This research made use of the International Variable Star Index (VSX) database, operated at AAVSO, Cambridge, Massachusetts, USA.

References

All Tables

All Figures

Fig. 1 Left plot: light curve of GJ3470b gathered with the blue channel of LBC in the Uspec band, plotted with the original sampling cadence (upper panel; the black line corresponds to the best-fit model adopted in Table 1) and residuals from the best-fit model (lower panel). Right plot: same as the left plot, but in the F972N20 band of the red channel of LBC. The two light curves are simultaneous.

In the text

Fig. 2 Upper panel: reconstructed transmission spectrum of GJ3470b. The Uspec and F972N20 data points are extracted from our LBC light curves (blue circles), the g′RcIcJ points from Fukui et al. (2013) (green circles), and the 4.5  μm point from Demory et al. (2013) (red circle). The model plotted with a dark green line is scaled from a Howe & Burrows (2012) model computed for a 10-M⊕ cloud-free, haze-free planet with Teq = 700 K and a pure water (H2O) atmosphere. The open squares represent the model points integrated over the instrumental passbands. Second panel: same as above, but with a cloud-free, haze-free, metal-poor atmosphere dominated by H and He, computed by the same authors (orange line). Abundances are scaled from the solar ones (Z = 0.3 × Z⊙). Third panel: same as above, but adding the contribution of a scattering haze from Tholin particles at 1 μbar. Lower panel: instrumental passbands employed for each data point.

In the text

	Fig. 3 Upper panel: same as Fig. 2, but zoomed on the optical region. The red line corresponds to the best weighted linear fit in the (lnλ,k) plane over all data points at λc < 1μm (ours and from Fukui et al. 2013) under the assumption of pure Rayleigh scattering; see text for details. Lower panel: instrumental passbands employed for each data point.
In the text

	Fig. 4 Upper panel: synthetic 1-R⊕ transit (black line) injected in the residuals of our F972N20 light curve (red points and error bars), showing that a short-period Earth around an M1.5V star would have been detectable from our photometry. Lower panel: same as above, binned over 20-min intervals.
In the text