© ESO, 2014

1. Introduction

With the exquisite photometric precision of the space missions, the field of exoplanets has reached several milestones in recent years. To enumerate a few, transiting planets with moderate temperatures like CoRoT-9b (Deeg et al. 2010), with periods of several months, are key objects for studying the effect of the incoming flux from the host stars on their atmospheres. The detection of transit time variations was achieved in the system Kepler-9 (Holman et al. 2010), and they have been used to measure the masses of several of the transiting planets found with the Kepler mission (e.g., Lissauer et al. 2011; Fabrycky et al. 2012; Steffen et al. 2013). All the currently known transiting planets with densities compatible with rocky compositions were detected from space (e.g., CoRoT-7b, Léger et al. 2009; Kepler-10b, Batalha et al. 2011; 55-Cnc e, Winn et al. 2011; Demory et al. 2011), because their shallow transits -of a few hundred ppm- are a challenge for ground-based observations.

While exoplanetary systems have been routinely discovered or confirmed with the radial velocity method, the improvements in photometric precision and transit detection capabilities with the space missions CoRoT (Baglin et al. 2006) and Kepler (Borucki et al. 2010) allow for detection of transit-like features equivalent to objects with the size of an Earth and smaller (e.g., Barclay et al. 2013). The weakest transit signals demand considerable resources, as shown, for instance, with the intense effort needed to measure the mass of CoRoT-7b (Queloz et al. 2009). Not only is the radial-velocity amplitude smaller, but the relative importance of the contribution by stellar activity is also higher, and new techniques need to be used to correct for this. This claim is supported by the more than six different analyses that have been published using the same radial velocity data set. When the transits are found around fainter host stars and/or with longer orbital periods, the current spectrographs do not reach the precision needed to measure the masses of the planets. In these cases, a study of the probabilities of the different scenarios that are compatible with all the follow-up observations performed on the target can result in a validation of a planet candidate. This is the approach started by Torres et al. (2004) and used in recent Kepler discovery papers (e.g., Fressin et al. 2011; Cochran et al. 2011; Ballard et al. 2011; Lissauer et al. 2014; Rowe et al. 2014). For small planets, it is also more difficult to exclude background eclipsing binaries, since fainter background stars could mimic the transit signal.

This paper presents the discovery of the first multiply transiting system found by CoRoT and, as such, is the first one that can be observed from the southern hemisphere observatories. We describe the CoRoT data in the next section. After the analysis, in Sect. 3, of the two transit signals that appear in the data, we describe the follow-up observations in Sects. 4−6. A study of the multiple stellar systems that are compatible with the data is presented in Sect. 7, and the probabilities of these systems are discussed in Sect. 8.

2. CoRoT observations

The target star CoRoT-24 was observed on the detector E1 of the CoRoT satellite during its second Long Run pointing at the direction opposite to the approximate galactic center (“anticenter” in the CoRoT jargon). The observations lasted from 16 Nov. 2008 until 8 March 2009, covering a total span of 111.64 days. The target is faint, having a r mag of 15.1 (see Table 3). A transit signal with a period of 11.76 d was detected in the framework of the alarm mode (Surace et al. 2008) close to the end of the observations, and thus the last 10.7 days of the run are sampled with a 32 s cadence, instead of the nominal 512 s of the rest of the observations. The EN2 data product¹ consists of 45 593 data points, out of which 28 672 are sampled with a 32 s cadence. The aperture used on-board is a monochromatic one, since the telemetry only allows for a recovery of chromatic information on a limited number of targets. Owing to the passage of the satellite over the South Atlantic Anomaly (SAA) and other instrumental effects, about 15% of the data are flagged as bad data points. A detailed description of the performance and pipeline of CoRoT can be found in Auvergne et al. (2009).

As in previous analyses, the power spectrum of the curve without the data points flagged as invalid shows several peaks at the satellite orbital period and its daily aliases. The period of the satellite orbit is about 1.7 h, and thus it might affect the determinations of the transit parameters, because the time scales are comparable. To filter out this signal, we corrected each orbit i with a model that was estimated by computing a Savitzky-Golay smoothed version of the 30 closest orbits to i, normalized independently (Alonso et al. 2008). While this process performs satisfactorily in removing the peak in the power spectrum at the orbital frequency, the daily aliases remain as a consequence of the gaps introduced in the data by removing the points that were flagged as bad. To avoid this side effect, we interpolated the data in the gaps by the use of an inpainting technique (Sato et al. 2010). The power spectrum before and after this correction is plotted in Fig. 1 and the filtered light curve is plotted in Fig. 2.

Fig. 1 Amplitude spectrum of the light curve, showing peaks at the frequencies of the satellite’s orbital signal and the daily harmonics, before and after the correction described in Sect. 2.

Fig. 2 Filtered light curve of CoRoT-24. The black dots are the 512 s sampled data, while the blue curve has been combined in 1 h bins only for display purposes. The inset shows a zoom, where individual transits can barely be detected. The vertical lines show the positions of the transits of CoRoT-24b (dashed light blue) and CoRoT-24c (solid dark blue).

3. Analysis of the CoRoT light curve

