© ESO, 2012

1. Introduction

Many statistical studies have been made using information coming from more than a decade of extensive searches for exoplanets, trying to answer questions either related to the properties of those objects, such as the mass, orbital period and eccentricity (Lineweaver & Grether 2003; Cumming et al. 2008), or to the relevance of the host star characteristics (mass, metallicity and binarity) on the final frequency and distribution of planetary systems (see Fischer & Valenti 2005; Santos et al. 2004; Johnson et al. 2007). Because the most successful techniques (radial velocity and transit) have focused on the inner (≤5 AU) environment of main-sequence solar-type stars, most of the available information on the frequency of planets concerns this class of stars.

Recent discoveries of young distant planetary mass objects with direct imaging (see e.g. Marois et al. 2008; Kalas et al. 2008; Lagrange et al. 2009) are giving us a first hint on the potential of the direct detections in the exploration of the outer region of the planetary systems, also raising many questions about how these objects could form (see Absil & Mawet 2009). This defines the niche of the next generation of high-contrast imaging instruments like the Gemini Planet Imager (GPI: Macintosh et al. 2007) and VLT/SPHERE (Spectro-Polarimetric High-contrast Exoplanet REsearch: Beuzit et al. 2008). These instruments will likely allow us to extend this systematic characterization to larger scales ( ≥ 10 AU). Owing to practical limitations (inner working angle, best contrast achievable), these instrument will focus on warm giant planets on orbits far away from their stars, paving the path for the ELTs facilities. A wide range of planetary masses and separations, down to the rocky planets (and, in very favourable cases reaching the habitable zone), will be explored with 30−40 meter-class telescopes, finally allowing an overlap between the discovery spaces of direct and indirect techniques.

In this context it is useful and crucial to predict the performances of the forthcoming instruments, not only in terms of number of expected detections, but also trying to figure out what the explored parameter space will be and even the possible synergies between different discovery techniques.

Here we present our Monte Carlo simulation code MESS, whose aim is to provide a flexible and reliable tool for the statistical analysis and prediction of the results of planet searches.

The code produces synthetic planet populations, deriving all physical parameters of these planets together with the observables that can be compared with the predicted capabilities of existing or planned instruments. These comparisons allow us to derive subsets of fully characterized detectable planets as well as a snapshot of what the evolution of the sample of detected planets would be in the next years.

A detailed description of the code and of all the assumptions that constitute its basis is given in Sect. 2, while in Sect. 3 we present the different operation modes of the code and their applications. Although MESS was built and has so far been applied only to analyse and/or predict the results of direct imaging surveys, an extension of the code to different techniques is planned. The first attempt in this direction is presented in Sec. 3.3. Conclusions and suggestions for further work will be finally drawn in Sect. 4.

2. MESS (multi-purpose exoplanet simulation system)

Over the past years, several groups (Kasper et al. 2007; Lafrenière et al. 2008; Chauvin et al. 2010; Nielsen & Close 2010) initiated statistical analyses to constrain the physical and orbital properties (mass, period, eccentricity distributions) of the giant planet population at large separations. These groups developed statistical analysis tools appropriate for exploiting the performances of deep imaging surveys. They also tested the consistency of various sets of parametric distributions of planet properties, using the specific case of a null detection. The first assumption of these tools is that planet mass and period distributions coming from the statistical results of radial velocity (RV) studies at short period (see e.g. Lineweaver & Grether 2003; Cumming et al. 2008) can be extrapolated and normalized to obtain information on more distant planets. Despite the model-dependency on the mass predictions, the approach is attractive for exploiting the complete set of detection performances of the survey and characterizing the outer portions of exoplanetary systems.

With all of this in mind, we aimed to go a step further by creating a multi-purpose exoplanet simulation system (hereafter MESS) that is to be applied also to techniques other than direct imaging and also uses the information from the planetary formation theories.

The code is written in IDLr¹, its source and some examples are available in electronic format at the CDS.

Uptadet versions of the code, and a step by step tutorial are also available online at www.messthecode.com.

The basic operations performed by the code are the following:

• 1.

  it generates a synthetic population of planets, including allorbital elements, either using the planet mass and perioddistributions coming the statistical results of RV studies or theoutcome of the planetary formation theories;

• 2.

  taking into account the characteristics of the host star and of the planetary orbit, it calculates all observable quantities needed for the comparison with the instrument performances, such as RV and astrometric signal, planet/star contrast, degree of polarization, etc.;

• 3.

  given the detection capability relation of an instrument, either already available or planned, it selects a sub-sample of fully characterized detectable planets, whose characteristics can then be analysed.

The code then assumes a given star population, a planet population with associated physical and orbital properties based on a theoretical or semi-empirical approach, the corresponding observables for different observing techniques, and finally generates a synthetic population of planets to be compared with the instrumental detection performances. Each step is described below.

2.1. Star population

