© ESO, 2017

1. Introduction

More than 3000 extra-solar planets have been detected to date, mostly with radial velocity (RV) and transit techniques¹. Beside these prolific indirect detection methods, about 30 young, wide-orbit, giant exoplanets or brown dwarfs have been discovered using direct imaging techniques (DI). Analyses of the DI or RV surveys have brought in particular new constraints on exoplanet formation and highlight correlations between giant planet formation with stellar mass or stellar metallicity (Biller et al. 2007; Lafrenière et al. 2007; Johnson 2008; Johnson et al. 2010; Chauvin et al. 2010; Bonfils et al. 2013; Delorme et al. 2012; Rameau et al. 2013; Vigan et al. 2012; Bowler et al. 2015). Most of the analyses rely on robust statistical tools based on Monte Carlo simulations or on Bayesian methods and fake planet injection, such as the MESS code described by Bonavita et al. (2012).

The imaged exoplanets generally populate age, mass, and semi-major axis domains (<1 Gyr, >2 M_Jup, >5 au) distinct from those of RV planets usually sought after (around mature stars, all masses, <5–10 au)². As a result, it has not been possible so far to explore the whole range of possible planetary orbits around a given star. It has been shown, though, that young stars can also be probed by RV (e.g., Lagrange et al. 2013) to search for giant planets, and that, combined with DI studies, can allow a full exploration of the parameter space, as far as giant planets are concerned. This brings new observational constraints at different separations, especially because the improving performances of adaptive optics allow now to probe regions closer to the star than previously (with notably SPHERE, GPI, HiCIAO, Beuzit et al. 2006; Macintosh et al. 2014; Suzuki et al. 2010), therefore allowing overlap between the different datasets. This combination of data should significantly improve the robustness of the statistical studies of planetary populations. However, the simple consideration of the best detection limits derived from RV and DI data is not optimal for this analysis. If we consider planets on eccentric orbits in the overlapping separation range, the RV method will preferably detect the planets near their periastron since the RV variations are important there, while DI will detect the planets near their apastron since their projected separation is the largest. So, the actual chance of detecting such planets is closer in this instance to the sum of the detection probability computed for each technique, rather than the best detection probability computed either for RV or for DI. A tool that self-consistently estimates the detection limits of a combined dataset is therefore necessary to derive the optimal detection limits for the growing number of stars currently investigated by RV and DI, and in future astrometry with Gaia (Sozzetti et al. 2014; Perryman et al. 2014).

We present in this article the MESS2 code, an advanced version of MESS that combines RV and DI data obtained at different epochs for individual stars. MESS2 can be used 1) to investigate the giant planet population at all separations in order to test the models of planet population synthesis; and 2) to help improvement of the observational strategy.

We detail the methods for the DI and RV data combination in Sect. 2. In Sect. 3, we apply the MESS2 code to AU Mic, a young and nearby M1 dwarf, and we present a quantitative analysis of systematic biases introduced by using independent analyses instead of our self-consistent tool. We present our conclusions and the associated new perspectives of our results in Sect. 4.

2. Multi-epochs multi-purposes exoplanet simulation system (MESS2): methods

MESS2 allows for derivation of statistical information on the presence of planets around individual stars, using both multi-epoch DI and RV data. It is basically an upgrade of the MESS code (Bonavita et al. 2012), which is optimized for a single dataset. The MESS2 principles are described in the following sections.

2.1. Detection probability with direct imaging data at different epochs

The MESS2 code uses a Monte Carlo approach, similarly to MESS. Step 1 is common to MESS and MESS2, while steps 2 and 3 are specific to MESS2.

• 1.

  Definition of the parameters. We generate a syntheticpopulation of planets using distributions predicted by theoreticalmodels (Mordasini et al. 2009) orsemi-analytical approaches, such as the one by Cumminget al. (2008). The code builds a grid of mass and semi-major axis (a), with a tunable step. The same set of N_gen orbital parameters (eccentricity, inclination, argument of periastron, time of periastron passage) is then generated for each point of the grid. For each of the N_gen sets the projected positions of the planets (x,y) on the plane perpendicular to the line of sight are calculated, using Eqs. (1) to (7) (Bonavita et al. 2012). where X and Y are the orbital coordinates, computed using the Thiele-Innes elements A,B,F, and G described in Eq. (4), ρ is the projected separation, E is the eccentric anomaly, e the eccentricity, ω the argument of periastron, Ω the longitude of the ascending node, i the inclination, M the mean anomaly, T₀ the time of periastron passage, ν the true anomaly, the orbital period of the planet, and t_obs the time of the observations. For each point of the [mass,a] grid, we calculate the position of each synthetic planet for all the available epochs of DI observation.