A box-fitting least squares (BLS) search (Kovács et al. 2002) shows a peak in the signal detection efficiency (SDE) with an amplitude 13.7 times greater than the standard deviation and no other peaks of comparable amplitude (Fig. 3). Its period of 11.75 d corresponds to the period of the transit identified by the alarm mode, and labelled as CoRoT-24c. This transit was modeled with the Giménez (2006) formalism, as done for previous CoRoT discoveries. As in the case of CoRoT-7b (Léger et al. 2009), the shallow depth of the transits prevents a precise determination of the mean stellar density of CoRoT-24  and we rely on the radius estimation from the spectroscopic analysis described in Sect. 6. The fit parameters were the transit center, the radii ratio k, and the inclination. The limb darkening coefficients² were estimated from Sing (2010) and an error bar of 0.04 was attributed to both u_a and u_b. The amoeba simplex minimization method (Press et al. 1992) was used to find the minimum χ² solution, and a residual permutation (prayer-bead) with a random variation of the limb darkening coefficients according to the values given before was used to estimate the uncertainties of each parameter. The eccentricity was fixed to zero. The phase folded light curve and the best fit model are plotted in Fig. 4. Once the best fit model was found, the transits of CoRoT-24c were removed from the light curve. A new BLS search on the residuals shows a peak at 5.11 days with a shallower depth (second panel of Fig. 3) and a peak with an amplitude 11.8 times larger than the standard deviation. This was fit using the same formalism as above. The phase-folded and best fit curves are plotted in Fig. 5 and the candidate was labelled as CoRoT-24b . As a final step, and to minimize the possible contribution of the transits of CoRoT-24b to the phase-folded curve of CoRoT-24c, we removed the best fit solution of CoRoT-24b from the light curve and repeated the process to estimate the final values for the parameters of CoRoT-24c that are given in Table 4.

Fig. 3 SDE of the BLS search of the original normalized curve, showing the peak due to CoRoT-24c (top); after removal of the transits of CoRoT-24c, showing the peak due to CoRoT-24b (middle), and after removal of b and c, showing no significant peaks (bottom).

Fig. 4 Phase-folded transit of the 11.75 d period CoRoT-24c.

Fig. 5 Phase-folded transit of the 5.11 d period CoRoT-24b.

4. Ground-based photometry

Due to the relatively large point spread function (PSF) of CoRoT that allows one to get color information on the bright targets, there might be several stars inside the PSF that will either cause slightly different parameters obtained from the transit fit if not taken into account, or in the worst case that could have deep eclipses which are diluted by the brightest star in the PSF. Thus, a photometric follow-up with telescopes that provide a narrower PSF is performed in the CoRoT framework. The instruments and techniques are detailed in Deeg et al. (2009), in this section we briefly describe the observations carried out on CoRoT-24.

The BEST-II telescope (Kabath et al. 2009) performed pre-characterization observations of the LRa02 field between November 2007 and February 2008, and none of the stars close to the target detected with BEST-II showed any variation with CoRoT ephemeris. The PSF of these observations has a typical full width at half maximum (FWHM) of about 6′′, thus higher angular resolution than the CoRoT dispersed PSF of 20−30′′. The transit itself is too shallow to be detected with this instrument. The star at about 28′′ SW in Fig. 6 was observed as LRa02-E1-5235 by CoRoT, and it does not exhibit eclipses or transits at the ephemeris of CoRoT-24 b nor c.

We obtained higher resolution images on both the 82 cm IAC-80 and the 2 m Faulkes Telescope North (FTN) telescopes. These observations revealed at least four faint stars at distances of about 10′′ from the target, as seen in Fig. 6, and labelled as A, B, C, and D. Only a small amount of the light from these targets is entering into the CoRoT photometric aperture. If we consider the worst – unrealistic – case in which all the flux from these stars is included into the aperture, they would have to show eclipses of at least 30% in depth to match the depth of the transits of CoRoT-24cand at least 17% to match that of CoRoT-24b. The required depth of the eclipses of each star is summarized in Table 1.

Fig. 6 Region around CoRoT-24, from the image of CoRoT taken at the beginning of the observations (left), imaged with the 1.2 m Euler Telescope (top right), and with the VLT/NACO (bottom right). The red lines delimit the aperture used to extract the photometry.

To check for this worst case scenario and in an effort to try to refine the ephemeris of CoRoT-24c, we observed the target in two sequences using the 1.2 m Euler Telescope, on 2 January 2010, from 00:41 to 09:37 UT, and with the new EulerCam detector on 12 November 2010, from 5:24 to 08:15 UT. None of the four stars showed variations larger than 20% for the first sequence, and 10% for the second one. While the second sequence only covers from slightly after the middle of the predicted transit to the 1σ upper limit of the fourth contact (using the ephemeris from Table 4), the January 2010 observations should have detected an ingress or egress of any of the four stars considering the error in the ephemeris. Since this was not the case and, furthermore, since none of the four stars are fully included into the CoRoT aperture, we conclude that the transits of CoRoT-24c take place on the target. Three new stars that were closer to the target were detected in the Euler images, but they are too faint to explain the transits of CoRoT-24 b or c.

We checked for eclipses on the four stars during one expected transit of CoRoT-24b on 29 December 2011, centered at 06:44 UT using the ephemeris of Table 4. We took images in four different 20-min observing windows centered at about 02:15, 04:15, 05:15 and 06:45 UT, also with the EulerCam detector, and built difference images with a similar observing sequence taken a few days after. None of the images reveals a change in flux of any of the four stars of more than 18%. Since the 1σ error in the ephemeris at this epoch was 2h50m and the transit duration about 2h51m, these observations discard background eclipsing binaries with the ephemeris of CoRoT-24b only partially, covering eclipses centered from about −2σ to +0.3σ from the ephemeris. Further observations are needed to extend this coverage to eclipses happening later than expected.

Finally, we explored the possibility of having other stars closer to the target that lie in the glare of CoRoT-24 and that would have not been detected in the other photometric observations. We used the instrument NACO at the VLT to obtain a deep image in J filter with high angular resolution (Guenther et al. 2013). This filter was selected in order to detect faint background stars that have red colors due to the extinction.The combined image is shown in Fig. 6, where the faint stars that were distinguishable in the Euler images now appear well resolved. A normalized PSF profile of CoRoT-24 is shown in Fig. 7, together with the flux limits that a star should reach in order to reproduce the transits of CoRoT-24 b and c, for the conservative cases of a 100% eclipse and a 50% one. To remain undetected in this high resolution image, potential eclipsing binaries should lie closer than 1.3′′ and 0.9′′ to the target to explain the targets of CoRoT-24 b and c, respectively³.

