Combining direct imaging and radial velocity data towards a full exploration of the giant planet population

TL;DR

This paper introduces MESS2, a Monte Carlo simulation tool that combines direct imaging and radial velocity data to improve detection and statistical analysis of giant exoplanets across all orbital separations.

Contribution

The paper presents MESS2, a novel tool that integrates different detection methods to accurately estimate giant planet populations around stars.

Findings

Detection limits are significantly improved at intermediate separations.

Independent analysis of RV and DI overestimates planet occurrence rates.

MESS2 enables correct planet occurrence rate estimation and optimizes future observations.

Abstract

Thanks to the detections of more than 3000 exoplanets these last 20 years, statistical studies have already highlighted some properties in the distribution of the planet parameters. Nevertheless, few studies have yet investigated the planet populations from short to large separations around the same star since this requires the use of different detection techniques that usually target different types of stars. We wish to develop a tool that combines direct and indirect methods so as to correctly investigate the giant planet populations at all separations. We developed the MESS2 code, a Monte Carlo simulation code combining radial velocity and direct imaging data obtained at different epochs for a given star to estimate the detection probability of giant planets spanning a wide range of physical separations. It is based on the generation of synthetic planet populations. We apply MESS2 on a young M1-type, the nearby star AUMic observed with HARPS and NACO/ESO. We show that giant planet detection limits are significantly improved at intermediate separations (~20au in the case of AUMic). We show that the traditional approach of analysing independently the RV and DI detection limits systematically overestimates the planet detection limits and hence planet occurrence rates. The use of MESS2 allows to obtain correct planet occurrence rates in statistical studies, making use of multi-epoch DI data and/or RV measurements. We also show that MESS2 can optimise the schedule of future DI observations.

Figures

Tables

Table 1. Table 1: Periods of observation for the AU Mic DI data.

Obs. period	run	number of observations
June, July, Sept. 2004	073.C-0834(A)	7
July 2010	085.C-0675(A)	1
Sept. 2010	085.C-0277(B)	1

Table 2. Table 2: AU Mic initial RV measurements with errors, and bisector velocity span with errors in km.s -1 , for the corresponding date of observation.

BJD-2453000	RV	RV errors	bisector velocity span	bisector velocity span errors
157.898424	0.022713	0.002816	-0.063852	0.007041
201.823450	-0.044461	0.003343	0.039650	0.008356
468.892370	-0.149082	0.003154	0.071819	0.007885
469.843534	-0.193213	0.002877	0.104851	0.007191
521.894368	-0.421558	0.003464	0.240802	0.008661
551.803998	-0.038124	0.003970	0.000000	0.009926
593.622139	0.145622	0.003777	-0.096523	0.009442
3568.505024	-0.290126	0.003997	0.229599	0.009993
3568.515706	-0.275093	0.004017	0.219286	0.010044
3569.494659	-0.032413	0.003461	-0.035353	0.008653
3569.505549	-0.056961	0.003481	0.030390	0.008703
3570.560455	0.198430	0.003022	-0.138663	0.007555
3570.571449	0.204715	0.003110	-0.152526	0.007775
3772.914293	-0.065683	0.003395	0.036329	0.008488
3772.924248	-0.069677	0.003422	0.063948	0.008556
3773.911762	-0.017799	0.002991	0.028629	0.007477
3773.926196	-0.015262	0.003064	0.022738	0.007661
3794.882288	-0.130789	0.003228	0.062227	0.008070
3795.885873	0.072862	0.003100	-0.034385	0.007749
3797.857541	0.018862	0.003558	-0.036264	0.008894
3844.801157	0.174065	0.003308	-0.114677	0.008270
3844.811169	0.172775	0.003392	-0.123908	0.008479
3982.534763	-0.080822	0.003252	0.032365	0.008129
3982.544392	-0.077453	0.003327	0.042436	0.008318
4333.530396	0.100939	0.002592	-0.085290	0.006480
4333.541066	0.104374	0.002721	-0.070044	0.006802

Equations

$$x=AX+FY$$

$$y=BX+GY$$

$$X=\cos E-e$$

$$Y=\sqrt{1-e^2}\sin E$$

$$\rho =\sqrt{x^2+y^2}$$

$$A=a(\cos \omega \cos \Omega -\sin \omega \sin \Omega \cos i)$$