The first input of the MESS is a sample of N_Star stars, that were targeted for planet searches or that are part of a sample for future observations. Various stellar parameters are assumed to be known, such as the apparent magnitude, the distance, the luminosity, the spectral type, the mass, the age, the metallicity, etc.

Figure 1 shows the characteristics of a sample of 600 nearby (d < 20 pc) stars selected from the Hipparcos catalogue (Perryman & ESA 1997), which were used to build the synthetic population showed in Fig. 3.

2.1.1. Binarity module

MESS also gives the possibility of taking into account the presence of one (or more) additional stellar companions in the analysis. If a star in the sample is flagged as a binary, the code uses the information about the binary orbit (if available) to compute the critical semi-major axis for the dynamical stability of the system. This sets the limiting value that the semi-major axis of a planet can attain and still maintain its orbital stability, as a function of the mass-ratio and orbital elements of the binary, as shown by Holman & Wiegert (1999).

Both the case of circumstellar (or satellite S-type) and circumbinary (or planet P-type) orbit are considered, and the critical semi-major axis is computed using Eqs. (1) and (2) respectively, from Holman & Wiegert (1999), In both equations, a_c is critical semi-major axis, μ = M₁/(M₁ + M₂), a_b and e_b are the semi-major axis and eccentricity of the binary, and M₁ and M₂ are the masses of the primary and secondary stars, respectively. If not available from literature, the eccentricity is assumed to be e_b = 0.36, reported as the mean value for the eccentricity of a binary system, by Duquennoy & Mayor (1991). If the value of the semi-major axis is not available, the code estimates it as a_b = 1.31   ρ (arcsec) d (pc) (see Fischer et al. 2002; Duquennoy & Mayor 1991)².

Note that in the first case (S-type orbit), a_c sets the maximum value that the semi-major axis of a planet can assume before compromising the stability, while it represents the minimum value of the semi-major axis of a stable planet for a P-type orbit.

2.2. Planet population

The core of the code is the generation of synthetic planets that are fully characterized both in terms of orbital parameter, and physical characteristics. Depending on the goal of the study, one can choose between a semi-empirical approach or a theoretical approach. These different approaches make the code suitable to constrain the planet properties under different assumptions, but also to test model predictions.

If the theoretical approach is chosen, masses and period values selected from a synthetic population provided by the output of the planetary formation models (see e.g. Mordasini et al. 2009) are given as input. In this case all orbital characteristics are provided as well, together with the physical properties of each planet, so that no random generation is needed, and the code only evaluates the observables and compares them to the provided detectability relations. Different populations of planets obtained assuming different stellar masses and metallicity values can be selected according to the characteristics of the real star in the sample, to take into account the effects of the stellar characteristics on the planet-formation processes³.

The semi-empirical approach uses the power-law distributions in Eqs. (3) and (4) for the mass and semi-major axis of the planets as retrieved from the statistical analysis of the properties of the planets discovered so far to generate a seed population of N_seed values of masses and periods (see Sect. 2.2.1) The user can also set a pre-determined grid of masses-periods and feed it to the code, without any assumption on the distributions. This would be the case if, for example, boundaries on the mass/semi-major axis space where planets can form are to be set using the outcomes of a formation model (see e.g. Mordasini et al. 2010), excluding from the sample those planets that are not compatible with the theory.

If the semi-empirical approach is used, mass, orbital parameters, temperature, and radius of the planets are obtained based on the assumption described in the next sections.

2.2.1. Mass-period seed generation

If the semi-empirical approach is chosen, the power-law distributions are fed to the Monte-Carlo core of the code, which randomly generates a fixed number of mass-period pairs. Both the planetary mass and period ranges can be given as inputs, together with the power-law exponents. In a typical setup, the power-law exponents are assumed to be α =  −1.31 and β =  −0.74 respectively, according to Cumming et al. (2008). The planetary masses span the range 0.6 M_Earth and 15 M_Jup, and the periods (P) are chosen to be between 2.5 days and 350 years (corresponding to 50 AU for 1   M_⊙ star).

A scaling of the planetary mass and even of the period with the stellar mass can be also introduced, according to recent results (e.g. Lovis & Mayor 2007). In addition, a dependence of the planet frequency on the stellar metallicity may also be considered (see Fischer & Valenti 2005).

2.2.2. Evaluation of the orbital parameters

For each mass-period pair in the seed generation, the code evaluates the semi-major axis computed using Kepler’s third law, using the mass of each star in the input sample. Then it generates N_gen values of all orbital parameters: eccentricity (e), inclination (i), longitude of periastron (ω), longitude of ascending node (Ω), and time of periastron passage (T₀). By default, all these parameters are randomly generated following uniform distribution⁴. The eccentricity distribution is cut at e = 0.6 as suggested by the results of the RV surveys (see Cumming et al. 2008). This also allows us to control a possible bias towards high-eccentricity planets that could affect the results of direct imaging surveys. A full discussion of the impact of the eccentricity distribution on the simulations results is given in Sect. 3.4.