Fig. 7 Normalised PSF profile around CoRoT-24, from a VLT/NACO combined image. The horizontal lines delimit the flux level that a contaminant included in the CoRoT ’s PSF should reach in order to produce the observed diluted transit depth. Two lines are plotted, for an unrealistic contaminant completely disappearing during eclipse (100% eclipse line), and a 50% eclipse. Dashed line is for CoRoT-24c, while dotted line is for CoRoT-24b. Realistic potential background eclipsing binaries should lie closer than 1.3′′ and 0.9′′ to the main star to explain the transit depths of CoRoT-24b and c, respectively.

5. Radial velocities

Two of the spectrographs that have shown to reach the precision required to discover super-Earth like planets were used in an effort to measure the masses of the planets in the system (see Santerne et al. 2011 for details on the CoRoT radial velocity follow up strategy). 37 spectra were taken with the Keck1/HIRES spectrograph (Vogt et al. 1994) as part of NASA’s key science project in support of the CoRoT mission, and 34 measurements were performed with the HARPS/ESO-3.6 spectrograph (Mayor et al. 2003). These measurements are summarized in Table 2, and plotted in Fig. 8. The data show a long term modulation that is compatible with the presence of a longer period object in the system, with a period of ~940 d, a K ~ 35 m s^-1 and an eccentricity close to 0.4. It would correspond to an object with a minimum mass of ~1.5 M_Jup and a semi-major axis of 1.8 A.U. The obtained parameters for this candidate planet rely strongly on the two single data points taken in 2010, which sample the phases of minimum observed radial velocity. Further observations are thus encouraged to confirm and better constrain the orbital parameters of this candidate signal, especially around the phases of the next predicted maximum radial velocity, which should happen around January, 2015.

The periodogram of the combined HARPS + HIRES data set, after the removal of the long period signal, is shown in Fig. 9. Several peaks are identified at around 27.3 d, 13.5 d, 9 d and 6.6 d, that could correspond to the rotational period of the star and its harmonics (P_rot/2, P_rot/3, P_rot/4). A rotational period of 27.3 d is compatible with the value of the v sin i_⋆ reported below assuming i_⋆ close to 90° and with the autocorrelation of the photometric curve.

Fig. 8 Radial velocities, showing a long-term trend compatible with a 1.5 MJup object with a period of 940 d. Blue points for HIRES data, red points for HARPS. More data are needed, especially at the phases of maximum RVs.

Fig. 9 Periodogram of the combined HARPS + HIRES data set, after the removal of the ~940 d signal. A low-significance excess of power around the period of CoRoT-24c (11.759 d) is detected, and zoomed in in the inset plot. The vertical lines show the positions of the photometrically determined period, and the dashed and dotted lines the locations of the harmonics due to the observational window.

The region close to the photometrically determined period of CoRoT-24c shows increased power, taking the observational RV window function into account. In fact, the peak with most power is consistent with the second harmonics of the window function, and it is found at ~11.038 d. Since the orbital period is well determined from the photometry, we evaluated the false alarm probability of this signal by constructing sets of mock data shuffling the RV data points, keeping the same time stamps. We computed the integrated power of the periodogram in the regions between the period of CoRoT-24c and three times the error listed in Table 4 and in each of the harmonics due to the window function (which shows two peaks in frequency at 0.00277 and 0.005535 c/d, corresponding to the annual and seasonal observing windows). The false alarm probability was estimated by counting the number of times the integrated power in these regions was higher than in the data, resulting in a value of 11%.

Under the assumption that the RV signal from the two transiting candidates is present in the data, we first estimated the respective K-amplitudes of CoRoT-24b and CoRoT-24c by performing a two-sine component fit to the data after removing the long-term trend. The periods and phases were fixed to the transit ephemerides. The orbital solution was calculated using the general nonlinear least squares fitting program GaussFit (Jefferys et al. 1988). This resulted in K_b = 0.80 ± 3.0 m s^-1 and K_c = 5.2 ± 3.2 m s^-1. This suggests that any short-term variations in the RV measurements are dominated by CoRoT-24c.

We then fitted the RV measurements using a variation of the floating offset fitting technique employed in detecting CoRoT-7b (Hatzes et al. 2011). The RV data were divided into eleven subset data, four for Keck and seven for HARPS. The time span of each data subset ranged from two to seven days, with a median value of five days. On several nights there were three Keck RV measurements and for the fit we used the median value. An orbital fit was made keeping the period and phase fixed to transit values, but allowing the K-amplitude and zero-point offsets of the subset data to vary.

The basis for this method is that during the short time interval spanning each subset the variations due to the long-term trend or stellar variations from rotational modulation is relatively constant. Although the time span of the subset data is a significant fraction of the rotational period (0.07 to 0.17) the star has a relatively low level of activity, so our assumption seems reasonable. This technique has the advantage of removing the long-term trend without a functional fit. It also provides a natural way of combining the HARPS and Keck data, each of which have their own velocity zero-point.

The floating offset fit to the data using the period of CoRoT-24c resulted in K_c = 8.31 ± 2.86 m s^-1. Applying the method using the period of CoRoT-24b resulted in K_b = 0.88 ± 2.46 m s^-1. We then tried fitting the residual velocities using the period of CoRoT-24b, but after first subtracting the contribution of the derived orbit of CoRoT-24c. This resulted in K_b = 0.0 ± 2.1 m s^-1. The floating offset method confirmed the results found by the two-component sine fitting, namely that CoRoT-24c is marginally detected (2.9σ) while CoRoT-24b is undetected.

Fits using the individual RV data sets gave consistent results. A fit using only the Keck data resulted in K_c = 6.64 ± 3.60 m s^-1. The fit to only the HARPS data resulted in K_c = 10.76 ± 4.62 m s^-1.

Figure 10 shows the RV variations phased to the transit ephemeris of CoRoT-24c after removing the calculate zero-point offsets from the respective subset data. To better see the RV variations we also showed binned values.

The obtained semi-amplitudes correspond to a 1σ upper limit for the mass of CoRoT-24b of 5.7 M_⊕ and a mass of 28$\begin{array}{c}\mathrm{+}\mathrm{11}\\ -11\end{array}$M_⊕ for CoRoT-24c  taking the host star mass estimation described in the next section into account.