$$B=a(\cos \omega \sin \Omega +\sin \omega \cos \Omega \cos i)$$

$$F=a(-\sin \omega \cos \Omega -\cos \omega \sin \Omega \cos i)$$

$$G=a(-\sin \omega \sin \Omega +\cos \omega \cos \Omega \cos i)$$

$$M=(\frac{t_{obs}-T_0}{\mathcal{P}})2\pi$$

$$E_0=M+e\sin M+\frac{e^2}{2}\sin 2M$$

$$M_0=E_0-e\sin E_0$$

$$E=E_0+\frac{M-M_0}{1-e\cos E_0}$$

$$\tan \frac{\nu }{2}=\sqrt{\frac{1+e}{1-e}}\tan \frac{E}{2}$$

$$K=\frac{28.432\times m_{\mathrm{planet}}\times sin(i)}{M_{\mathrm{star}}^{2/3}\times {\mathcal{P}}^{1/3}\times \sqrt{1-e^2}}$$

$$V(t)=V_0(t)+K(t)\times (cos({\omega }_{\ast }(t)+v_{\ast }(t))+e\times cos({\omega }_{\ast }(t)))$$

$$\textrm{RMS}_{\textrm{generated RV}}>$X$\times\textrm{RMS}_{\textrm{observed RV}}$$

$$\mathcal{P}(P)={\mathcal{P}}_{\mathrm{ref}}(P)\times (\frac{m_{\mathrm{p}}}{m_{\mathrm{ref}}}\times (\frac{m_{\mathrm{star}}+m_{\mathrm{ref}}}{m_{\mathrm{star}}+m_{\mathrm{p}}}{)}^{2/3}{)}^2$$

$$K_1=(\frac{2\pi G}{\mathcal{P}}{)}^{\frac{2}{3}}\frac{m_1\sin i}{(m_{\ast }+m_1{)}^{\frac{2}{3}}}(1-e^2{)}^{-\frac{1}{2}}$$

$$K_2=(\frac{2\pi G}{\mathcal{P}}{)}^{\frac{2}{3}}\frac{m_2\sin i}{(m_{\ast }+m_2{)}^{\frac{2}{3}}}(1-e^2{)}^{-\frac{1}{2}}$$

$$\frac{K_1}{K_2}=\frac{m_1}{m_2}(\frac{m_{\ast }+m_2}{m_{\ast }+m_1}{)}^{\frac{2}{3}}$$

$${\mathcal{P}}_1(P)\propto K_1^2$$

$${\mathcal{P}}_2(P)\propto K_2^2$$

$$\frac{{\mathcal{P}}_1(P)}{{\mathcal{P}}_2(P)}=(\frac{m_1}{m_2}\times (\frac{m_{\ast }+m_2}{m_{\ast }+m_1}{)}^{\frac{2}{3}}{)}^2$$

Full text

*Aims. *We wish to develop a tool that combines direct and indirect methods so as to correctly investigate the giant planet populations at all separations.

*Methods. *We developed the MESS2 code, a Monte Carlo simulation code combining radial velocity and direct imaging data obtained at different epochs for a given star to estimate the detection probability of giant planets spanning a wide range of physical separations. It is based on the generation of synthetic planet populations.

*Results. *We apply MESS2 on a young M1-type, the nearby star AU Mic observed with HARPS and NACO/ESO. We show that giant planet detection limits are significantly improved at intermediate separations ($\approx$20 au in the case of AU Mic). We show that the traditional approach of analysing independently the RV and DI detection limits systematically overestimates the planet detection limits and hence planet occurrence rates. The use of MESS2 allows to obtain correct planet occurrence rates in statistical studies, making use of multi-epoch DI data and/or RV measurements. We also show that MESS2 can optimise the schedule of future DI observations.

††offprints: J. Lannier,

1 Introduction