The date of observation is also required. If not available from the real data, an epoch of observation, t_obs, is generated over a time-span chosen according to the considered instrument.

The code also offers the possibility to fix each orbital parameter to known or predicted values for all the planets in the population.

The coordinates x and y of the projected orbit on the plane perpendicular to the line of sight are finally computed using the ephemeris formulae of Heintz (1978), reported in Eqs. (5) to (7). where X and Y are the coordinates of the orbit (Eq. (6)), ρ is the projected separation, and A,B,F,G are the Thiele-Innes elements, which can be obtained from the classical ones (the semi-major axis a, ω, Ω, and i) using Eq. (8): (8)In these equations, E is the eccentric anomaly (obtained from the mean anomaly M (Eq. (9)) using Eq. (10)) and ν the true anomaly (Eq. (11)): The projected separation, ρ (in arcsec), can be obtained either using Eq. (7) or Eq. (12) (which also gives an estimate of the radius vector r), then dividing by the star distance, (12)

2.2.3. Planet temperature

Because we aim at considering both the thermal and reflected flux of the planets, we need two different estimates of the temperature. The first one is the internal temperature, T_int, obtained from the evolutionary models (see e.g. Baraffe et al. 2003). The second one is the equilibrium temperature, T_eq, obtained through Eq. (13) (from Sudarsky et al. 2003) (13)where L _∗  is the star luminosity. The Bond albedo A_B is assumed to be 0.35 in the J-band (Jupiter value, see Hanel et al. 1981) and it is randomly generated between 0.3 and 0.52 in the visible (the latter being the Jupiter albedo in V-band, see Sudarsky et al. 2003).

Our final assumed value for the effective temperature of the planet T_eff is given by (14)

2.2.4. Planet radius

MESS uses the approach developed by Fortney et al. (2007) to evaluate the planetary radius. Practically the radius is assumed to depend on the planet mass, with the following recipes:

• 1.

  For Jupiter-like planets (M ≥ 100   M_Earth), an interpolation is performed within the published values given by Fortney et al. (2007). Values of age and distance of each star are entered, yielding a value for R_Gas. A core mass of 10   M_Earth is assumed.

• 2.

  Equations (15) and (16) from Fortney et al. (2007) are used for the smallest planets (M ≤ 10   M_Earth). These are either (15)or (16)for ice/rock and rock/iron planets, respectively. In these equations, R is in R_Earth and M is in M_Earth, while imf is the ice mass fraction (1.0 for pure ice and 0.0 for pure rock) and rmf is the rock mass fraction (1.0 for pure rock and 0.0 for pure iron). In the typical MESS setup, the ice/rocky or rocky/iron fraction is set to 0.3 (50% of chance for a planet being mainly icy or rocky).

• 3.

  Finally, predictions are uncertain for the Neptune-like planets, where the transition between the two relations described above should occur. The most sensible approach seems to be to fit the mass-radius relation of the solar system in the same mass-range (10 − 40   M_Earth). This procedure provides a good agreement with the radii of Uranus and Neptune and of the few transiting Neptunes confirmed so far (as listed by The Extrasolar Planet Encyclopaedia⁵ see Fig. 2).

The resulting mass-radius relations are shown in Fig. 2, with the data corresponding to the planets discovered with the transit technique and the planets from our solar system overplotted for comparison.

2.3. Predicted observables

After obtaining the full set of orbital and physical parameters of the planets, the code then provides an estimate of observable quantities such as the luminosity contrast or the degree of polarization, which are needed for direct observations, but also quantifies the indirect effects of the presence of the planet, providing a measure of the semi-amplitude of radial velocity (RV) and the astrometric signal.

2.3.1. Planet/star contrast

MESS gives an estimate of both the intrinsic and reflected flux in the selected band for each planet. Throughout the paper we will refer to the planets whose luminosity is dominated by the intrinsic contribution as self-luminous or warm planets, as opposed to the cold planets, for which the reflected light provides most of the contribution to the planet/star contrast.

The intrinsic emission is estimated using the prediction of evolutionary models at the age of the star (assumed to be also the age of the system). To this purpose two classes of models can be considered, based on different assumptions on the initial conditions: Hot start models (Chabrier et al. 2000; Baraffe et al. 2003; Saumon & Marley 2008), which consider an initial spherical contracting state; and core accretion models (Marley et al. 2007; Fortney et al. 2008), which couple planetary thermal evolution to the predicted core mass and thermal structure of a core-accretion planet-formation model.