Fig. 10 Radial velocity measurements folded at the period of CoRoT-24c, and best fit orbital solution with the parameters given in the text. Blue dots are for HIRES, red triangles for HARPS, and black squares data binned for displaying purposes.

6. Spectral analysis

Seven HARPS spectra of CoRoT-24that were not contaminated by the moon light were co-added once set in the rest frame of the star. This provides us a spectrum with a signal-to-noise ratio of 88 per element of resolution. We determined the photospheric parameters of the host star using the same methodology as in past CoRoT papers: the spectroscopic analysis was done independently using the semi-automated package vwa (Bruntt et al. 2010, and references therein) and sme (Valenti & Piskunov 1996). The parameters obtained are listed in Table 4. The mass and radius of the star were then estimated using the T_eff, metallicity and logg derived from the spectroscopic analysis, to fit the evolutionary tracks of starevol (A. Palacios, priv. comm.; see Turck-Chièze et al. 2010; Lagarde et al. 2012 for a description of the code). It yields a mass of ${\mathit{M}}_{\mathrm{\ast }}\mathrm{=}\mathrm{0.9}{\mathrm{1}}_{-0.07}^{\mathrm{+}\mathrm{0.08}}$ M_⊙  and ${\mathit{R}}_{\mathrm{\ast }}\mathrm{=}\mathrm{0.8}{\mathrm{6}}_{-0.04}^{\mathrm{+}\mathrm{0.03}}$ R_⊙  for the star. The errors are formal statistic errors, which are known to provide unrealistically low values due to, among other things, the intrinsic errors of the stellar models. We thus adopted in Table 4 error values of 10% for both mass and radius of the star, to account for these extra errors.

We used the derived photospheric parameters and the available photometry to fit and scale a model spectrum, in order to determine the interstellar extinction and the distance to the system. The spectral energy distribution (SED) nicely follows the scaled model, as shown in Fig. B.3 of Guenther et al. (2013), and shows no evidence of any mid-infrared excess.

The projected stellar rotation velocity was determined by fitting the profile of several well isolated spectral metal lines. We found v sin i_⋆ = 2.0$\begin{array}{c}\mathrm{+}\mathrm{1.0}\\ -1.5\end{array}$ km s^-1. The autocorrelation of the light curve shows a peak at about 29 days (Fig. 11), which is consistent with a rotation with the v sin i_⋆  found above. As in previous cases, we found models that reproduce the spectroscopic parameters corresponding to a pre-main sequence evolutionary phase. Due to the v sin i_⋆  and the estimated rotation period, we did not consider these solutions, because the star should be rotating at higher velocity. The age of the star in the main sequence solutions is around 11 Gyr.

Fig. 11 Autocorrelation of the normalized light curve.

In order to search for late-type stars close to CoRoT-24, whose brightness difference is much smaller in the infrared than in the optical, we obtained a high resolution infrared spectrum using the CRIRES spectrograph which is mounted on the VLT-UT1 (Antu), as part of a program for follow up of CoRoT candidates. The strategy, instrumental configuration and results for several candidates, including CoRoT-24, is described in Guenther et al. (2013, the case of CoRoT-24 is described in their Appendix B.6). These observations exclude non-resolved physically bound companions with spectral types earlier than M2.5V.

7. Blend simulations

Without a very significant detection of the radial velocities, we need to carefully check different blend scenarios that might explain all the observations described in the previous sections. In the next paragraphs, we describe the technique we used to simulate and evaluate the probabilities of the transit signals to be due to different types of stellar systems involving eclipsing binaries or stars with bigger transiting planets. The method is quite similar and strongly inspired by the blender software described in Torres et al. (2011) and Fressin et al. (2011). We apply it to systems with either four or five bodies, because we consider the two periodic transit signals real, and there is no stellar system that can remain stable with such short periods. The simulations of systems with four bodies account for those configurations in which there is one transiting planet around one star, and an eclipsing binary that can be physically associated or in the line of sight. Simulations with five bodies are for systems where the star 1 (the one detected in the spectra) is blended with two stars, each in an eclipsing system or containing a transiting planet, also at a wide variety of distances.

As a starting point, we used the cross-correlation function (CCF) of the HARPS spectra to put flux limits to other stars that are included in the spectrograph’s fiber. The main peak of the CCF was subtracted from the individual exposures, and fake CCFs with different widths were inserted at different velocities. We quantify a 3σ level for non-detection of CCFs corresponding to flux ratios of 0.01 of the main star. This condition was included in the simulations described below.

We constructed light curves with the contribution of one of the planets removed according to the results of the fits described in the previous section⁴. For each of these curves, we built a phase folded light curve with the period of the remaining transit. The goal of having a whole orbital phase curve for each of the transiting objects is to check for potential secondary eclipses at any of the orbital phases, since eccentric false positives will also be checked. To construct these curves, we normalized each of the planetary orbits with a 7th order polynomial, thus removing most of the variations due to the stellar activity or uncorrected instrumental effects. We checked that these low order polynomials would not be able to remove the short (compared to the orbit) secondary eclipses. Finally we combined all the orbits into a normalized phase curve binned in 0.0012 in phase (total of 834 points). To check for systems whose period is twice the one of the potential planet, which correspond to stellar systems having similar primary and secondary eclipses, we also built curves folded at twice the period. In total, four light curves were built this way, two for each of the transiting signals to be studied.

These light curves will be compared to many different models of four or five body systems, as described before. In the five body scheme, we will call star 1 the one whose parameters were extracted in the spectroscopic analysis of Sect. 6, stars 2 and 3 are the components of one of the eclipsing systems, and stars 4 and 5 the components of the other. In the four body scheme, the only eclipsing binary is composed of stars labelled 3 and 4.