• 2.

  Definition of the detectability. For every observational epoch, the DI detectability of each planet is tested by comparing its mass to the detection capability of the considered instrument. Here, we use detection limit maps (for further detail see Delorme et al. 2012; Lannier et al. 2016). If a planet is detectable at at least one epoch, it is considered as detectable by DI.

• 3.

  Probability derivation. The probability of detecting a planet of a given mass with a given semi-major axis is derived from the number of detectable planets over the N_gen generated ones. Then, we build a probability map.

2.2. Deriving detection probabilities with RV data only

MESS2 uses a root mean square-based method (RMS) as well as a local power analysis (LPA, Meunier et al. 2012). The RMS method is based on the comparison of the dispersion of synthetic planet RVs with the dispersion of the observed RVs. The LPA method is based on the generation of periodograms of synthetic planet RV time series, which are compared with the periodogram of the observed RV data within given orbital periods.

2.2.1. RMS-based method

• 1.

  Definition of the parameters. The same set of planets used forthe DI data analysis is used. We restrict however the semi-majoraxis range to values that correspond to periods less than half theavailable time baseline.

• 2.

  Definition of the detectability. MESS2 computes synthetic RV curves for each synthetic planet. For a specific planet, given its mass m_planet (in M_Jup), its orbital period (in year), its eccentricity e, and its inclination i, MESS2 computes the half amplitude K (in m s^-1) of the star RV variations (given the stellar mass M_star, in M_Sun), as shown in Eq. (8): (8)Then, the RV time series V of the star due to the presence of the simulated planet orbiting the star are computed at each epoch when RV data were available, using: (9)V₀ is a white noise derived from the standard deviation of the observed temporal series of the star, ω_∗ is the argument of the periastron of the orbit of the star, v_∗ is the true anomaly of the star’s position from the center of mass, and e is the planet eccentricity. N_gen such temporal series are generated for each point of the [mass,a] grid. To test the detectability of the generated planets, we compare the generated time series with the observed one and we apply the condition that no more than 0.15% of the planet detections are actually false positive signals (this condition corresponds to a 3-sigma criterion). In practice, to fulfil the 3-sigma criterion, we generate a large number N of RV time series without planet signal, by adding a white noise on the observed RV measurements, and we count how many times the RMS of the planet-less signal is above X times the RMS of the observed RV; we choose X so that the number of false alarm detections does not exceed 0.15% of N. We therefore consider that a planet is detected if the RMS of its simulated RV time series is higher than the observed one, using the threshold X: (10)

• 3.

  Probability derivation. The probability of detecting planets of a given mass at a given semi-major axis is computed after counting how many planets are detected among the N_gen generated ones, for each considered planetary mass and semi-major axis.

2.2.2. LPA-based method

• 1.

  Definition of the parameters. The free parameters for thegeneration of the periodograms are: i) the period range and thetotal number of periods for the computation of each periodogram;and ii) the orbital and physical parameters of the planets asdescribed in Sect. 2.1. The width of the periodwindows used to compare the periodograms locally to test thedetectability of the planets of period P is fixed to [; ] (the choice of the width of the period window is discussed in Meunier et al. 2012).