In the following, we only consider the results obtained using the hot start models for the nearby sample. However, the problem of the initial condition and the uncertainties on the stellar ages are among the main limitations for young stellar samples, not only for our code, but also for any kind of study that uses the same kind of approach (see e.g. Chauvin et al. 2010; Bonavita et al. 2010, for a detailed discussion). These limitations also apply to the theoretical approach if the evolutionary models are used to evaluate the planet intrinsic luminosity and radii produced by the models, as in Mordasini et al. (2010, 2011).

For the evaluation of the reflected light we scaled the Jupiter value according to the planet radius (expressed in Jupiter radii), semi-major axis, albedo and illuminated fraction of the planet. This last contribution was computed through a phase-dependent term, Φ(β), which is given by Eq. (17) (see Brown 2004), where β is the phase angle (angle at companion between star and the observer) and z = rsin⁽ν + ω⁾ is the radial coordinate of the radius vector. (17)The Jupiter/Sun contrast is obtained using Eq. (18), which gives an estimate of the fraction of stellar light captured by a planet, depending on the values of the planet radius, semi-major axis and geometrical albedo, being Φ(β) = 1 (at opposition) (18)where A_Jup = 0.35 is the value of the Jupiter albedo in the the J-band, (see Hanel et al. 1981).

We conclude with a final value of the contrast in reflected light given by Eq. (19). (19)As a consequence of Eq. (19), the results of MESS will be sensitive to the choice of A_λ, especially for the cold planets, in which the contribution of the reflected light is dominant. Following the outcomes of Jupiter observations and theoretical models (see e.g. Burrows 2004), we decided to uniformly generate the values of the albedo between 0.2 and 0.7. The code offers the option to fix the value of the albedo to a chosen value for all planets in the generation. A test of the impact of the choice of the albedo value on the redults of the simulations is presented in Sect. 3.4.

2.3.2. RV and astrometric signal

The indirect effects of the presence of the planet, such as the semi-amplitude of RV variations and the astrometric signal can be inferred, knowing all the orbital characteristics for each planet.

2.3.3. Degree of polarization

The degree of polarization Π is assumed to be of the form (see e.g. Stenflo 2005) (20)where Π_max is the maximum polarization value (which is assumed to be randomly generated between 0.1 and 0.3), and β is the same as in Eq. (17). Then the contrast due to the polarized light of the planets is Π times the contribution in reflected light evaluated with Eq. (19).

2.4. Planet population synthesis

Depending on the purpose of the analysis, the code can generate the planet population in two different ways:

• a)

  Full population: the value ofN_seed sets the spacing of the mass-period grid, and for each point on it N_gen planets are generated, ending with N_seed × N_gen planets per each star. The population for each star is saved in an independent file. This approach is useful for the statistical analysis of existing data, because in this case MESS provides the fraction of detectable planets per star, which can be used to derive the general probability of finding a planet over the whole target list. This can then be compared with the real results.

• b)

  Reduced population: only one orbit is generated for each point in the mass-period grid. N_gen in this case sets the number of planet in a planetary system associated with each star⁶. The final population is then composed by N_star × N_seed planets, and all planets are saved together in one file. Then the predicted detection performances of a given instrument can be used to derive the population of objects that are expected to be detected around each star if the whole input sample is observed.

As an example, we generated a reduced population (assuming five planets per star) of planets around the stars of the nearby sample described in Sect. 2.1. We chose the semi-empirical approach, and used the typical setup we discussed in Sect. 2.2.1⁷.

Figure 3 shows the position of the planets in the mass vs semi-major axis plane.

The planets are separated into the three classes, using different colours:

• Giant (or Jupiter-like) planets(M_planet > 40   M_Earth). A distinction between cold Jupiters (orange dots) and warm Jupiters (red dots), as defined in Sect. 2.3.1, is also made.

• Neptune-like planets (10   M_Earth ≤ M_planet ≤ 40   M_Earth: green dots).

• Rocky planets (M_planet < 10   M_Earth: blue dots).

The distribution of the observable quantities for the planet showed in Fig. 3 are summarized in Figs. 4 − 6.

2.5. Instrument detection performances

The last step is the comparison of the observables of the generated synthetic planets with the detection limits of different observing techniques with the possibility to actually combine them. In this framework, and in particular for comparative studies, it is important to make sure that the detection performances that are used as inputs for the code were estimated by correctly taking into account each instrumental biases, specific to each technique, and the stellar characteristics.

In the context of the MESS applications, the code has been extensively used considering two possible inputs for the detection performances:

• The 1D mode, which selects the detectable planets using athreshold or a curve giving the lower detection limits (RV,astrometric precisions or contrast performances) as a function ofthe period, the semi-major axis, the angular separations etc.,defined by the instrumental capabilities

• The 2D mode, which is especially built for the analysis of the performances of the deep imaging instruments. This mode takes advantage of the knowledge of all orbital elements of the planets and places them on a 2D detection map. This mode allows using the entire spatial information stored in the images. Using the whole 2D map not only allows us to take into account possible peculiar characteristics of the circumstellar environment, such as the presence of disks, but also prevent under/overestimation of the contrast curve that depends on the method chosen for the extraction itself (see Bonavita et al. 2010).