For each of the models, the stellar parameters of star 1 are fixed to the results of the spectral analysis. We interpolate the ATLAS9 model atmospheres of Castelli & Kurucz (2004) to the values of the T_eff, log g and metallicity, scale it with the stellar radius and its distance, apply an extinction law of Fitzpatrick (1999) to the E(B − V) of the galactic interstellar extinction model of Amôres & Lépine (2005), and finally apply the spectral response function of CoRoT given in Auvergne et al. (2009) to get the observed flux of star 1.

7.1. Four-body simulations

The star 1 is orbited by one planet (2), and we simulated eclipsing binaries (systems of star 3 and 4, where 4 is also allowed to be a planet as a particular case). We explored systems where the period of the planet 2 is that of 5.11 days, and systems where it was 11.76 d. In each case, the system 3-4 will be assigned the other period, or twice this value. This is done to account for the possibility of star 3 and 4 being similar, and thus having no significant differences between primary and secondary eclipses. We extracted one value for the age and mass of the star 3-star 4 eclipsing system in the model from an uniform random distribution. We also chose a value for the eccentricity and the longitude of the periastron randomly from an uniform distribution, and a random distance for the eclipsing system. This distance is allowed to vary from 0 to 10 kpc for the observations towards the anticenter direction. Once these random values have been generated for stars 3 and 4, we interpolate the isochrones of Marigo et al. (2008)Girardi et al. (2010)⁵, assuming solar metallicity, to get the values of the T_eff and log g. Then the ATLAS9 model atmospheres are interpolated to these values and scaled with the stellar radius and distance coming from the isochrones. We apply the same extinction law, model and CoRoT response function to these spectra, and obtain the observed fluxes of stars 3 and 4. The semi-major axis is calculated from the masses and the period of stars 3 and 4 and the orbital inclination is estimated from the primary eclipse duration with the associated error bars (measured in the light curve), using Eq. (7) in Tingley & Sackett (2005). We can then compute the model light curve of the star 1-planet 2 and of the system star 3-star4, using the jktebop code (Southworth et al. 2004a,b, based on ebop Popper & Etzel 1981; Etzel 1981) and the limb-darkening coefficients of Sing (2010), dilute the system using the observed flux of star 1, and calculate the χ² with respect to the CoRoT-24 phase light curves. This process is repeated for many (6 × 10⁷) different random realizations of the parameters of the star 3 and 4, and of planet 2. In order to speed-up the process, several conservative conditions are plugged in the system in order not to compute unnecessarily the eclipsing binary models: for instance, if there are no eclipses that will match the observed duration, or if the undiluted eclipse depth will be in any case smaller than the observed transit depth. We also impose the condition to have the star 1 brighter than the total of star 3 + star 4. In this case, our models have 13 free parameters, five for the planet around star 1 (mass, radius, inclination, eccentricity and omega), and eight for the star 3-star 4 system (age, initial masses of 3 and 4, distance, inclination, eccentricity, omega, and possibility of having 2× period). In the case of the body 4 being a planet the number of free parameters remains unchanged (age, initial mass of 3, radius and mass of the planet, distance, inclination, eccentricity, omega).

7.2. Five-body simulations

We simulated two eclipsing systems, namely star 2-star 3, and star 4-star 5. In each of the cases the period is allowed to be one of the periodic signals found in the light curve, or twice this value. Solutions where the system 3 or 5 is a planet were also included. The total number of free parameters is 16, i.e., 8 for each of the systems, the same as in the case of the star 3-star 4 scenario of the four body simulations.

7.3. F-tests

To compare the blend and multiple transiting planetary system scenarios we used the F-test, that is based on the ratio of two reduced χ² . The χ² maps were built in the parameter space mass of star 2 (or 3, or 4, depending on the case) vs distance modulus ($\mathrm{\Delta }\mathit{\delta }\mathrm{=}\mathrm{5}{\log }_{\mathrm{10}}{}^{\mathrm{\right(}}\frac{{\mathit{D}}_{\mathrm{blend}}}{{\mathit{D}}_{\mathrm{star}\mathrm{1}}}^{\mathrm{\left)}}$) for each scenario. In each box of the two-dimensional grid, we took the minimum χ² of all the models inside that box, and we computed $\mathit{F}\mathrm{=}{\mathit{\chi }}_{\mathrm{reduced}\mathit{,}\mathrm{blend}}^{\mathrm{2}}\mathit{/}{\mathit{\chi }}_{\mathrm{reduced}\mathit{,}\mathrm{planetarysystem}}^{\mathrm{2}}$.

The results for all the different configurations explored are plotted in Fig. 12. Since the radial velocity results provide a clearer detection of the CoRoT-24c signal, we describe in more detail the results of the simulations in which one of the stars has a transiting planet with the period of CoRoT-24c and explore the possible scenarios that are compatible with the data. There are typically two regions in the maps that are not excluded above a 3σ level. One corresponds to eclipsing binaries with the primaries being hot and massive stars, with masses higher than 1.7 M_⊙. In these cases, a CCF might have not been detected in the radial velocity observations, but they need to lie at large distances in order not to become the brightest member of the system, what would give results that are inconsistent with the spectral analysis. The other set of solutions is for eclipsing binaries with masses of the primary below 0.9 M_⊙, at a distance in which the diluted eclipses could be mistaken with the transits of a planet, but not so close as to have been detected in the CCF of the spectra.

The five body results are more restrictive, as they require none of the primaries of the two blended eclipsing systems to be detected in the CCFs, and each of the systems needs to lie at a particular physical distance with respect to the brightest object (star 1).

The first row of Fig. 12 shows the case of two planets orbiting around star 2, that after dilution created by star 1 would appear as small planets. Most of the allowed solutions are for sizes of the planets of >1 R_Jup. Such a system with two Jupiter-sized planets orbiting that close to their star has not been discovered yet, and multiple transiting systems are known to rarely harbour a transiting giant planet (Latham et al. 2011).

8. Discussion

8.1. Probabilities for the different scenarios