More than 3000 extra-solar planets have been detected to date, mostly with radial velocity (RV) and transit techniques (http://www.exoplanet.eu). Beside these prolific indirect detection methods, about 30 young, wide-orbit, giant exoplanets or brown dwarfs have been discovered using direct imaging techniques (DI). The 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; Chauvin et al. 2010; Johnson 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 MJup, ¿5 au) distinct from those of RV planets usually searched (around mature stars, all masses, ¡5-10 au)111Transiting exoplanets are found on average on shorter separations.. 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 showed, though, that young stars can also be probed by RV (eg, 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 to detect 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, DI and in the 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, 2) to help improving the observational strategy.

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

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

MESS2 allows to derive 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 optimised 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.

Definition of the parameters. We generate a synthetic population of planets using distributions predicted by theoretical models (Mordasini et al. 2009) or semi-analytical approaches as the one by Cumming et al. (2008).

The code builds a grid of mass and semi-major axis (a), with a tunable step. A 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 the equations (1) to (7) (Bonavita et al. 2012).

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

$X$ and $Y$ are the orbital coordinates, computed using the Thiele-Innes elements $A,B,F,G$ described in Eq. (4). $\rho$ is the projected separation. $E$ is the eccentric anomaly, $e$ the eccentricity, $\omega$ the argument of periastron, $\Omega$ the longitude of the ascending node, $i$ the inclination, $M$ the mean anomaly, $T_{0}$ the time of periastron passage, $\nu$ the true anomaly, ${\cal P}$ 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. 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 least at one epoch, it is considered as detectable by DI. 3. 3.

Probability derivation. The probability to detect 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 solely

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 planets RVs with the dispersion of the observed RVs. The LPA method is based on the generation of periodograms of synthetic planet RV time series, that are compared with the periodogram of the observed RV data within given orbital periods.

2.2.1 RMS-based method

Definition of the parameters. The same set of planets used for the DI data analysis is used. We restrict however the semi-major axis range to values that correspond to periods less than half the available time baseline. 2. 2.

Definition of the detectability. MESS2 computes synthetic RV curves for each synthetic planet. For a specific planet, given its mass $m_{\rm planet}$ (in $M_{Jup}$), its orbital period ${\cal P}$ (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_{\rm star}$, in $M_{\rm Sun}$), as shown in Eq. 8:

[TABLE]

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:

[TABLE]

$V_{0}$ is a white noise derived from the standard deviation of the observed temporal series of the star. $\omega_{*}$ is the argument of the periastron of the orbit of the star, $v_{*}$ the true anomaly of the star’s position from the centre of mass, and $e$ 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, and 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$:

[TABLE] 3. 3.

Probability derivation. The probability to detect 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

Definition of the parameters. The free parameters for the generation of the periodograms are: i) the period range and the total number of periods for the computation of each periodogram; and ii) the orbital and physical parameters of the planets as described in Section 2.1. The width of the period windows used to compare the periodograms locally to test the detectability of the planets of period $P$ is fixed to [0.75$\times{\cal P}$;1.25$\times{\cal P}$] (the choice of the width of the period window is discussed in Meunier et al. 2012). 2. 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 [0.75$\times{\cal P}$;1.25$\times{\cal P}$]. The periods P are calculated from the input semi-major axes using the approximate third Kepler law222Note 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. Note that the rate of false positives is very difficult to estimate: as discussed in Meunier et al. (2012), the use of false alarm probability, which is sometimes used, has a limited interest, as it cannot take into account the temporal structure of the noise333In order to compare our threshold with the 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 computed only for a reference planetary mass $m_{\rm ref}$ for each semi-major axis and for each input orbital parameter (we call it the reference periodogram, and we arbitrarily choose $m_{\rm ref}=1$ MJup). We therefore get $N_{gen}$ reference periodograms for each value of semi-major axis. The periodograms corresponding to the other masses $m_{\rm p}$, but to the same orbital parameters, are scaled from the reference periodograms following the formula that we demonstrate in appendix A:

[TABLE]

where ${\cal P}$ are the periodograms corresponding to the other masses $m_{\rm p}$, and ${\cal P}_{{\rm ref}}$ are the reference periodograms computed for 1 Mjup.