3. Applications

Once the synthetic population of planets has been created, the next step is to compare the characteristics of the generated planets with the detection limits appropriate for the instrument under consideration.

MESS offers three different operation modes (OM), depending on which kind of analysis is needed.

• 1.

  The hybrid mode (MESS_HM) which is the most flexible oneand can be used to probe the physical and orbital properties of aputative companion around one given system based on thecombination of different techniques, a priori information on thepossible orbit given the presence of other planets or circumstellardisk.

• 2.

  The statistical analysis mode (MESS_SAM), which is built for the analysis of real data and uses the full population defined in Sect. 2.4. It enables one to test different sets of planet populations or constrain the maximum occurrence of planets for a given population that would be consistent with the results of detection and/or null-detection of a complete survey of a large target sample.

• 3.

  The predictive mode (MESS_PM), which starts from the reduced population (see Sect. 2.4), and given the predicted performances of a planned instruments, can be used to select the most suitable targets given the science goals of the instrument itself, to test the results of different observing strategies and finally to predict possible synergies with other instruments.

3.1. Single object characterization

The first and more versatile MESS mode is the so-called hybrid mode. This mode can be used for the study of particularly interesting targets, or to test specific hypotheses. It allows for example to take into account all available informations about the orbit of a planet already discovered around the target, to put constraints on the planet generation. A preliminary version of this mode has been used to put constraints on the presence of a planetary companion embedded in the disk surrounding the T Tauri star LkCa15 (see Bonavita et al. 2010).

We present here an analogous analysis made for TWA 11. This star has been found to be surrounded by a debris disk by Schneider et al. (2009). Using STIS, Schneider et al. (2009) provided a full characterization of the disk geometry, and suggested a possible unseen companion responsible for some of the observed properties. We then decided to use MESS_HM to verify which constraints can be established using the VLT-NACO observations of this star.

A pixel-to-pixel 2D noise map was estimated from the reduced NACO images using a sliding box of 5 × 5 pixels over the whole field of view (FOV). We then considered a 6σ threshold to build the final detection limit maps to be used for the statistical analysis. These maps were also converted in terms of minimum mass map using the evolutionary model predictions at the age of the system. Figure 7 shows an example of the resulting sensitivity map⁸.

We considered only circular orbits coplanar with the disk, with an inclination and a longitude of the ascending node fixed by the disk properties reported in Schneider et al. (2009): i_Disk = 75.88 ± 0.16, Ω = PA ± 90 = (27.1 ± 90). Schneider et al. (2009) also estimated the inner and outer boundary of the disk to be 0.515 and 2.114 arcsec respectively, corresponding to 37 and 154 AU at the target distance (72.8 pc, see Van Leeuwen 2007). TWA11 is also known to have a stellar companion at ρ = 7.7′′ (Jura et al. 1995). As pointed out by Schneider et al. (2009), the value of the outer boundary of the disk is consistent with the presence of the companion. Using Eq. (1) we indeed obtained a value for the critical semi-major axis for the planet stability (a_crit) of about 165 AU. Taking into account these constraints, we set the range of explored semi-major axes to 35 − 160 AU. The results of our simulations in terms of non detection probability maps as a function of the companion mass and semi-major axes are shown in Fig. 8. The disk boundaries are also shown, as reported by Schneider et al. (2009).

Clearly with the NACO images we are not able to put strong constraints on planetary-mass objects, but surely low-star companions and brown dwarfs more massive than 30   M_Jup can be excluded at a > 35 AU and 20   M_Jup ones for a > 100 AU.

3.2. Statistical analysis of a survey

3.2.1. Testing the planet population assumption

The MESS_SAM operation mode allows one to test the consistency of various sets of (mass, eccentricity, semi-major axes) parametric distributions of a planet population with observational data. Given the detection performances of a survey, the frequency of detected simulated planets (over the complete sample) enables a derivation of the probability of non-detection of a given planet population associated with a normalized distribution set. Then the comparison with the survey results directly tests the disagreement with observations at an appropriate level of confidence.

As an example of the use of SAM/MESS statistical analysis mode, we present the analysis of a small sample of young neighbourhood stars that are part of the TWA association and were observed with NACO/VLT. These stars are part of a bigger sample for which the observations and statistical analysis performed with a preliminary version of MESS have been presented by Chauvin et al. (2010). The characteristics of the stars in the sample are listed in Table 1. We present a new analysis of these targets made with the 2D module of MESS_SAM.