We used an average image from the second Euler transit sequence in order to estimate the density of stars in the direction of the CoRoT-24 observations. The FWHM of the PSF in this image is about 1′′. Figure 13 shows the cumulative number of stars per square arcsec, and the upper magnitude limits that a star can have in order to reproduce the observed depths, neglecting color effects that might move these lines slightly if the colors of the contaminants are different between the CoRoT and the RG bandpasses. From this figure, we can (conservatively) estimate the odds for a random alignment of a background/foreground star, that would remain undetected in the NACO high resolution images, and whose brightness is enough to explain one of the two transiting signals, to about 1:180 (=1:1/0.0055, corresponding to the crossing of the vertical line on the right of the plot and the solid line). To this, we need to add the fact that the non resolved star should be an eclipsing binary. The fraction of eclipsing binaries found by CoRoT in previous fields (IRa01: 145 out of 9872 sources, 1.47% Carpano et al. 2009, LRc01: 158 out of 11 408 sources, 1.38% Cabrera et al. 2009) is similar in the center and anticenter directions, and comparable to the results of the Kepler mission (2165 out of 156 453 sources, 1.38% Slawson et al. 2011). Out of these, close to 50% are not “appropriate” to be confused with a transiting planet, as they belong to the contact or semi-detached systems. Thus, the fraction is on the order of 0.7−0.8%, and the odds to have a star that is an eclipsing binary non resolved in the NACO image drop to 1:22 000−1:26 000.

Fig. 13 Cumulative number of stars per square arcsec, computed from a combined image of the field of CoRoT-24 that is complete down to RG ~25. The two vertical lines correspond to the extremes of the horizontal lines in Fig. 7, i.e., the highest magnitude a star could have in order to reproduce the observed transit depth of CoRoT-24b while completely disappearing during eclipses, and the magnitude of a star with 50% eclipses needed to explain the transit depth of CoRoT-24c.

According to the results of the previous section, the eclipsing binary has also to lie at a particular physical distance from the main target, and have an adequate primary mass and q-ratio in order to reproduce the photometric observations. The constraints imposed by the CCF and the photometry exclude most⁶ binaries with a primary with a mass from 0.9 to 1.7 M_⊙. Slightly different distances or q-ratios could produce blended transits with depths that would be confused with that of a planet of different size as that of CoRoT-24b. Using the distribution of detached eclipsing binaries of Slawson et al. (2011) we consider that we can safely exclude up to 50% of the potential eclipsing binaries because their primaries are in the mass range mostly excluded by the CCF. We estimate the odds then to have a star that is an adequate eclipsing binary, non resolved in the NACO image, and not detected in the CCF of 1:44 000−1:52 000.

Because the radial velocities of the brightest star do not show variations at the typical levels for stellar companions (km s^-1) and since stellar eclipsing triple systems with periods of five and eleven days are dynamically not stable, we have to estimate the probability to have two adequate eclipsing binaries inside the CoRoT PSF. With the photometric follow-up we have excluded the case where the CoRoT-24c eclipses happen in any of the stars that are at about 10′′ in Fig. 6. But, due to the non complete verification for eclipses with the ephemeris of CoRoT-24b (excluded in a window from −2σ to +0.3σ of the predicted ephemeris), we can currently only estimate the probability that one of the five potentially worrying stars is an adequate eclipsing binary as ~5:130 (5 stars – the 4 labelled in Fig. 6, plus a hypothetical star non resolved in the NACO images – and a 0.7−0.8% fraction of appropriate eclipsing systems). Thus, the probability to have a blend of two eclipsing binaries and a single target star can be currently estimated as <8.8 × 10^-7 (=1:44 000 × 5:130). If future photometric campaigns are able to definitely reject the possibility of having eclipses in the detected stars at projected distances of about 10′′ from the target (see Sect. 4) then this probability could drop to <5.2 × 10^-10 (=1:44 000 × 1: 44000).

This is to be compared with the probability of having a genuine Neptune transiting system. Using the results from the Kepler mission (Borucki et al. 2011)⁷, we count 162 planetary candidates with sizes between 3.3 and 5.5 times the size of the Earth, out of the ~190 751 stars surveyed by the mission. Assuming that all these are true planets, we get that a fraction of 0.85 × 10^-3 of the stars host a planet with a size in the range of the candidates around CoRoT-24. Alternatively, we can also estimate the fraction of multiple transiting systems with two transiting components using the results of Lissauer et al. (2014), as 272 candidates out of the observed stars with Kepler, representing a fraction of 1.4 × 10^-3. This later number includes also systems with sizes smaller than those detectable in the CoRoT-24 light curve, and thus constitutes an upper limit.

Thus, we estimate the five body scenario as >966 times less probable than having at least one transiting planet around the main target (=0.85 × 10^-3/8.8 × 10^-7).

Out of the four body scenarios, in which there is one genuine transiting planet and an eclipsing binary inside the CoRoT PSF, the blend scenario having the highest probability is CoRoT-24c being a genuine transiting planet, while one of the five stars at about 10′′ is an adequate eclipsing binary. Using the same assumptions as above, we estimate the probability to have such a configuration as 3.3 × 10^-5 (=5:130 × 0.85 × 10^-3). If the close targets are excluded as eclipsing binaries (either with a photometric campaign of with a few radial velocities taken in each of the stars), then the eclipsing binary should be non resolved in the NACO image, with a probability of <1:44 000 = 2.3 × 10^-5 estimated above.

Taken together, the scenario in which one of the transiting signals is due to a planet is at least 966 times more plausible than a double-blend system. The scenario of a multiple transiting system in which both transit signals are due to planets is at least 26× more probable than a transiting planet + eclipsing binary blend (=0.85 × 10^-3/3.3 × 10^-5). This can rise to 37× if the five stars at a distance of about 10′′ from the target are definitely excluded as eclipsing binaries.

8.2. Dynamical stability

