# MESS (multi-purpose exoplanet simulation system). A Monte Carlo tool for the statistical analysis and prediction of exoplanet search results

M. Bonavita, G. Chauvin, S. Desidera, R. Gratton, M. Janson, J. L. Beuzit, M. Kasper, C. Mordasini

<figure><img src="image_1.png" /></figure>

To cite this version:

M. Bonavita, G. Chauvin, S. Desidera, R. Gratton, M. Janson, et al., MESS (multi-purpose exoplanet simulation system). A Monte Carlo tool for the statistical analysis and prediction of exoplanet search results. Astronomy and Astrophysics - A&A, 2012, 537, 10.1051/0004-6361/201116852. insu-03612457

HAL Id: insu-03612457

https://insu.hal.science/insu-03612457

Submitted on 17 Mar 2022

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# MESS (multi-purpose exoplanet simulation system)*,**

## A Monte Carlo tool for the statistical analysis and prediction of exoplanet search results

M. Bonavita<fnref n="1" />,<fnref n="2" />, G. Chauvin<fnref n="3" />,<fnref n="4" />, S. Desidera<fnref n="1" />, R. Gratton<fnref n="1" />, M. Janson<fnref n="2" />, J. L. Beuzit<fnref n="3" />, M. Kasper<fnref n="5" />, and C. Mordasini<fnref n="4" />

Received 8 March 2011 / Accepted 21 October 2011

## ABSTRACT

*Context.* The high number of planet discoveries made in the last years provides a good sample for statistical analysis, leading to some clues on the distributions of planet parameters, such as masses and periods, at least in close proximity to the host star. We likely need to wait for the extremely large telescopes (ELTs) to have an overall view of the extrasolar planetary systems. Those facilities will finally ensure an overlap of the discovery space of direct and indirect techniques, which is desirable to completely understand the nature of the discovered objects, obtaining both orbital parameters and physical characterization.

*Aims.* In this context it would be useful to have a tool that can be used for the interpretation of the present results that are obtained with various observing techniques, and also to predict what the outcomes would be of the future instruments.

*Methods.* For this reason we built MESS: a Monte Carlo simulation code that uses either the results of the statistical analysis of the properties of discovered planets or the results of the planet formation theories to build synthetic planet populations that are fully described in terms of frequency, orbital elements and physical properties. These synthetic planets can then be used to either test the consistency of their properties with the observed ones given different detection techniques (radial velocity, imaging and astrometry) or to predict the expected number of planets for future surveys, as well as to optimize the future multi-technique observations for their characterization down to telluric masses.

*Results.* In addition to the code description, we present here some of its applications to probe the physical and orbital properties of a putative companion within the circumstellar disk of a given star and to test or constrain the orbital distribution properties of a potential planet population around the members of the TW Hydrae association. Finally, we investigated the synergy of future space and ground-based telescopes instrumentation with the predictive mode of the code, to identify the mass-period parameter space that will be probed in future surveys for giant and rocky planets.

**Key words.** brown dwarfs – methods: statistical – methods: data analysis – planetary systems

## 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.

<footnote n="1">¹ INAF – Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, 35122 Padova, Italy</footnote>

<footnote n="2">² Department of Astronomy and Astrophysics, University of Toronto, 50 St. George Street M5S 3H4 Toronto ON, Canada</footnote>

<footnote n="3">³ UJF-Grenoble 1/CNRS-INSU, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG) UMR 5274, 38041 Grenoble, France</footnote>

<footnote n="4">⁴ Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="5">⁵ European Southern Observatory (ESO), Karl-Schwarzschild-Str. 2, 85748 Garching, Germany</footnote>

<footnote>* www.messthecode.com</footnote>

<footnote>** The code is only available at CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/537/A67</footnote>
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 IDL®<fnref n="1" />, 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 all orbital elements, either using the planet mass and period distributions coming the statistical results of RV studies or the outcome of the planetary formation theories;

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Principal characteristics of the sample of nearby stars used to build the example synthetic population. *Upper left:* apparent magnitude in the J band vs. distance in pc. *Upper right:* stellar mass ($M_\odot$) vs. distance in pc. *Lower left:* histogram of stellar ages (Myrs). *Lower right:* histogram of stellar masses ($M_\odot$).</figcaption></figure>

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_{\text{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).

<footnote n="1">¹ http://www.ittvis.com</footnote>
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),