2D minimum mass maps were obtained with the same method used for the analysis of TWA 11 (see Sect. 3.1) for all stars in the sample. We then used MESS_SAM to calculate the detection probability (P_D) of companions of various masses and orbital parameters (semi-major axis a, eccentricities e, inclination i, longitude of the ascending node Ω, longitude of periastron ω and time of periastron passage T_p). We used the empirical approach, generating a full population of 10 000 planets for each target, with a mass range spanning between 0.3 and 30 M_Jup and a cut-off in semi-major axis of 100 AU.

Each simulated companion was placed on the 2D minimum mass map according to its position on the projected orbit to test its detectability, comparing its mass with the minimum value achievable at the same position in the FoV.

Only circumbinary planets were considered around TWA 22, adopting the total mass of the system as M_Star. The binary is so close (ρ = 0.1′′ see Bonnefoy et al. 2009) the value of the critical semi-major axis for circumstellar planets a_CS is only 0.456 AU and of 8.395 AU for the circumbinary ones.

Two sets of indices for the power-law distribution were tested:

• 1.

  Those derived by Cumminget al. (2008, CM08):α =  −1.31, β =  −0.74.

• 2.

  Those derived by Lineweaver & Grether (2003, LW03): α =  −1.81, β =  −0.30.

Finally, fixing α and β to the CM08 values, we also introduced different values for the scaling of the planetary mass with the primary mass.

The results of these simulations are summarized in Fig. 9.

3.2.2. Estimate of the frequency of giant planets

A second more general use is to constrain the exoplanet fraction f within the physical separation and mass probed by a survey for null or positive detections. Contrary to what was assumed before, f becomes an output of the simulation, which actually depends on the assumed (mass, period, eccentricity) distributions of the giant planet population. This statistical analysis aims at determining f within a confidence range as a function of mass and semi-major axis, given a set of individual detection probabilities p_j directly linked to the detection limits of each star observed during the survey and to the considered giant planet distributions.

The probability of planet detection for a survey of N stars can be described by a binomial distribution, given a success probability fp_j, with f being the fraction of stars with planets, and p_j the individual detection probabilities of detecting a planet, if present around the star j. Each individual p_j can be replaced by  ⟨ p_j ⟩ , the mean survey detection probability of detecting a planet if present. Finally, assuming that the number of expected detected planets is small compared to the number of stars observed (f ⟨ p_j ⟩  ≪ 1), the binomial distribution can be approximated by a Poisson distribution to derive a simple analytical solution for the exoplanet fraction upper limit f_max for a given level of confidence CL: (21)Fig. 10 shows the results obtained applying this module at the sample of stars listed in Table 1.

Although the significance of our results is not really high, given the small size of the sample, they still agree with the results of the whole analysis presented by Chauvin et al. (2010) and with the results of the other deep imaging surveys (see e.g. Nielsen & Close 2010; Lafrenière et al. 2007)⁹.

3.2.3. Theoretical approach

The MESS_SAM can also be used to test the predictions of specific planet formation theories. An extensive use of this OM has been made to analyse a sample of massive stars (B-type and early A-type) observed with the Gemini Near-infrared Imager and spectroscoper (NIRI Hodapp et al. 2003), to test the applicability of planet formation by disk-instability in those systems. Starting from a uniform mass versus semi-major axis grid with a sampling of 5 AU in semi-major axis and 1   M_jup in mass, 10⁴ orbits were generated for each grid point. Models of disc instability (Bell et al. 1997; Mordasini et al. 2011; Klahr et al., in prep.) were then used to provide boundaries in the mass versus semi-major axis space, within which sub-stellar companions can form by this mechanism. These boundaries were dependent on the stellar properties, and so appropriate values should be used for each target in the sample. The planets falling within the allowed range were subsequently evaluated against the 1D detection limits from the high-contrast images of the survey. In this way, by testing a range of planet distributions within the set boundaries, meaningful limits could be placed on the frequency of planet and brown dwarf formation by disk-instability in massive disks. The full analysis is presented in detail in the survey paper (Janson et al. 2011). An example of a detection probability map in mass versus semi-major axis space is shown in Fig. 11.

3.3. Predictive mode

In addition to the analysis of the real data, MESS can also be used to predict the output of forthcoming searches, the goal being to provide information about the capabilities of future planet search instruments. With this mode, the flexibility of the code reaches its maximum, providing a wide range of possible applications.

Once the synthetic planet population has been created, and assuming the characteristics of a given instrument, MESS_PM allows predicting the number of detections expected from a future facility. This provides informations on

• 1.

  the expected frequency of planets;

• 2.

  the properties of these objects;

• 3.

  the kind of constraints that their observation can put on the planet formation theories.

Furthermore, it also allows one to test different instrumental configurations and observational strategies that can be adopted, thus providing a tool to tune the instrument characteristics to fulfil the requirements needed to access a certain domain in the parameter space, and reach the proposed science goals.

3.3.1. Comparison of future direct imaging instrument capabilities