• 2.

  Definition of the detectability. The Lomb-Scargle periodograms (Lomb 1976; Scargle 1982; Press & Rybicki 1989) for each generated planet with a given mass and semi-major axis (hereafter called “synthetic periodograms”) are computed, as well as the periodogram using the observed data. To determine the detectability of a planet, the maximum power in a synthetic periodogram is compared with the maximum power in the periodogram of the observed RV data, within the considered period window [; ]. The periods P are calculated from the input semi-major axes using the approximate third Kepler law. We note that since the conversion semi-major axis to period is made using the approximate third Kepler’s law, the planet mass is neglected in this calculation. To mitigate the effects of this approximation when the mass ratio between the star and the planet is high, the value of the exact period (using the non-approximated formula of the third Kepler’s law) is compared to the values of the approximated periods. The closest reference periodogram in terms of period is therefore used to derive the N_gen periodograms for each point of the [mass,a] grid. We estimate that this period allocation leads to semi-major axis shifts of less than 4%. A planet is detected if its maximum power is at least 1.3 times higher than the maximum power in the periodogram of the observed RV data. The choice of this threshold results from an empirical compromise between our need to detect as many small planets as possible and our need to reduce as much as possible the rate of false positives. In order to compare our threshold with the false alarm probability (FAP) formalism, we performed tests on a series of RV made of white noise and with the same calendars as the actual observations. We computed a “local FAP” at each period range, the level corresponding to the 0.1%, 1% etc. “detectability” level. We see that the threshold is equivalent to a local FAP in the range 0.1–1%, depending on the considered periods. In the case of non white noise (e.g. with peaks of stellar origine), the local FAP could be below the peaks in the periodogram that are signal of stellar origine. The 1.3 threshold will be higher than the local FAP in this range, which is more conservative and more realistic. To speed up the computation, the periodograms are only computed for a reference planetary mass m_ref for each semi-major axis and for each input orbital parameter (we call this the reference periodogram, and we arbitrarily choose m_ref = 1M_Jup). We therefore get N_gen reference periodograms for each value of semi-major axis. The periodograms corresponding to the other masses m_p, but to the same orbital parameters, are scaled from the reference periodograms following the formula that we demonstrate in Appendix A: (11)where are the periodograms corresponding to the other masses m_p, and are the reference periodograms computed for 1 M_jup. For both the LPA-based and RMS-based method, only planets with periods lower than two times the time baseline are considered.

• 3.

  Probability derivation. A probability map is computed by adding the number of detected planets among the N_gen generated ones, for each planetary mass and semi-major axis of the grid.

2.2.3. Comparison of the two RV methods used in MESS2

Meunier et al. (2012) compared different methods for determining the detection limits from RV data. They conclude that the LPA-based method is very robust and “provides the most significant improvement on the RMS method”. Indeed, conversely to the RMS method, the LPA method takes into account the temporal structures of the stellar noise, which is in most cases the dominant source of noise. The RMS-based method can however be used to get a first estimate of the detection probability as it requires a much smaller computational time. This is especially useful for cases where large sets of RV measurements are available.

The RV time series are often dominated by the effects of the stellar activity and/or pulsation. Correction for these effects should be applied before using MESS2. In the case of magnetic activity it is possible to perform corrections using the bisector or other activity criteria (as in the case of AU Mic, see below). In the case of pulsations, the RV can be averaged over appropriate timing (see an illustration in Lagrange et al., in prep., in the case of the pulsating A-type star β Pictoris).

2.3. Combining RV measurements and DI data obtained at different epochs

For each point of the (mass, a) grid, the detectability of each generated planet is determined by checking if each planet can be detectable at least in one DI epoch (see Sect. 2.1) or in the RV ones (see Sect. 2.2). The MESS2 code derives the final probability map by counting, for each [mass,a], how many of the N_gen generated planets are detected with either technique.

Notably, we set the detectability criterion at X = 1.44 for the RMS-based method in the case of AU Mic. This criteria corresponds to an equivalent 3-sigma that is appropriate for RV detections. For DI, we adopt the usual 5-sigma criterion; this criterion is realistic far from the star halo (in background limited region) but is probably optimistic in the speckle-dominated region (e.g., Delorme et al. 2012; Rameau et al. 2013; Lannier et al. 2016).

3. Applications: the case of AU Mic

AU Mic (HIP 102409/HD 19748/GJ 803) is a M1V-type star belonging to the β Pictoris moving group (20–26 Myr; in this paper we use 21 Myr, Binks & Jeffries 2014), located at 10 pc (van Leeuwen 2007). Its magnitude is K = 4.5 and its mass is estimated to be 0.61 M_Sun using BT-Settl models (Allard et al. 2012). A debris disk was detected around AU Mic from ≈10–17 to 210 au (Liu et al. 2004; Metchev et al. 2005; Krist et al. 2005), with an edge-on configuration whose inner part is asymmetric (Kalas et al. 2004; Liu 2004; Boccaletti et al. 2015). The gas to dust mass ratio of the disk is low (6:1 Roberge et al. 2005) compared to the typical mass ratio 100:1 of protoplanetary disks, indicating that most of the gas has dissipated. Neither planets nor brown dwarfs have yet been detected in the inner hole of the disk or within it.