We calculated the orbital separation of CoRoT-24b and CoRoT-24c in Hill-units, to get an idea of the stability of the system. With the upper limit estimated for CoRoT-24b and the measurement of the mass of CoRoT-24c  we estimate the orbital separations in Hill units of around eight. A two-planet system is supposed to be stable in Gyr time scales if this number is larger than seven (Funk et al. 2010), so we conclude that dynamical stability is not a problem for the CoRoT-24 system.

9. Conclusion

The light curve of CoRoT-24 is compatible with the transits of two Neptune-sized planets (with radii of 3.7 and 5.0 R_⊕) at 5.11 and 11.76 d periods. We performed different photometric and spectroscopic measurements to try to secure the planetary nature of the objects and measure their masses and densities. The relative faintness of the target and the expected low amplitude of the radial velocity variations do not allow a precise mass measurement of CoRoT-24b and we provide an upper mass limit of 5.7 M_⊕, while for CoRoT-24c we obtain a tentative mass of 28 ± 11 M_⊕. The deduced low densities, of <0.9 g cm^-3 and 1.3$\begin{array}{c}\mathrm{+}\mathrm{0.5}\\ -0.4\end{array}$ g cm^-3 are indicative planets with significant gaseous envelopes.

We have studied the different configurations of four and five body systems that might reproduce the photometric light curve, which allowed a rejection of several configurations. A combination of these simulations with a study of the stellar number density in the direction of the target shows that the planetary scenario is at least 900 times more probable than a blend involving two eclipsing binary systems, and the multiple transiting planets scenario is at least 26 times more probable than a single planet + blend involving one eclipsing binary. In our reasoning, we have not considered other factors that might contribute to the multiple planet scenario, as the fact that there is a tentative detection of a ~940 d, slightly eccentric ~1.5 M_Jup object orbiting around the main target. As there is not enough statistics in the current sample of known multiplanetary systems, we have not used this information to favor or penalize the planet scenario.

Thus CoRoT-24 is most probably a multiple transiting system consisting of two Neptune-sized planets and a third more massive object at a larger orbital separation. They orbit a K1V star with a higher metallicity than the Sun, that shows photometric activity at a level of about 2% with time scales of about 29 d. We interpret this as an indication of the stellar rotation period. The system is thus a good example of the Neptune-sized multiple transiting candidates found in the Kepler mission Borucki et al. (2011), and shows the intense ground-based effort needed to confirm these objects and measure their masses, due in part to the relative faintness, but also to the additional difficulty of having non-transiting or longer period planets in the system that complicate the analysis.

To proceed further, more precise light curves of the transiting objects would allow a narrower range of scenarios that might reproduce all the observables. A photometric search for deep eclipses in the companions at about 10′′ of the target should be performed with the ephemeris of the more challenging candidate CoRoT-24b. Alternatively, these companions could be detected spectroscopically. Despite being challenging for ground-based observations, the progress in the recent years have shown that it might be possible to detect the transits of CoRoT-24c (and probably also CoRoT-24b) and thus refine the ephemeris.

A more intense radial velocity campaign devoted to CoRoT-24 might result in a more precise measurement of the masses of the planets, and thus their densities. These data could serve to better model the long period signal seen in the radial velocities, and to evaluate the effects of additional planetary companions in the system not considered in this work. Since the star shows photometric activity at the 2% level and a period of around 29 days, the observing strategy should be designed in order to minimize or to understand better the effect of the activity on the measured masses. The star is currently the only multiple transiting system observable from several observatories at the southern hemisphere. Thus, CoRoT-24 could be a good target for future high resolution spectrographs planned for bigger aperture telescopes, such as ESPRESSO/VLT.

Online material

Fig. 12 Maps showing the regions where the F-test resulted in solutions not excluded to a level above 3σ with the multiple transiting system solution (solutions excluded above 3σ appear white in the plot). The different configurations explored are sketched at the left of each plot. The color scale is for the different q-ratios of the binary star or for the radii of the planets, and the color intensities are for different levels of exclusion (more intense meaning less excluded). The ordinate axis shows the distance modulus ($\mathrm{\Delta }\mathit{\delta }\mathrm{=}\mathrm{5}{\log }_{\mathrm{10}}{}^{\mathrm{\right(}}\frac{{\mathit{D}}_{\mathrm{blend}}}{{\mathit{D}}_{\mathrm{star}\mathrm{1}}}^{\mathrm{\left)}}$). Bound, or closer to the observer, second stars are excluded with CRIRES spectra for spectral types later than M2.5V (Guenther et al. 2013). Regions labelled as “partially excluded by HARPS" are those in which a second CCF would have been detected in the HARPS data. In the case of rapidly rotating stars, this second CCF might have been unnoticed. In every case, second stars with masses above 1.7 M⊙ have been assumed to be undetectable in the CCFs, which explains the vertical cuts at this value. The upper left parts of the plots are typically excluded because the hypothetical system would produce eclipse depths below the observed values. The diagonal solid lines labelled with different J values are plotted to provide a direct comparison with the exclusion values from the NACO J-band images given by Guenther et al. (2013): a J = 17 star would have been detected at (and above) 0.18′′ from the primary, a J = 18 at 0.23′′, and a J = 19 at 0.58′′.

Acknowledgments