For both the LPA-based and RMS-based method, only planets with periods lower than two times the time baseline are considered. 3. 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 to determine 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, 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. 2016, in prep., in the case of the pulsating A-type star $\beta$ 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 Section 2.1) or in the RV ones (see Section 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.

Note: 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 which 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 (eg 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, that belongs to the $\beta$ Pictoris moving group (20-26 Myr, in this paper we will 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 $\approx$10-17 to 210 au (Liu et al. 2004; Metchev et al. 2005; Krist et al. 2005), with an edge-on configuration whose the 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 been detected in the inner hole of the disk or within it yet.

3.1 Observations

AU Mic has been observed 10 times by NaCo/ESO between 2004 and 2010, in L’ band (see Tab. 1), and 26 times by HARPS/ESO over 11 years (from 2004 to 2015). HARPS provides high precision RV data (¡1m/s). 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 years), 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 has been taken between 2004 and 2005, and a second set of data was obtained 8 years later. AU Mic is a very active M1V star with a rotation period P = 4.5 d, a $v\times\sin i=9.3\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}$ (Torres et al. 2006), and peak-to-peak V-band light curve amplitudes up to $\Delta$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 $\approx$-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}-\alpha BIS$, where $\alpha$ is the slope of the [BIS, RV] data. This correction decreases the RV RMS from 156.6 m/s to 33.9 m/s. We then use the RV corrected time series shown in Fig. 3 (bottom). The RV measurements and bisector velocity spans are given in Table 2.

3.2 Simulations Setup

We set the semi-major axis and mass ranges to respectively [0.05;130] au and [0.05,80] $M_{Jup}$. 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 just above or underneath the disk plane within 7∘*, following the example of the planets of our Solar System. We use $N_{gen}=10000$.

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] MJup) that we obtain after following the process described in Section 2. Figure 2 (middle) illustrates the mass detection limit improvement when using all the available DI data. Note that the low mass detection limits of less than 1 MJup 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 MJup.

Compared to the use of a 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 MJup. In practice, significant improvements will be obtained if the planet period is $\approx 4$ 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 will 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, combining self-consistently DI and RV data rather than using the standard combination, dramatically improves the detection limits at shorter separations compared to the use of a 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 $\approx$8 MJup using the LPA approach only, while MESS2 provides a mass limit down to $\approx$3.5 MJup 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 under-estimates 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 was 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 MJup 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 higher 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, that is due to the use of multi-epoch DI data (see Fig. 4).

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 during more than 10% of their orbit.

Note also that we 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 too 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 optimisation 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 a second DI data is taken 5 and 25 years later. Note 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 years, 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 years, over different time baselines, for planet between 0.5 and 14 MJup 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 a 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 than observing the target more frequently, if we search for planets with masses 0.5-14 MJup at 4-8 au, with appropriate time sampling. Indeed, the plateau that we observe around an improvement of $\approx$90% at observing the star several times, is reachable with a 10-year time baseline, for DI observations separated by less than 5 years.

4 Concluding remarks and perspectives

We presented the MESS2 code (Multi-epochs multi-purposes exoplanet simulation system), 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 optimise 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, $\beta$ Pictoris, in a coming paper (Lagrange et al. 2016, in prep.).

In the case of AU Mic, located at 10 pc, we could bridge the gap between RV and DI data and constrain the presence of giant planets with masses $\geq$2 MJup (using probabilities of 60%) from a fraction of au to several tens of AU, and $\leq$1 MJup 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 showed that more than 90% of planets down to 2 MJup are detectable. On the opposite, 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 overestimate the planet occurence rate. Additional RV data and/or future GAIA astrometric data as described by Sozzetti et al. (2014) could better constrain the intermediate separations, i.e. less than a few astronomical units.

The detection probabilities computed 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 together greatly improves the probability to detect 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 gain at combining all DI data. The intermediate cases between those extremes are more common, like HR8799 (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 HD141569 (51∘, van der Plas et al. 2015). Fig. 8 shows that the different inclinations affect differently the detection limits. 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 $\beta$ Pictoris data is a good test for constraining the position angle of the orbits of the synthesized planets (see Lagrange et al. 2016, 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, but this is not necessary true for RV (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 that would produce this trend. Indeed, either the planet is detected within the DI data and thus its parameters are constrained, either 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. 2016).

Acknowledgements.

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.

Appendix A

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. Let’s consider two planets respectively with masses $m_{\rm 1}$ and $m_{\rm 2}$, orbiting the same star of mass $m_{\rm*}$ at the same period ${\cal P}$ and having the same orbital inclination $i$. The half-amplitudes of their RV curves are respectively:

[TABLE]

Then, the ratio of these half-amplitudes is:

[TABLE]

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:

[TABLE]

Then, the ratio between the values of two periodograms for the period $P$ can be written:

[TABLE]