3.1. Observations

AU Mic has been observed 10 times by NaCo/ESO between 2004 and 2010, in L′ band (see Table 1), and 26 times by HARPS/ESO over 11 yr (from 2004 to 2015). HARPS provides high-precision RV data (<1 m s^-1). The data were reduced by the instrument data reduction software. Setups of simultaneous Thorium-Argon exposures were used. Thanks to its mass, its close distance (10 pc) and its long RV time base data (11 yr), AU Mic is a good test for MESS2, as it is one of the few cases available for which the parameter space probed by RV and DI significantly overlap. The time sampling of the RV data is not homogeneous: a first set of measurements was taken between 2004 and 2005, and a second set of data was obtained 8 yr later. AU Mic is a very active M1V star with a rotation period P = 4.5 d, a v × sini = 9.3 km s^-1 (Torres et al. 2006), and peak-to-peak V-band light curve amplitudes up to ΔV = 0.10 mag (Messina et al. 2010). Its RV are strongly impacted by the presence of long-lived spots/plages. Yet, the RV can be very efficiently corrected thanks to a very good correlation (Pearson’s coefficient ≈–0.97) between the RV and bisector (see Fig. 3, and other examples of correction of the RV time series using the correlation with the bissector span in Lagrange et al. 2013). Basically, the corrected RV_corr are RV_obs−αBIS, where α is the slope of the [BIS, RV] data. This correction decreases the RV RMS from 156.6 m s^-1 to 33.9 m s^-1. We then use the RV-corrected time series shown in Fig. 3 (bottom). The RV measurements and bisector velocity spans are given in Table A.1.

3.2. Simulations setup

We set the semi-major axis and mass ranges to [0.05; 130] au and [0.05, 80] M_Jup, respectively. The semi-major axis step is logarithmic between 0.05 and 2 au, then linear between 2 and 130 au with a smaller step between 2 and 8 au. We chose this composite sampling to better probe the closer-in region where the RV data are more sensitive. We choose a uniform distribution for the eccentricity between 0 and 0.6, in order to test the detectability of weakly to moderately eccentric planets that would more likely maintain the debris disk dynamically stable. We also choose the inclination range to be [83°, 97°] assuming then that the planets orbit within the plane of the edge-on disk, with the possibility to generate planets immediately above or underneath the disk plane within 7°, following the example of the planets of our Solar System. We use N_gen = 10 000.

3.3. Combining the DI data

We use the giant planet detection probability maps that we had previously derived for each observation epoch (for further detail, see Delorme et al. 2012; Lannier et al. 2016). Figure 1 (top) shows the detection probability map ([0.05, 20] au and [0.5, 14] M_Jup) that we obtain after following the process described in Sect. 2. Figure 2 (middle) illustrates the mass detection limit improvement when using all the available DI data. We note that the low mass detection limits of 1 M_Jup obtained at large separations are in fact upper limits. Indeed, the evolution models by Baraffe et al. (2003) and BT-Settl model atmosphere Allard et al. (2012) do not provide values for masses less than 1 M_Jup.

Compared to the use of one-epoch DI data, our detection probabilities are improved at relatively short separations and for low-mass planets, typically for semi-major axis range 3–15 au and masses under 10 M_Jup. In practice, significant improvements will be obtained if the planet period is ≈4 times the time baseline or less, so that the planets move enough on their orbits to be detectable at one epoch.

As the combination of multiple epochs mostly improves the detection probability for low-mass planets on short separations, in the following we focus on companions with masses between 0.5 and 14 M_Jup and within 20 au.

3.4. Combining the AU MIC DI and RV data

Figure 1 (middle and bottom) shows the probability maps obtained when combining the RV and DI data using the RMS and LPA approaches (middle and bottom, respectively). Figure 2 shows the associated detection limits. Figure 4 shows the comparison of the best detection limits that we could obtain with MESS2 (solid line) and without (dashed line). For the calculation of the detection limits represented by the dashed line, we compute the detection limits obtained with RV data on the one hand and with the most sensitive DI dataset on the other hand, taking the best mass limit between the two of them for each probed separation (in the text we call this approach “the standard combination”, while the “self-consistent combination” refers to the DI and RV data combination within MESS2).