We wish to thank the ESO staff for their support and for their contribution to the success of the HARPS project and operation at ESO La Silla Observatory. Some of the data presented are from observations made with the IAC80 operated at Teide Observatory on the island of Tenerife by the Instituto de Astrofísica de Canarias. The team at the IAC acknowledges support by grants ESP2007-65480-C02-02 and AYA2010-20982-C02-02 of the Spanish Ministry of Science and Innovation (MICINN). The CoRoT/Exoplanet catalogue (Exodat) was made possible by observations collected for years at the Isaac Newton Telescope (INT), operated on the island of La Palma by the Isaac Newton group in the Spanish Observatorio del Roque de Los Muchachos of the Instituto de Astrofísica de Canarias. The German CoRoT team (TLS and University of Cologne) acknowledges DLR grants 50OW0204, 50OW0603, and 50QM1004. The French team wishes to thank the CNES and the French National Research Agency (ANR-08- JCJC- 0102-01) for their continuous support to our planet search program. The Swiss team acknowledges the ESA PRODEX program and the Swiss National Science Foundation for their continuous support on CoRoT ground follow-up. R. Alonso acknowledges support from Spanish Ministry through its Ramón y Cajal program RYC-2010-06519. A.S. Bonomo acknowledges a CNES grant. S. Aigrain acknowledges STFC grant ST/G002266. M. Gillon acknowledges support from the Belgian Science Policy Office in the form of a Return Grant. M. Endl, W. D. Cochran and P. J. MacQueen were supported by NASA Origins of Solar Systems grant NNX09AB30G. T. Mazeh acknowledges support from the Israeli Science Foundation (grant No. 655/07). A. Santerne acknowledges support from the European Research Council/European Community under the FP7 through Starting Grant agreement number 239953. This research has made use of the SIMBAD database, operated at the CDS, Strasbourg, France, and of NASA’s Astrophysics Data System. Part of the data presented herein were obtained at the W. M. Keck Observatory from telescope time allocated to the National Aeronautics and Space Administration through the agency’s scientific partnership with the California Institute of Technology and the University of California. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation.

References

All Tables

All Figures

	Fig. 1 Amplitude spectrum of the light curve, showing peaks at the frequencies of the satellite’s orbital signal and the daily harmonics, before and after the correction described in Sect. 2.
In the text

	Fig. 2 Filtered light curve of CoRoT-24. The black dots are the 512 s sampled data, while the blue curve has been combined in 1 h bins only for display purposes. The inset shows a zoom, where individual transits can barely be detected. The vertical lines show the positions of the transits of CoRoT-24b (dashed light blue) and CoRoT-24c (solid dark blue).
In the text

	Fig. 3 SDE of the BLS search of the original normalized curve, showing the peak due to CoRoT-24c (top); after removal of the transits of CoRoT-24c, showing the peak due to CoRoT-24b (middle), and after removal of b and c, showing no significant peaks (bottom).
In the text

	Fig. 4 Phase-folded transit of the 11.75 d period CoRoT-24c.
In the text

	Fig. 5 Phase-folded transit of the 5.11 d period CoRoT-24b.
In the text

	Fig. 6 Region around CoRoT-24, from the image of CoRoT taken at the beginning of the observations (left), imaged with the 1.2 m Euler Telescope (top right), and with the VLT/NACO (bottom right). The red lines delimit the aperture used to extract the photometry.
In the text

Fig. 7 Normalised PSF profile around CoRoT-24, from a VLT/NACO combined image. The horizontal lines delimit the flux level that a contaminant included in the CoRoT ’s PSF should reach in order to produce the observed diluted transit depth. Two lines are plotted, for an unrealistic contaminant completely disappearing during eclipse (100% eclipse line), and a 50% eclipse. Dashed line is for CoRoT-24c, while dotted line is for CoRoT-24b. Realistic potential background eclipsing binaries should lie closer than 1.3′′ and 0.9′′ to the main star to explain the transit depths of CoRoT-24b and c, respectively.

In the text

	Fig. 8 Radial velocities, showing a long-term trend compatible with a 1.5 MJup object with a period of 940 d. Blue points for HIRES data, red points for HARPS. More data are needed, especially at the phases of maximum RVs.
In the text

	Fig. 9 Periodogram of the combined HARPS + HIRES data set, after the removal of the ~940 d signal. A low-significance excess of power around the period of CoRoT-24c (11.759 d) is detected, and zoomed in in the inset plot. The vertical lines show the positions of the photometrically determined period, and the dashed and dotted lines the locations of the harmonics due to the observational window.
In the text

	Fig. 10 Radial velocity measurements folded at the period of CoRoT-24c, and best fit orbital solution with the parameters given in the text. Blue dots are for HIRES, red triangles for HARPS, and black squares data binned for displaying purposes.
In the text

	Fig. 11 Autocorrelation of the normalized light curve.
In the text

Fig. 13 Cumulative number of stars per square arcsec, computed from a combined image of the field of CoRoT-24 that is complete down to RG ~25. The two vertical lines correspond to the extremes of the horizontal lines in Fig. 7, i.e., the highest magnitude a star could have in order to reproduce the observed transit depth of CoRoT-24b while completely disappearing during eclipses, and the magnitude of a star with 50% eclipses needed to explain the transit depth of CoRoT-24c.

In the text

Fig. 12 Maps showing the regions where the F-test resulted in solutions not excluded to a level above 3σ with the multiple transiting system solution (solutions excluded above 3σ appear white in the plot). The different configurations explored are sketched at the left of each plot. The color scale is for the different q-ratios of the binary star or for the radii of the planets, and the color intensities are for different levels of exclusion (more intense meaning less excluded). The ordinate axis shows the distance modulus ($\mathrm{\Delta }\mathit{\delta }\mathrm{=}\mathrm{5}{\log }_{\mathrm{10}}{}^{\mathrm{\right(}}\frac{{\mathit{D}}_{\mathrm{blend}}}{{\mathit{D}}_{\mathrm{star}\mathrm{1}}}^{\mathrm{\left)}}$). Bound, or closer to the observer, second stars are excluded with CRIRES spectra for spectral types later than M2.5V (Guenther et al. 2013). Regions labelled as “partially excluded by HARPS" are those in which a second CCF would have been detected in the HARPS data. In the case of rapidly rotating stars, this second CCF might have been unnoticed. In every case, second stars with masses above 1.7 M⊙ have been assumed to be undetectable in the CCFs, which explains the vertical cuts at this value. The upper left parts of the plots are typically excluded because the hypothetical system would produce eclipse depths below the observed values. The diagonal solid lines labelled with different J values are plotted to provide a direct comparison with the exclusion values from the NACO J-band images given by Guenther et al. (2013): a J = 17 star would have been detected at (and above) 0.18′′ from the primary, a J = 18 at 0.23′′, and a J = 19 at 0.58′′.

In the text