$$a_c/a_b = 0.464 - 0.38 \mu - 0.361 e_b + 0.586 \mu e_b + 0.150 \mu^2 - 0.198 \mu e_b^2 \quad (1)$$

$$a_c/a_b = 1.60 + 4.12 \mu + 5.10 e_b - 4.27 \mu e_b - 5.09 \mu^2 - 2.22 e_b^2 + 4.61 e_b^2 \mu^2. \quad (2)$$

In both equations, $a_c$ is critical semi-major axis, $\mu = M_1/(M_1 + M_2)$, $a_b$ and $e_b$ are the semi-major axis and eccentricity of the binary, and $M_1$ and $M_2$ 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 \rho$ (arcsec) d (pc) (see Fischer et al. 2002; Duquennoy & Mayor 1991)$^2$.

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$^3$.

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_{\text{seed}}$ values of masses and periods (see Sect. 2.2.1)

$$\frac{dN}{d(M_{\text{p}})} \propto (M_{\text{p}})^{\alpha} \quad (3)$$

$$\frac{dN}{dP} \propto P^{\beta}. \quad (4)$$

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 $\alpha = -1.31$ and $\beta = -0.74$ respectively, according to Cumming et al. (2008). The planetary masses span the range $0.6 M_{\text{Earth}}$ and $15 M_{\text{Jup}}$, and the periods (P) are chosen to be between 2.5 days and 350 years (corresponding to 50 AU for $1 M_{\odot}$ 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_{\text{gen}}$ values of all orbital parameters: eccentricity ($e$), inclination ($i$), longitude of periastron ($\omega$), longitude of ascending node ($\Omega$), and time of periastron passage ($T_0$). By default, all these parameters are randomly generated following uniform distribution$^4$. 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_{\text{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).

$$x = AX + FY \quad (5)$$

$$y = BX + GY \quad (6)$$

$$X = \cos E - e \quad (6)$$

$$Y = \sqrt{1 - e^2} \sin E \quad (7)$$

$$\rho = \sqrt{x^2 + y^2}, \quad (7)$$

where $X$ and $Y$ are the coordinates of the orbit (Eq. (6)), $\rho$ is the projected separation, and $A, B, F, G$ are the Thiele-Innes

<footnote>$^2$ Note that Duquennoy & Mayor (1991) refer to solar-type star multiplicity.</footnote>

<footnote>$^3$ The results discussed in this paper were obtained using mainly the semi-empirical approach. An extensive use of the theoretical approach, using as input the newest outcomes of the Bern formation models (Mordasini et al. 2011), will be the subject of a forthcoming paper.</footnote>

<footnote>$^4$ Note that in the case of the inclination, $\cos i$ and not $i$ itself is uniformly generated by the code.</footnote>
elements, which can be obtained from the classical ones (the semi-major axis *a*, $\omega$, $\Omega$, and *i*) using Eq. (8):

$$
\begin{align}
A &= a(\cos \omega \cos \Omega - \sin \omega \sin \Omega \cos i) \nonumber \\
B &= a(\cos \omega \sin \Omega + \sin \omega \cos \Omega \cos i) \tag{8} \\
F &= a(-\sin \omega \cos \Omega - \cos \omega \sin \Omega \cos i) \nonumber \\
G &= a(-\sin \omega \sin \Omega + \cos \omega \cos \Omega \cos i). \nonumber
\end{align}
$$

In these equations, *E* is the eccentric anomaly (obtained from the mean anomaly *M* (Eq. (9)) using Eq. (10)) and *v* the true anomaly (Eq. (11)):

$$
\begin{align}
M &= \left( \frac{t_{\text{obs}} - T_0}{p} \right) 2\pi \tag{9} \\
E_0 &= M + e \sin M + \frac{e^2}{2} \sin 2M \nonumber \\
M_0 &= E_0 - e \sin E_0 \nonumber \\
E &= \frac{E_0 + (M - M_0)}{1 - e \cos E_0} \tag{10} \\
\tan v/2 &= \sqrt{\frac{(1+e)}{1-e}} \tan E/2. \tag{11}
\end{align}
$$

The projected separation, $\rho$ (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,

$$
\begin{align}
\rho &= r \cos(v + \omega) \sec(\theta - \Omega) \tag{12} \\
r &= a(1 - e^2)/(1 + e \cos v). \nonumber
\end{align}
$$

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 tem-
perature. The first one is the internal temperature, Tint, obtained
from the evolutionary models (see e.g. Baraffe et al. 2003). The
second one is the equilibrium temperature, Teq, obtained through
Eq. (13) (from Sudarsky et al. 2003)

$$
T_{\text{eq}} = \left[ \frac{(1 - A_{\text{B}}) L_{*}}{16\pi\sigma a^2} \right], \qquad (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_{\text{eff}}$ is given by

$$
T_{\text{eff}}^4 = T_{\text{int}}^4 + T_{\text{eq}}^4. \qquad (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 \ge 100 M_{\text{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_{\text{Gas}}$. A core mass of 10 $M_{\text{Earth}}$ is assumed.

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Summary of the planetary mass-radius relations adopted for the different mass ranges. All model computation are made assuming a host star of 1 $M_{\odot}$, and the semi-major axis value is fixed to 5 AU. Filled symbols correspond to known transiting planets; open symbols are for solar system planets.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Mass semi-major axis distribution of the synthetic planets in the populations generated by MESS using the semi-empirical approach. The different classes of planets (see text) are plotted using different colours: red/orange for the warm/cold Jupiters, green for the Neptune-like planets, blue for the rocky planets.</figcaption></figure>

2. Equations (15) and (16) from Fortney et al. (2007) are used for the smallest planets ($M \le 10 M_{\text{Earth}}$). These are either

$$
\begin{equation}
\begin{split}
R ={}& (0.0912 \quad \mathit{imf} + 0.1603)(\log M)^2 \\
      & + (0.3330 \quad \mathit{imf} + 0.7387) \log M \\
      & + (0.4639 \quad \mathit{imf} + 1.1193)
\end{split}
\tag{15}
\end{equation}
$$

or

$$
\begin{equation}
\begin{split}
R ={}& (0.0592 \quad rmf + 0.0975)(\log M)^2 \\
      & + (0.2337 \quad rmf + 0.4938) \log M \\
      & + (0.3102 \quad rmf + 0.7932)
\end{split}
\tag{16}
\end{equation}
$$

for ice/rock and rock/iron planets, respectively. In these
equations, *R* is in *R*<sub>Earth</sub> and *M* is in *M*<sub>Earth</sub>, 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 plan-
ets, 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_{\text{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<fnref n="5" /> 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, $\Phi(\beta)$, which is given by Eq. (17) (see Brown 2004), where $\beta$ is the phase angle (angle at companion between star and the observer) and $z = r \sin(\nu + \omega)$ is the radial coordinate of the radius vector.

$$ \Phi(\beta) = [\sin\beta + (\pi + \beta) \cos\beta] / \pi, \quad (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 $\Phi(\beta) = 1$ (at opposition)

$$ (L_{\text{Jup}}/L_{*})_{\text{Ref}} = A_{\text{Jup}} \frac{R_{\text{Jup}}^2}{a_{\text{Jup}}^2} = 2.5 \times 10^{-9}, \quad (18) $$

where $A_{\text{Jup}} = 0.35$ is the value of the Jupiter albedo in the J-band, (see Hanel et al. 1981).

We conclude with a final value of the contrast in reflected light given by Eq. (19).

$$ (L_p/L_*)_{\text{Ref}} = (L_{\text{Jup}}/L_*)_{\text{Ref}} \Phi(\beta) \frac{(R_p/R_{\text{Jup}})^2}{(a/a_{\text{Jup}})^2}. \quad (19) $$

As a consequence of Eq. (19), the results of MESS will be sensitive to the choice of $A_\lambda$, 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 $\Pi$ is assumed to be of the form (see e.g. Stenflo 2005)

$$ \Pi = \Pi_{\max} \times (1 - \cos^2\beta)/(1 + \cos^2\beta), \quad (20) $$

where $\Pi_{\max}$ is the maximum polarization value (which is assumed to be randomly generated between 0.1 and 0.3), and $\beta$ is the same as in Eq. (17). Then the contrast due to the polarized light of the planets is $\Pi$ 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 of $N_{\text{seed}}$ sets the spacing of the mass-period grid, and for each point on it $N_{\text{gen}}$ planets are generated, ending with $N_{\text{seed}} \times N_{\text{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_{\text{gen}}$ in this case sets the number of planet in a planetary system associated with each star<fnref n="6" />.

<footnote n="6">⁶ Note that no consideration on the planet stability is made, and for the purpose of the analysis each planet is considered separately.</footnote>

<footnote n="5">⁵ www.exoplanet.eu</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Distribution of radial velocity vs. period of the synthetic planets for the population shown in Fig. 3.</figcaption></figure>

The final population is then composed by $N_{\text{star}} \times N_{\text{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[^7].

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_{\text{planet}} > 40 M_{\text{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_{\text{Earth}} \le M_{\text{planet}} \le 40 M_{\text{Earth}}$: green dots).

- Rocky planets ($M_{\text{planet}} < 10 M_{\text{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 a threshold or a curve giving the lower detection limits (RV, astrometric precisions or contrast performances) as a function

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Distribution of the astrometric signal vs. period of the synthetic planets for the population shown in Fig. 3.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Planet/star contrast vs. projected separation of the same planets shown in Fig. 3.</figcaption></figure>

of the 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 one and can be used to probe the physical and orbital properties of a putative companion around one given system based on the combination of different techniques, a priori information

<footnote>[^7]: Note that the whole calculation of the physical characteristics and observables described in Sect. 2.2.3 to 2.3.2 can be skipped (with considerable gain in computing speed), the code then provides only the orbital elements.</footnote>
on the possible orbit given the presence of other planets or circumstellar disk.

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<fnref n="8" />.

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_{\text{Disk}} = 75.88 \pm 0.16$, $\Omega = \text{PA} \pm 90 = (27.1 \pm 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 $\rho = 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_{\text{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

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** 2D map giving the values of the minimum mass of detectable companions (6σ) as a function of projected separation around TWA 11.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Non detection probability map of a faint companion around TWA11 as a function of its mass and semi-major axis, in the case of a circular orbit. Inclination and longitude of the ascending node have been fixed using the disk properties $i_{\text{Disk}} = 75.88 \pm 0.16$, $\Omega = \text{PA} \pm 90 = (27.1 \pm 90)$.</figcaption></figure>

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_{\text{Jup}}$ can be excluded at $a > 35$ AU and 20 $M_{\text{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.

<footnote n="8">⁸ Note that the decreasing values of the non-detection probability at separations lower than 30–50 AU are caused by systematic errors. In fact the detection limit drops to unrealistic low values really close to the star.</footnote>
<table><caption>**Table 1.** Sample of TWA stars considered in our analysis.</caption><thead><tr><th>Name</th><th>α<br>[J2000]</th><th>δ<br>[J2000]</th><th>SpT</th><th>Mass<br>M<sub>☉</sub></th><th>d<br>(pc)</th><th>Age<br>(Myr)</th><th>V<br>(mag)</th><th>K<br>(mag)</th><th>Mode</th><th>Notes</th></tr></thead><tbody><tr><td>TWA22</td><td>10 17 26.9</td><td>-53 54 28</td><td>M5</td><td>0.15</td><td>18</td><td>8</td><td>13.2</td><td>7.69</td><td>DI/CI, S27, Ks</td><td>Bin (ρ = 0.1")</td></tr><tr><td>TWA14</td><td>11 13 26.3</td><td>-45 23 43</td><td>M0</td><td>0.55</td><td>63</td><td>8</td><td>13.8</td><td>8.50</td><td>DI/CI, S27, Ks</td><td></td></tr><tr><td>TWA12</td><td>11 21 05.6</td><td>-38 45 16</td><td>M2</td><td>0.30</td><td>32</td><td>8</td><td>13.6</td><td>8.05</td><td>DI/CI, S27, Ks</td><td></td></tr><tr><td>TWA19</td><td>11 47 24.6</td><td>-49 53 03</td><td>G5</td><td>1.50</td><td>104</td><td>8</td><td>9.1</td><td>7.51</td><td>DI/CI, S13, H</td><td></td></tr><tr><td>TWA23</td><td>12 07 27.4</td><td>-32 47 0</td><td>M1</td><td>0.40</td><td>37</td><td>8</td><td>12.7</td><td>7.75</td><td>DI/CI, S13, H</td><td></td></tr><tr><td>Twa25</td><td>12 15 30.7</td><td>-39 48 42</td><td>M5</td><td>0.15</td><td>44</td><td>8</td><td>11.4</td><td>7.31</td><td>DI/CI, S27, Ks</td><td></td></tr><tr><td>TWA11</td><td>12 36 01.0</td><td>-39 52 10</td><td>A0</td><td>2.10</td><td>67</td><td>8</td><td>5.8</td><td>5.77</td><td>DI/CI, S27, H</td><td>Bin (ρ = 7.7"), Star with disk</td></tr><tr><td>Twa17</td><td>13 20 45.4</td><td>-46 11 38</td><td>K5</td><td>1.00</td><td>133</td><td>8</td><td>12.6</td><td>9.01</td><td>DI/CI, S13, H</td><td></td></tr></tbody></table>

**Notes.** In addition to name, coordinates, galactic latitude (b), spectral type, distance, V and K photometry, we also list the observing mode direct imaging (DI) or coronagraphy (COR) and the status of the primary (single, binary Bin, triple).

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 Cumming et 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

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Non-detection probability for the stars listed in Table 1 based on various sets of period and mass distribution. Mass and period distribution are extrapolated and normalized from RV studies. Top: variation of the non-detection probability using two different sets of power-law distributions (see text). Bottom: variation of the non-detection probability fixing $α = -1.31$ and $β = -0.74$ (Cumming et al. 2008) and different scaling the planet mass with the primary mass.</figcaption></figure>

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 $\langle p_j \rangle$, 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\langle p_j \rangle \ll 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:

$$f_{max} = \frac{-\ln(1 - \text{CL})}{N\langle p_j \rangle}, \quad (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
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Top: survey mean detection-probability derived as a function of the semi-major axis assuming the parametric mass and period distributions derived by Cumming et al. (2008). The results are reported for individual masses: 1.5, 3, 6, 9, 12, 15, 30 $M_{\text{Jup}}$. The integrated probability for the planetary mass regime is shown with the thick green line. Bottom: planet fraction upper limit derived as a function of semi-major axis given the same mass and period distributions.</figcaption></figure>

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)<fnref n="9" />.

### 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_{\text{Jup}}$ in mass, $10^4$ 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

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Non-detection probability map for 41 Ari (HIP 13209).</figcaption></figure>

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 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

<footnote n="9">⁹ An extensive analysis with MESS_SAM of the results of the major deep-imaging surveys published in the last decade is ongoing, and will be presented in a forthcoming dedicated paper.</footnote>
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Planets expected to be discovered by SPHERE (representative of planet finders on 8 m class ground-based telescopes), JWST-MIRI, 1.5 m space coronagraphs, and EPICS/E-ELT (representative of 30–40 m class telescopes) in the mass vs separation plane. Different colours are used for warm giants (orange), cold giants (red), Neptune-like (green), and rocky planets (blue), respectively.</figcaption></figure>

<table><caption>**Table 2.** Instruments for direct imaging of exoplanets considered in our analysis.</caption><thead><tr><th>Instrument</th><th>Contrast</th><th>Wavel.<br>(µm)</th><th>IWA<br>(")</th><th>Year</th><th>Ref.</th></tr></thead><tbody><tr><td colspan="6">8 m ground-based telescopes</td></tr><tr><td>VLT-SPHERE<br>Gemini-GPI</td><td>10<sup>-7</sup></td><td>0.9–1.7</td><td>0.08</td><td>2011</td><td>B10<br>G07</td></tr><tr><td colspan="6">JWST</td></tr><tr><td>NIRCAM</td><td>10<sup>-5</sup></td><td>2.1–4.6</td><td>0.30</td><td rowspan="2">2014</td><td>G05</td></tr><tr><td>MIRI</td><td>10<sup>-4</sup></td><td>5–25</td><td>0.35</td><td>R10</td></tr><tr><td colspan="6">1.5 m space coronagraphs</td></tr><tr><td colspan="6">10<sup>-9</sup>–10<sup>-10</sup><br>ELT class instruments</td></tr><tr><td>E-ELT-EPICS</td><td>10<sup>-8</sup>–10<sup>-9</sup></td><td>0.9–1.7</td><td>0.03</td><td rowspan="2">&gt;2018</td><td>K10</td></tr><tr><td>E-ELT-METIS</td><td>10<sup>-5</sup></td><td>2.5–20</td><td>0.08</td><td>S10</td></tr></tbody></table>

**References.** B10 Beuzit et al. (2010); G05 Green et al. (2005); G07 Graham et al. (2007); K10 Kasper et al. (2010); R10 Rieke et al. (2010); S10 Stuik et al. (2010); T10 Trauger & Moody (2010).

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

<table><caption>**Table 3.** Summary of expected detections from imagers in the next decade.</caption><thead><tr><th>Instrument</th><th>Year</th><th>Young<br>Giants</th><th>Old<br>Giants</th><th>Nept.</th><th>Rocky</th></tr></thead><tbody><tr><td>Gr. based 8 m</td><td>2011</td><td>tens</td><td>few</td><td></td><td></td></tr><tr><td>JWST</td><td>2014</td><td>tens</td><td>few</td><td></td><td></td></tr><tr><td>1.5 m Space Coro.</td><td>?</td><td>tens</td><td>tens</td><td>tens</td><td>few</td></tr><tr><td>ELT's</td><td>&gt;2018</td><td>hundr.</td><td>hundr.</td><td>tens</td><td>few</td></tr></tbody></table>

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).
<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Planets expected to be detected by EPICS (nearby sample) in the RV signal vs. period plane, compared with detection limits for RV instruments (HARPS, ESPRESSO and CODEX). The colour code is the same as in Fig. 1.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Planets expected to be detected by EPICS (nearby sample) in the astrometric signal vs. period plane, compared with detection limits for astrometric satellites GAIA. Colour code is the same as in Fig. 1.</figcaption></figure>

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

<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Fraction of detected planet as a function of the semi-major axis value (AU) for different values of the planet eccentricity. The black solid line shows the results obtained with the *standard setup*.</figcaption></figure>

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 a 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_\lambda$ being the only free parameter. We used a hypothetical G2V star ($J = 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_\lambda$ randomly generated between 0.2 and 0.7, black solid line), together with those obtained by fixing $A_\lambda$ 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).
<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Fraction of detected planet as a function of the semi-major axis value (AU) for different values of the planet albedo ($A_{\lambda}$). The black solid line shows the results obtained with the *standard setup*.</figcaption></figure>

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 of each star and their effect on either the characteristics of the planets or the instrument capabilities. The binarity aspects are also included to take into account the possible effects of a stellar companion 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 constrain 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, was presented in Sect. 3.1. It can be used to probe the physical and orbital properties of a putative companion around a given system based on the combination of different techniques, and possibly a priori information on the orbit given the presence of other 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.
**Acknowledgements.** 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

Absil, O., & Mawet, D. 2009, A&Ar, 16

Baraffe, I., Chabrier, G., Barman, T. S., Allard, F., & Hauschildt, P. H. 2003, A&A, 402, 701

Bell, K. R., Cassen, P. M., Klahr, H. H., & Henning, T. 1997, ApJ, 486, 372

Beuzit, J., Feldt, M., Dohlen, K., et al. 2008, in Ground-based and Airborne Instrumentation for Astronomy II. ed. I. S. McLean, & M. M. Casali, Proc. SPIE, 7014, 701418

Beuzit, J., Boccaletti, A., Feldt, M., et al. 2010, in Pathways Towards Habitable Planets, ed. V. Coudé Du Foresto, D. M. Gelino, & I. Ribas, ASP Conf. Ser., 430, 231

Boley, A. C. 2009, ApJ, 695, L53

Bonavita, M., Chauvin, G., Boccaletti, A., et al. 2010, A&A, 522, A2

Bonnefoy, M., Chauvin, G., Dumas, C., et al. 2009, A&A, 506, 799

Brown, R. A. 2004, ApJ, 610, 1079

Burrows, A. 2004, in Planetary Systems in the Universe, ed. A. Penny, IAU Symp., 202, 255

Casertano, S., Lattanzi, M. G., Sozzetti, A., et al. 2008, A&A, 482, 699

Chabrier, G., Baraffe, I., Allard, F., & Hauschildt, P. 2000, ApJ, 542, 464

Chauvin, G., Lagrange, A., Bonavita, M., et al. 2010, A&A, 509, A52

Cumming, A., Butler, R. P., Marcy, G. W., et al. 2008, PASP, 120, 531

Dodson-Robinson, S. E., Veras, D., Ford, E. B., & Beichman, C. A. 2009, ApJ, 707, 79

Duquennoy, A., & Mayor, M. 1991, A&A, 248, 485

Fischer, D. A., & Valenti, J. 2005, ApJ, 622, 1102

Fischer, D. A., Marcy, G. W., Butler, R. P., et al. 2002, PASP, 114, 529

Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, ApJ, 659, 1661

Fortney, J. J., Marley, M. S., Saumon, D., & Lodders, K. 2008, ApJ, 683, 1104

Graham, J. R., Macintosh, B., Doyon, R., et al. 2007, in BAAS 38, Am. Astron. Soc. Meet. Abstracts, 134.02

Green, J. J., Beichman, C., Basinger, S. A., et al. 2005, in Techniques and Instrumentation for Detection of Exoplanets II, ed. D. R. Coulter, Proc. SPIE, 5905, 185

Hanel, R., Conrath, B., Herath, L., Kunde, V., & Pirraglia, J. 1981, J. Geophys. Res., 86, 8705

Heintz, W. D. 1978, Geophys. Astrophys. Monographs, 15

Hodapp, K. W., Jensen, J. B., Irwin, E. M., et al. 2003, PASP, 115, 1388

Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621

Janson, M., Bonavita, M., Klahr, H., et al. 2011, ApJ, 736, 89

Johnson, J. A., Butler, R. P., Marcy, G. W., et al. 2007, ApJ, 670, 833

Jura, M., Ghez, A. M., White, R. J., et al. 1995, ApJ, 445, 451

Kalas, P., Graham, J. R., Chiang, E., et al. 2008, Science, 322, 1345

Kasper, M., Apai, D., Janson, M., & Brandner, W. 2007, A&A, 472, 321

Kasper, M., Beuzit, J., Verinaud, C., et al. 2010, in Ground-based and Airborne Instrumentation for Astronomy III ed. I. S. McLean, S. K. Ramsay, & H. Takami, Proc. SPIE, 7735, 77352E

Lafrenière, D., Doyon, R., Marois, C., et al. 2007, ApJ, 670, 1367

Lafrenière, D., Jayawardhana, R., & van Kerkwijk, M. H. 2008, ApJ, 689, L153

Lagrange, A., Gratadour, D., Chauvin, G., et al. 2009, A&A, 493, L21

Lineweaver, C. H., & Grether, D. 2003, ApJ, 598, 1350

Lovis, C., & Mayor, M. 2007, A&A, 472, 657

Macintosh, B., Graham, J. R., Palmer, D., et al. 2007, in BAAS 38, 782

Marley, M. S., Fortney, J. J., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2007, ApJ, 655, 541

Marois, C., Macintosh, B., Barman, T., et al. 2008, Science, 322, 1348

Mordasini, C., Alibert, Y., & Benz, W. 2009, A&A, 501, 1139

Mordasini, C., Klahr, H., Alibert, Y., Benz, W., & Dittkrist, K. 2010 [arXiv:1012.5281]

Mordasini, C., Alibert, Y., Benz, W., & Klahr, H. 2011 [arXiv:1102.4146]

Nielsen, E. L., & Close, L. M. 2010, ApJ, 717, 878

Perryman, M. A. C., & ESA, eds. 1997, The HIPPARCOS and TYCHO catalogues. Astrometric and photometric star catalogues derived from the ESA HIPPARCOS Space Astrometry Mission, ESA, 1200

Rieke, G., Wright, G., Glasse, A., Ressler, M., & MIRI Science Team 2010, in BAAS, 42, Am. Astron. Soc. Meet. Abstracts, 215, 439.04

Santos, N. C., Israelian, G., & Mayor, M. 2004, A&A, 415, 1153

Saumon, D., & Marley, M. S. 2008, ApJ, 689, 1327

Schneider, G., Weinberger, A. J., Becklin, E. E., Debes, J. H., & Smith, B. A. 2009, AJ, 137, 53

Stamatellos, D., & Whitworth, A. P. 2009, MNRAS, 400, 1563

Stenflo, J. O. 2005, A&A, 429, 713

Stuik, R., Jolissaint, L., Kendrew, S., et al. 2010, in Adaptive Optics Systems II, ed. B. L. Ellerbroek, M. Hart, N. Hubin, & P. L. Wizinowich, Proc. SPIE, 7736, 77363

Sudarsky, D., Burrows, A., & Hubeny, I. 2003, ApJ, 588, 1121

Trauger, J. T., & Moody, D. C. 2010, in Am. Astron. Soc. Meet. Abstracts, 216, 425.01

Van Leeuwen, F. 2007, A&A, 474, 653