As expected, self-consistently combining DI and RV data rather than using the standard combination, dramatically improves the detection limits at shorter separations compared to the use of one-epoch DI data only, since this parameter range is very well investigated by RV but not yet reachable by the DI instrumental sensitivity. On the contrary, using RV data does not improve the detection of wider-orbit planets, since the RV technique is not sensitive to large separations.

Second, we find that the combination of both RV and DI observations significantly improves the detection limits up to 20 au. For instance, at 5 au, without MESS2, we find a detection limit at 90% probability of ≈8 M_Jup using the LPA approach only, while MESS2 provides a mass limit down to ≈3.5 M_Jup at this separation (LPA method). When an overlap between the RV and DI data is possible, using a standard DI and RV independent detection limits combination becomes obsolete, because this standard combination underestimates the number of generated planets that could be detected, resulting in worse detection limits. This has a direct impact on the planet frequency derived from surveys, since the planet occurrence rate is therefore systematically overestimated within the parameter domain overlapped by DI and RV data. MESS2 gives more robust detection limits and probabilities in this parameter range. For instance, with the LPA approach, twice as many synthetic planets of 2 M_Jup at 4 au are detected than using the standard combination of DI and RV detection limits. The difference between using MESS2 instead of a standard combination of RV and DI detection limits is greater for smaller planets close to the detection limits that are usually used to define the best sensitivity of a data set. At larger separations, there is still an important difference also close to the detection limits, due to the use of multi-epoch DI data (see Fig. 4).

We note that for separations typically between 7 and 16 au, the poor detection limits at 90% probability derived from the standard combination (see Fig. 4) are due to projection effects, even bright companions would be hidden behind the star for more than 10% of their orbit. We also show in Fig. 5 the detection limits that we obtain by testing circular orbits instead of eccentric ones, as the debris disk around AU Mic can prevent the presence of overly eccentric, wide-orbit planets. In the following, we still use the eccentricity range [0, 0.6].

3.5. Scheduling future DI observations

Another independent use of MESS2 is the optimization of the schedule of new DI observations, more precisely to predict the best time intervals and time baselines to obtain new DI data of a star, depending on the semi-major axis and mass ranges.

In order to illustrate this feature of the code we constructed several sets of simulated observations, by combining two copies of the same DI observation (we use the first epoch observation of AU Mic, taken on 2004-06-10) and changing the time span between the two simulated data. This allowed us to assess the changes in detection probability as a function of the time between the single epochs. Figure 6 shows the probability maps when further DI data is taken 5 and 25 yr later. We would like to highlight that we see poor-gain vertical lines that move towards larger semi-major axis when the time lapse increases. They are associated to the revolution periods of the planets. For instance, planets located at 2.5 au and 7.3 au cover their entire orbit respectively in 5 and 25 yr, so there is no benefit to re-observe them after a full revolution. Figure 7 illustrates the detection probability gain obtained when combining different DI data separated by 0.5, 1, 2, 3, 5 or 8 yr, over different time baselines, for planets between 0.5 and 14 M_Jup located in the restrictive semi-major axis range [4–8] au. The gain is defined as the median value of the detection probability after the combination of several DI data, over the median value of the detection probability using one-epoch data, within the planet masses and semi-major axes restrictions. Of course, observing a star as often as possible gives the best gains, but a very high observational frequency is not necessary; observing only three times AU Mic over ten years gives almost the same gain as observing the target more frequently, if we search for planets with masses 0.5–14 M_Jup at 4–8 au, with appropriate time sampling. Indeed, the plateau that we observe, around an improvement of ≈90% after observing the star several times, is reachable with a 10-yr time baseline, for DI observations separated by less than 5 yr.

4. Concluding remarks and perspectives

We present the multi-epochs multi-purposes exoplanet simulation system (MESS2) code, a statistical tool that combines multi-technique data. Our aim with MESS2 was 1) to properly constrain giant planet populations from the inner to outer regions of stars (using RV and DI techniques) to eventually constrain the models of planet population synthesis; and 2) to optimize the scheduling of DI observations. A MESS2 application is presented in this paper in the case of the late-type star AU Mic. Another application will be presented in the case of a pulsating early-type star, β Pictoris, in a forthcoming paper (Lagrange et al., in prep.).