As an example of the application of MESS_PM, we report the results of a comparison of the capability of a set of instruments for the direct imaging of exoplanets, planned for the next decade, which are briefly described in Table 2.

Because the purpose of the presented analysis is purely illustrative, we adopted for each instrument an averaged detectability relation, taken from the reference indicated in Table 2, and only the 1D approach is used. The sample of stars used is the one described in Sect. 2.2.1, whose properties are summarized in Fig. 1. This sample was originally selected as a preliminary sample for the planet search survey to be performed with with SPHERE, the next-generation planet finder of VLT (Beuzit et al. 2008), and it is therefore optimized for this kind of instrument, which possibly introduces some biases against some of the other instruments analysed. The analysis was made using the reduced population, assuming five planet per star.

The results of the analysis, shown in Fig. 12 and also summarized in Table 3, predict the enormous progress that can be expected in the next decade. The available measurements are already giving us indirect information on far-away planets around young stars, but passing through the intermediate step of next-generation image and finally with the advantage of ELT instruments we will have a wide view on planetary systems at different stages of their evolution.

3.3.2. Predicting the synergies between different techniques

Once the RV and astrometric modules will be completed, MESS_PM will provide an estimate of both the direct and indirect signatures of the presence of the planets, and can therefore be used to compare the outcomes of imaging with dynamical methods. These are interesting, because the latter allow determining the planet masses, thus eliminating the degeneracy with age, which is currently one of the major problems affecting direct detections. Moreover, possible synergies between different discovery methods are becoming more and more likely, ELT’s instruments representing the ideal link between direct and indirect detections, covering both young, nearby systems discovered by next generation imagers and also meant to provide the first images of planets already detected by RV.

Figures 13 and 14 summarize the results of the preliminary version of the RV and astrometric modules of MESS_PM. The planets shown are the same as in the lower right panel of Fig. 12.

If confirmed, these results would suggest that the discovery space for EPICS at E-ELT overlaps well with those from RV instruments (HARPS at ESO 3.6 m telescope, ESPRESSO at VLT, and especially CODEX at E-ELT) as well as with that of GAIA (Casertano et al. 2008).

Because the RV module is still under test, and we do not have enough data to perform a consistent and accurate analysis of the performances and comparison between the instruments under scrutiny, this analysis is not meant to tell which instrument is going to provide the highest number of detections, but just to show the potential of future versions of the code.

3.4. Testing the influence of the physical inputs

In this last section we present the results of some tests whose goal was to show how MESS can be used to investigate the influence of the various physical parameters considered as input for the planet generation. In particular we focused on the eccentricity distribution and on the value of the planetary albedo.

Eccentricity distribution

Direct-imaging surveys are, by definition, mostly sensitive to planets in wide orbits. Also, planets on highly eccentric orbits could be preferred targets, because they are more likely to be found farther out with respect to planets on a circular orbit with the same semi-major axis. This could lead to a bias towards high-eccentricity planets in our results. As mentioned in Sect. 2.2, the eccentricity distribution of the planets generated by MESS is uniform, and cuts at e = 0.6.

As an additional check of the impact of the eccentricity on the DI results, we repeated the analysis of the TWA sample performed in Sect. 3.2, by fixing the eccentricity for all generated planets to a given value. The results are shown in Fig. 15. The black solid line shows the results of the standard setup, for which the uniform eccentricity distribution cut at e = 0.6. The red, green, blue, purple, and light blue lines show the outcomes of the simulations made by fixing e = 0, e = 0.2, e = 0.4, e = 0.6 and e = 0.8, respectively. As expected, the higher eccentricity values can lead to an higher fraction of detected planets for a given semi-major axis value.

This simple exercise shows not only that the standard setup of the MESS does not introduce any systematic bias towards high eccentricity, but also that the code allows us to easily take this kind of bias into account, if they are proven to be real, by changing the simulation parameters.

As a final remark, the effect of the eccentricity is important only in the case of warm planets, such as those that could be found around our TWA targets. As the age of the stars increases, the reflected light contribution to the planet contrast becomes more and more important, which counterbalances the effect of the eccentricity.

Albedo distribution

As mentioned in Sect. 2.3, the albedo of the planets in the synthetic population is randomly generated between 0.2 and 0.7. Especially for the cold planets, the value of the albedo can be a critical parameter for the planet detection. We therefore decided to check how big the impact on the simulation results is. With an approach similar to the one used to test the eccentricity effect (see Sect. 3.3.3), we performed different sets of simulation A_λ being the only free parameter. We used a hypothetical G2V star (J = , age = 4.5 Gyrs) at 20 parsecs as target, and the detection limits of EPICS (see Table 2).

Figure 16 shows the results of the standard setup (A_λ randomly generated between 0.2 and 0.7, black solid line), together with those obtained by fixing A_λ to 0.2, 0.35 (the Jupiter value, see Sect. 2.3.1), 0.5 and 0.7 (red, green, blue, and purple line, respectively).

As expected, the fraction of detectable planets is higher for higher albedo, all other parameters being the same.

This test confirms that the use of a uniform distribution prevents us to favour planets with high albedo, and shows that the results are similar to the one obtained by using the Jupiter value.

4. Summary and conclusions

We presented MESS (multi-purpose exoplanet simulation system), a Monte Carlo tool for the statistical analysis and prediction of survey results for exoplanets.

Our aim was to build an extremely versatile code that could be used to test the outcomes of any instrument/technique for the detection of planets. We considered several assumptions on

• The star population, and how to take into account the properties ofeach star and their effect on either the characteristics of the planetsor the instrument capabilities. The binarity aspects are alsoincluded to take into account the possible effects of a stellarcompanion to the planet formation.

• The planet population, providing the complete set of orbital elements and a large number of physical parameters of the planets (radius, temperature, luminosity, etc.), either generated using the information from the analysis of the planets confirmed up to now (semi-empirical approach) or using the results of the planet formation theories (theoretical approach).

• The predicted observables (luminosity and polarimetric contrast, RV semi-amplitude, astrometric signal).

• The synthesis of a planet population that can be easily adapted to the purpose of the investigation.

• The final comparison with the detection limits, with the possibility to combine the informations from different observing techniques, to select a sub-sample of detectable planets whose characteristics can then be investigated.

The code is such that each and every one of these assumptions can be released and/or changed. This not only provides a tool that is independent of the models (e.g. the planet-formation theory chosen if the theoretical approach is used, or the evolutionary models used to estimate the planet luminosity and radius) but also makes it relevant to test model prediction, and to constraint the properties of the known planets under different initial conditions.

So far only the direct-imaging module of the code has been extensively used, but the combination of various techniques is under test and will offer rich perspective for future combined studies of exoplanets.

Three main applications of the MESS code have been shown:

• 1.

  The hybrid mode, built for the analysis of single objects, waspresented in Sect. 3.1. It can be used to probe thephysical and orbital properties of a putative companion around agiven system based on the combination of different techniques,and possibly a priori information on the orbit given the presence ofother planets or of a circumstellar disk.

• 2.

  The SAM mode (Sect. 3.2), optimized for the analysis of a large sample of stars, showed its full potential in Sect. 3.2 by providing a detailed statistical analysis of a sample of stars observed with direct imaging. Both the agreement of the observations with the observed parameter distributions (Sect. 3.2.1) and the planet-formation theories (Sect. 3.2.3) were tested, using the semi-empirical and theoretical approach, respectively.

• 3.

  The PM mode finally aims at the prediction of the outcomes of future searches, and can be used to tune not only the main instrument parameters, but even the observing strategy.

However, an extensive use of the code requires a complete knowledge of the instrument under test, of all error sources and of the detection capabilities. Then, to really extend the use of MESS to other facilities one should first properly set all needed parameters. As already mentioned before, both the RV and the astrometric part are currently included in a very simplistic way. A better treatment of the dependence of the detectability with astrometry from the orbital parameters should be included. A rigorous treatment of the stellar jitter evaluation must be implemented to allow a better comparison between the imaging and radial velocity capabilities. Especially in the case of E-ELT instruments, this would allow one to better define the synergies between the various channels, for a more focused observing strategy.

Moreover, a precise measure of the stellar characteristics is also needed, to minimize the effects that errors on parameters such as the age of the system or the presence of stellar companions can have on the analysis.

Finally, the inclusion of an analysis of the planet stability in case of multiple objects is planned, together with an extensive use of the theoretical approach, using the outcomes of the most recent Bern models (Mordasini et al. 2010).

Each technique performance varies with the star properties (age, mass, distance...), have different observables (luminosity, minimum mass, radius...), different observing strategies. It is therefore extremely important to take this into account to achieve the maximum constrains, first on the properties of giant planets (physical, orbital parameter space) that will entirely shape the planetary system architecture, then possibly on the telluric planets. A better characterization of the giant and telluric planet orbital and physical properties, including their dependency with the host properties, is critical for a better understanding of their formation processes because various mechanisms may be at play (Boley 2009; Dodson-Robinson et al. 2009; Stamatellos & Whitworth 2009), but also for their architecture and dynamical evolution. In the end, one additional and important point is to understand the required physical conditions that will lead to the formation of telluric planets in the habitable zone within planetary systems shaped by giant planets, and which will possibly lead to the development of life.

---

Acknowledgments

This paper was part of M. Bonavita’s Ph.D. thesis, which was funded by the Italian Institute for Astrophysics (INAF).

The authors would like to thank Prof. W. Benz and Dr. Y. Alibert for providing essential input for the development of the theoretical approach of the MESS.

We also acknowledge support from the French National Research Agency (ANR) through project grant ANR10-BLANC0504-01.

References

All Tables

All Figures