In the case of AU Mic, located at 10 pc, we were able to bridge the gap between RV and DI data and constrain the presence of giant planets with masses ≥2 M_Jup (using probabilities of 60%) from a fraction of au to several tens of AU, and ≤1 M_Jup outside the 2–4 au range. The detection limits obtained after the self-consistent combination of RV and DI data with MESS2 are significantly better than those obtained with the standard combination of RV and DI independent detection limits. In particular, between 7 and 16 au, we show that more than 90% of planets down to 2 M_Jup are detectable. On the other hand, the standard combination of RV and DI independent detection limits cannot exclude, with a 90% probability, the presence of brown dwarfs of any mass at these separations. In addition to the fact that MESS2 provides more constrained and more robust detection limits, not using a self-consistent tool such as MESS2 for AU Mic leads to an overestimation of the planet occurrence rate. Additional RV data and/or future Gaia astrometric data as described by Sozzetti et al. (2014) could better constrain the intermediate separations, that is, to less than a few astronomical units.

The computed detection probabilities depend on the configuration of the planets’ orbits (i.e., the inclination with respect to the line of sight, and the position angle of the orbit with respect to the North). First, the RV and DI methods are sensitive in opposite ways to the planets inclinations. If the system is observed edge-on (as AU Mic), then combining sets of DI and RV data greatly improves the probability of detecting planets at all separations. On the contrary, if the system is observed pole-on, then the RV measurements do not constrain the detection probabilities. Moreover, if one assumes circular and pole-on orbits, there will be no advantage of combining all DI data. The intermediate cases between those extremes are more common, like HR 8799 (25° with respect to the plane of the sky, Contro et al. 2015), TW Hydrae (7° for the outer disk and 4.3° for the inner disk, Setiawan et al. 2008; Pontoppidan et al. 2008), or HD 141569 (51°, van der Plas et al. 2015). Figure 8 shows that the different inclinations affect the detection limits in different ways. Second, it is also possible, within the MESS2 code, to take into account the position angle of the system by constraining the longitude of ascending mode. Constraining this orbital parameter is relevant in cases where a structure such as a circumstellar disk is detected in the considered DI data and affects the data sensitivity. Taking into consideration the position angle of the system does not affect the detectability of RV planets, while the impact on the detectability of DI planets can be strong. Applying MESS2 on β Pictoris data is a good test for constraining the position angle of the orbits of the synthesized planets (see Lagrange et al., in prep.).

The impact of each detection technique also depends on the stellar properties. For instance, the younger the stars are, the more sensitive DI is. This is not necessary true for RV, however; it depends on the stellar activity, pulsation and rotation levels.

Finally, the combination of RV and DI data when RV trends are observed can be used to constraint the orbital and physical parameters of the potential planet thatwould produce this trend. Indeed, either the planet is detected within the DI data and thus its parameters are constrained, or it is not imaged and then the DI detection limits can be used to put upper limits on its physical properties. This application of MESS2 will be the object of a future paper (Bonavita et al., in prep.).

---

Acknowledgments

We thank the staff of ESO-VLT and ESO-La Silla for their support at the telescope. We acknowledge support from the French National Research Agency (ANR) through the GuEPARD project grant ANR10-BLANC0504-01 and GIPSE project grant ANR-14-CE33-0018. We thank Pascal Rubini for his work on the software SAFIR.

References

Appendix A: Additional table

Appendix

We demonstrate here the formula of proportionality that connects the values of periodograms computed using two different planet masses, at the same semi-major axis. Considering two planets with masses m₁ and m₂, respectively, orbiting the same star of mass m_∗ at the same period and having the same orbital inclination i, the half-amplitudes of their RV curves are respectively: Then, the ratio of these half-amplitudes is: (B.3)

The Lomb-Scargle periodograms give the spectral power of the RV measurements as a function of the planet periods. Mathematically, the Lomb-Scargle periodograms are the square of the Fourier Transform of the RV measurements, we have: Then, the ratio between the values of two periodograms for the period P can be written: (B.6)

All Tables

All Figures