ASTROMETRIC EXOPLANET DETECTION WITH GAIA

Gaia, as for its predecessor Hipparcos, utilizes a small number of key measurement principles (observations above the atmosphere, two widely separated viewing directions, and a uniform "revolving scanning" of the celestial sphere) to create a catalog of star positions, proper motions, and parallaxes of state-of-the-art accuracies (Perryman et al. 2001). Crucially, both missions provide absolute trigonometric parallaxes, rather than the relative parallaxes accessible to narrow-field astrometry from the ground. The observations are reduced to an internally consistent and extremely "rigid" catalog of positions and proper motions, but whose system orientation and angular rate of change are essentially arbitrary, since the measured arc lengths between objects are invariant to frame rotation. Placing both positions and proper motions on an inertial system corresponds to determining these six degrees of freedom (three orientation and three spin components). For Gaia, they will be derived using the large numbers of observed quasars (Claeskens et al. 2006; Perryman et al. 2014).

On board detection ensures that objects brighter than G ∼ 20 at that measurement epoch will be detected and observed astrometrically and photometrically, the latter through low-resolution spectrophotometry at the trailing edge of the astrometric field (see Jordi et al. 2010, Figures 1–2). The highest photometric accuracy will come from the unfiltered G-band astrometric field photometry, which will range (per field crossing) from 1 mmag or better for G < 14 mag to ∼0.2 mag at G = 20 mag (Jordi et al. 2010, Figure 19).

Final astrometric accuracies (in positions, parallaxes, and annual proper motions) should be roughly constant at ∼10 μas (micro-arcsec) between V ∼ 7–12, degrading according to photon statistics to ∼20–25 μas at V = 15, and to ∼300 μas at V = 20 (precise values depend on photometric passband, star color, and astrometric parameter). These final accuracies result from the combination of the one-dimensional measurements throughout the mission, assembled using a global iterative adjustment (Perryman et al. 2001; O'Mullane et al. 2011; Lindegren et al. 2012).

A single star at finite distance and with rectilinear space motion can be described by just five astrometric parameters, representing its position (α, δ), proper motion (μ_α, μ_δ), and parallax (ϖ). Any orbiting companions, including those of planetary mass, will perturb the stellar motion and result in deviations of the individual ("intermediate") astrometric data from a simple five-parameter model. Detectability will depend on the amplitude of the deviations (Section 2.2), and the number and coverage of the individual measurements.

The number of individual field of view crossings, N_fov, along with the final mission accuracies from which they are constructed, are primarily dependent on ecliptic latitude, β. This results from the satellite "scanning law," which is optimized to maintain a constant (solar) thermal payload illumination, while maximizing separability of the astrometric parameters. N_fov is independent of magnitude, and ranges between about N_fov ∼ 60 at β < 10° to about N_fov ∼ 80 at β > 80°, with a maximum of about N_fov ∼ 150 at intermediate ecliptic latitudes, β ∼ 45°, where the scanning density is highest (Table 1). The high values of N_fov around β ∼ ±45° do not necessarily improve planet detection substantively, adding little to the number of distinct epochs and projection geometries.

Table 1. Average Number of Field of View Passages per Star, 〈N_fov〉, versus Ecliptic Latitude

| b |	〈Nfov〉	$\langle N_{\rm fov}^\prime \rangle$
(°)
0–5	64.8	51.9
5–10	65.4	52.3
10–15	67.2	53.7
15–20	69.5	55.6
20–25	73.1	58.5
25–30	78.3	62.6
30–35	86.4	69.1
35–40	99.0	79.2
40–45	134.3	107.4
45–50	138.3	110.6
50–55	106.3	85.0
55–60	95.2	76.2
60–65	88.8	71.0
65–70	84.6	67.7
70–75	81.5	65.2
75–80	79.4	63.5
80–85	78.5	62.8
85–90	77.6	62.1
0–90	86.1	68.9

Note. The column $\langle N_{\rm fov}^\prime \rangle$ includes a mission-level contribution of 20% "dead time."

Download table as:  ASCIITypeset image

Our simulations of detection and orbit reconstruction require estimates of σ_fov, the along-scan accuracy per field of view crossing as a function of G magnitude (Table 2). In terms of the centroiding accuracy for each of the nine astrometric CCDs, σ_η, we adopt

$$\begin{equation} \sigma _{\rm fov} =\left({\sigma _\eta ^2 \over 9} + \sigma _{\rm att}^2 + \sigma _{\rm cal}^2\right)^{0.5}, \end{equation} \tag{ 1 }$$

where σ_att is the contribution from (both random and systematic modeling) attitude errors, and σ_cal is that from calibration errors. Both are assumed constant over the field crossing; we adopt σ_att = 20 μas (Risquez et al. 2013), and similarly for σ_cal (Lindegren et al. 2012). Evidently, both are provisional pending results from the global iterative solution.

For bright stars, G ≲ 12, signal saturation is avoided by CCD "gating," activated according to the star's measured brightness, allowing a reduced number of active TDI lines, and designed to result in a more or less constant measurement precision over the range G ∼ 3–12 mag.

In terms of the inverse relative number of photons in the image

$$\begin{equation} z = 10^{\,0.4(\max [G,12]-15)}, \end{equation} \tag{ 2 }$$

normalized to z = 1 at G = 15, we use

$$\begin{equation} \sigma _\eta = ({53{,}000}\,z + 310\,z^2)^{0.5}, \end{equation} \tag{ 3 }$$

which is a fit to the values 92, 230, 590, 960, 1600, 2900 μas quoted for G = 13, 15, 17, 18, 19, 20 by Lindegren et al. (2012, Table 1).

We complete these forms with an expression for the sky-averaged parallax accuracy, which can be approximated by

$$\begin{equation} \sigma _\varpi = 1.2 \times 2.15\; \sigma _{\rm fov} / \sqrt{68.9} \,= 0.311\, \sigma _{\rm fov}, \end{equation} \tag{ 4 }$$

where 2.15 is the geometric factor linking the (sky-averaged) parallax accuracy with the error per field crossing, 68.9 is the (sky-averaged) number of field crossings per star over the nominal 5 yr mission including dead time (Table 1), and 1.2 is a margin (Lindegren et al. 2012).

Values of z, σ_η, σ_fov, and σ_ϖ, as a function of G magnitude, are given in Table 2.

As a planet detection and characterization technique, astrometry aims to measure the influence of an orbiting planet in addition to the two other classical astrometric effects: the linear path of the system's barycenter projected on the sky (the star's proper motion), and the reflex motion (the star's parallax) resulting from Earth's orbital motion around the Sun. Both the star and the planet orbit the star–planet barycenter and, after accounting for the parallax and proper motion terms, the orbit of the primary therefore appears projected on the plane of the sky as an ellipse with semi-major axis given by

$$\begin{equation} a_\star =\left(\frac{M_{\rm p}}{M_\star }\right) a_{\rm p}, \end{equation} \tag{ 5 }$$

where M_p and M_⋆ are the planet and star mass, respectively, and a_p is the semi-major axis of the planet orbit with respect to the barycenter.

The observable for astrometric planet detection is the corresponding quantity in angular measure, generally referred to as the astrometric signature, given by

$$\begin{equation} \alpha = \left(\frac{M_{\rm p}}{M_\star }\right) \left(\frac{a_{\rm p}}{\rm 1\, AU}\right) \left(\frac{d}{\rm 1 \,pc}\right)^{-1} \hspace{-5.0pt} \,{\rm arcsec}, \end{equation} \tag{ 6 }$$

where d is the distance, and M_p and M_⋆ are in common units. The definition may also be adopted for e ≠ 0, but with detectability dependent on orbital phase (Section 4). The effect is linearly proportional to a_p and, importantly, applies equally to hot or rapidly rotating stars. However, while the technique is most sensitive to massive planets at large a_p, measurement timescales must be proportionally long (of the order of the orbital period).

The size of the effect calculated for all confirmed exoplanets to date (2014 September 1) is shown in Figure 1(a) as a function of orbit period. Vertical lines illustrate the period limits between which Gaia will be most efficient in its discovery space (0.2 ≲ P ≲ 6 yr; see Section 5). On the assumption that an astrometric signature of ∼1–3 times the parallax standard error could be detected (see Section 4), a sizeable fraction of known systems will have their exoplanet-induced photocentric motion determined at some level by Gaia.

Zoom In Zoom Out Reset image size

Figure 1(b) restricts the plot to the known transiting planets, and demonstrates that Gaia astrometry will provide little orbital information for the majority of known transiting planets. No transiting planets have α > 30 μas, and the great majority have α ≪ 1 μas. Indeed, known exoplanets with large α are almost exclusively those discovered by radial velocity measurements. Even the nearest hot Jupiters will be undetected astrometrically, and the same applies to the P ≲ 6 day planets that might be discovered in the Gaia photometric data (Dzigan & Zucker 2012). Gaia astrometry may nonetheless clarify the existence of massive outer companions of hot Jupiters, which have been invoked to explain their inward migration (e.g., Bakos et al. 2009; Neveu et al. 2013; Knutson et al. 2014).

Various astrophysical noise sources will contribute to the accuracy of astrometric measurements in principle, but appear to lie below relevant limits in practice, and have been ignored here. These include the effects of variable stellar surface structure (star spots, plages, granulation, and non-radial oscillations) on the observed photocenter (e.g., Reffert et al. 2005; Ludwig 2006; Eriksson & Lindegren 2007; Lanza et al. 2008; Makarov et al. 2009), relativistic modeling at ∼1 μas (e.g., Anglada-Escudé et al. 2007), and possible effects at optical wavelengths of interstellar and interplanetary scintillation and stochastic gravitational wave noise (Perryman et al. 2001).

A three-dimensional (3D) Keplerian orbit is described by seven parameters, for example, the classical elements a, e, P, t_p, i, Ω, ω. The semi-major axis a and eccentricity e specify the size and shape of the orbit. The period P is related to a and the component masses through Kepler's third law, while t_p specifies the position of the object along its orbit at some reference time, generally with respect to a specified pericenter passage. The three angles (i, the orbit inclination to the plane tangent to the celestial sphere; Ω, the longitude of the ascending node; and ω, the argument of pericenter) give the projection of the true orbit into the observed (apparent) orbit; they depend solely on the orientation of the observer with respect to the orbit.

Both radial velocity and astrometry measure the host star's barycentric motion rather than that of the planet directly, and somewhat different information is provided by each measurement technique.⁸

From astrometry, all seven Keplerian elements are accessible in some form, irrespective of the orbit inclination to the line-of-sight. Specifically, the orbit solution (including the planet location along the orbit as a function of time) gives i and α. From Equation (6), a_⋆ can be determined from α if d is known, with a_⋆ + a_p (and hence a_p) obtained from Kepler's third law assuming that (M_⋆ + M_p) ≃ M_⋆ can be estimated from the star's spectral type or from evolutionary models. Then M_p is determined from Equation (5). If the planet is invisible (the case for all but a few more massive long-period planets that have been imaged) the orbital motion of the star around the system barycenter is correctly determined by astrometry only if the star position is measured with respect to an "absolute" reference frame, which is the case for Gaia (Perryman et al. 2014). Astrometric measurements alone are, however, unable to identify which of the nodes is ascending, i.e., where the planet moves away from the observer through the reference plane, an ambiguity resolved by radial velocity observations.

For multiple exoplanet systems, and if the orbital contributions from each can be separated, astrometry can also establish the relative inclination between pairs of orbits (e.g., van de Kamp 1981, Equation (16.5); Casertano et al. 2008; McArthur et al. 2010).

Four orbit elements (a_⋆, e, ω, t_p) are in common between astrometric and spectroscopic orbit solutions. Combined observations therefore further constrain and improve the 3D orbit, as well as the individual component masses (e.g., Wright & Howard 2009).

As an input to a new set of simulations, we used the population synthesis Galaxy star count model TRILEGAL (TRIdimensional modeL of thE GALaxy; Girardi et al. 2005, 2012). This is based on a theoretical stellar luminosity function $\phi (M,\pmb {r},\lambda)$ (i.e., as a function of absolute magnitude M, Galactic position $\pmb {r}=(\ell,b,r)$, and photometric passband λ), derived from a set of evolutionary tracks, together with suitable distributions of stellar masses, ages, and metallicities. TRILEGAL (version 1.6) includes five distinct Galaxy components: the thin and thick disks, the halo, the bulge, and the disk extinction layer (Girardi et al. 2005, Section 3.6). The model has been calibrated with respect to a variety of observational counts, including multi-passband catalogs from very deep galaxy surveys (including CDFS, DMS, and SGP), the "intermediate-depth" near-infrared point source catalog Two Micron All Sky Survey (2MASS), and the local stellar sample derived from Hipparcos.

A run of TRILEGAL is formally a Monte Carlo simulation in which stars are generated according to specified probability distributions. The number of stars in each bin of distance modulus is predicted according to (Girardi et al. 2005, Equation (1))

$$\begin{equation} N(m_\lambda, \ell, b) = {\rm d}m_\lambda \int_0^\infty {\rm d}r\, r^2 \,\rho (\pmb {r})\, \phi (M_\lambda,\pmb {r})\, {\rm d}\Omega. \end{equation} \tag{ 7 }$$

For each simulated star, the star formation rate, age–metallicity relation, and initial mass function are used to derive the stellar age, metallicity, and mass. Absolute photometry is derived via interpolation in the grids of evolutionary tracks (or isochrones), and converted to the apparent magnitudes using the appropriate values of bolometric corrections, distance modulus, and extinction. All relevant stellar parameters can be retained from the simulations, including the initial and current mass, age, metallicity, surface chemical composition, surface gravity, luminosity, and effective temperature.

The photometric system is an option in the simulation input. We selected Sloan ugriz combined with 2MASS JHK_s. The broadband Gaia G magnitudes are then estimated as a function of Sloan g and z as (Jordi et al. 2010, Table 5)

$$\begin{equation} G \!=\! g - 0.1154 - 0.4175\,(g - z) - 0.0497\,(g - z)^2 + 0.0016\,(g - z)^3, \end{equation} \tag{ 8 }$$

with σ = 0.08 mag over a broad range of colors and extinction. After a number of investigation runs, the magnitude limit of our simulations was set to r = 17.5 mag as a compromise between retrieving a representative selection of low-luminosity stars with large astrometric signature (for reference, a magnitude limit of r = 15 at d = 200 pc corresponds to M = 0.38M_☉), while keeping the computational demands at a reasonable level.

We used the current version of TRILEGAL (version 1.6), with a perl script provided by L. Girardi to automate the interaction with the www interface (stev.oapd.inaf.it/trilegal). Default model values and normalizations were used for the various Galaxy components: (1) an initial mass function (IMF) given by the Chabrier log normal distribution, and a binary fraction of 0.3 with mass ratio in the range 0.7–1; (2) extinction according to an exponential disk of scale height 110 pc, scale length 100 kpc, and A_V(∞) = 0.0378 mag; (3) solar position R_☉ = 8.7 kpc, z_☉ = 24.2 pc; (4) thin and thick disks given by ρ∝exp (− R/h_R) f(z), where the vertical distribution is given by f(z) = sech²(0.5z/h_z), and where the vertical scale height of the thin disk is assumed to increase with stellar age as h(t) = z₀(1 + t/t₀)^α; (5) the halo is defined by an oblate r^1/4 spheroid, and the Galaxy is assumed to comprise a triaxial bulge; (6) default selections for the star-formation rate and age-metallicity relation were also used for each Galaxy component (thus for the thin-disk: two-step star-formation rate, with the age–metallicity relation from Fuhrmann (1998), and α-enhancement; for the thick-disk: 11–12 Gyr constant star formation, z = 0.008 with 0.1 dex standard deviation, and solar-scaled abundances; for the halo: 12–13 Gyr age, with [M/H] distribution from Ryan & Norris (1991); for the bulge: 10 Gyr age, with [M/H] distribution from Zoccali et al. (2008) enhanced by 0.3 dex).

Simulations were run for a field size of 1 deg² in 10° steps of Galactic longitude ℓ = 0°–180° (being symmetric for 180°–360°) and Galactic latitude b = −90° to +90° (excluding b = 0° where the simulations did not complete within the time imposed by the server, but including b = ±5° and b = ±15°). For the main part of our simulations, we restricted our selection to stars with log g > 3.0 and log T_eff < 4.0, thus excluding giant and high-mass stars (this approximation is discussed further in Section 7.2). This yielded a total of ∼915,000 stars in the mass range 0.07–3.27M_☉. The full-sky sample was then generated by interpolating these numbers over a mesh of 0.01 rad in ℓ and b, finding the nearest simulated field to each grid node, and drawing the specified number of samples, with replacement, from that file. Each sample is then a simulated star.

The result is a list of N_⋆ ∼ 260 × 10⁶ stars, including binaries, whose distance distribution is shown in Table 4 ($N_\star ^{\rm P14}$). Fluctuations between 300–500 pc are due to the TRILEGAL model using discrete steps in distance modulus. Because the simulations are magnitude-limited, simulated stars actually extend out to 16 kpc. Although we consider all sample stars when simulating planets, only those within 700 pc (for α > 2σ_fov) yield detections. Of the 260 × 10⁶ stars simulated, ∼20 × 10⁶ are within 700 pc (Table 4).

In simulations using the Besançon Galaxy model (version 2003, as available on the web site) to replicate the detection numbers reported by Casertano et al. (2008), we found that the TRILEGAL star counts were some 30% lower out to ∼200 pc than those returned by the Besançon galaxy model. A. Robin (2014, private communication) has attributed this to the fact that the star formation rate in the Besançon model was assumed constant over the thin disk life time; the revised model described by Czekaj et al. (2014) is expected to correct much of this discrepancy. We have accordingly adopted the results from TRILEGAL, noting that some 30% more exoplanets would be detected astrometrically should the star counts (and the treatment of extinction) follow more closely that predicted by the Besançon model.

Recent developments in characterizing exoplanet occurrence frequencies as a function of planetary and stellar properties are directly relevant to our re-assessment. These include (1) improved statistics of the longer-period radial velocity discoveries, allowing a more secure extrapolation for large M_p and P; (2) larger numbers of Neptune-mass discoveries, which augments the contribution of lower-mass planets at large a_p (while the distant limit for detection will decrease compared with those of Jupiter mass, this will be compensated by the increased numbers at smaller d); (3) improvements in characterizing the occurrence of planets around M dwarfs. At the same time, we have introduced some very simple, and intentionally conservative, assumptions in cases where occurrence rates are unknown or at best poorly known empirically.

More specifically, for each of the 260 million stars returned by the TRILEGAL simulations to r < 17.5 mag, we simulate the exoplanet occurrence and properties according to the following dependencies.

• 1.  
  Binary stars: the secondary stars of binaries were ignored, and planets simulated around the primaries according to the occurrence distributions for single stars (this simplification is discussed in Section 7.2).
• 2.  
  Host star mass and metallicity: for giant planets around host stars with M_⋆ > 0.6M_☉, we used the relationship between giant planet occurrence and host star mass and metallicity given by Johnson et al. (2010). This relation was determined for stars with 0.5M_☉ < M_⋆ < 2.0 M_☉ and [Fe/H] < 0.4. We assume that the occurrence rate of giant planets around stars with M_⋆ > 2.0 M_☉ is equal to that at M_⋆ = 2.0 M_☉, and that the rate for stars with [Fe/H] > 0.4 is equal to that at [Fe/H] = 0.4. We extrapolate occurrence rates to stars off the main sequence (including white dwarfs), while excluding from consideration giant and high-mass stars (as noted in Section 3.1 and discussed further in Section 7.2). Our treatment of host star mass and metallicity for smaller planets and for planets around stars with M < 0.6 M_☉ is discussed below.
• 3.  
  Planet mass and orbital period: we assume that the planet mass and period distributions are independent of metallicity. For each star we draw planets from the joint mass–period distribution given by Cumming et al. (2008), as determined from radial velocity surveys for GK stars. For planet masses >0.3M_J and periods <2000 d, they obtained a power-law fit dN = C M^α P^β dln M dln P with α = −0.31 ± 0.2, β = 0.26 ± 0.1, and the normalization constant C such that 10.5% of solar-type stars have a planet with mass in the range 0.3–10M_J and orbital period 2–2000 days. For M_⋆ > 0.6 M_☉, we extrapolated this power law to cover the extended mass range 0.1–15M_J, and to orbital periods up to 10 yr (for very wide orbits around A stars, the best constraints are currently based on direct imaging, but the semi-major axes probed lie well beyond the range of orbit period relevant for Gaia). For P < 418 days, the outer period to which the Kepler distributions were determined, we did not extrapolate the Doppler-based distributions below 0.3M_J, and instead used the Kepler results for planets of lower mass (Table 3).
• 4.  
  Occurrence around low-mass stars (M dwarfs): we examined two dependencies. The first is a simplified extension of the Johnson–Cumming distributions, for which we assumed that the occurrence rate of gas giant planets around stars with M_⋆ < 0.5 M_☉ was equal to that at M_⋆ = 0.5 M_☉. The second (which we have adopted for all subsequent steps) follows the conclusions of Montet et al. (2014), who found that the Cumming et al. mass–period power law does not agree with the microlensing results for M dwarfs (see also Clanton & Gaudi 2014 for a recent analysis reconciling the radial velocity and microlensing planet yields for M dwarfs). From Doppler measurements of 111 M dwarfs, they found a lower occurrence rate compared to higher-mass stars, finding only 6.5% ± 3.0% host one or more massive companions with 0 < a_p < 20 AU and 1M_J < M_p < 13M_J. Accordingly, we adopt as our baseline, and more conservative estimate for M_⋆ < 0.6 M_☉, their double power law in stellar mass and metallicity

  $$\begin{equation} f(M_\star, {\rm [Fe/H]}) = 0.039\, M_\star ^{0.8} \, 10^{3.8\,{\rm [Fe/H]}}, \end{equation} \tag{ 9 }$$

  with a dependency on M_p (and flat in log a_p), consistent with both their observations and the microlensing observations, given by

  $$\begin{equation} dN \propto \ M_{\rm p}^{-0.94} \ d\ln M_{\rm p}\, d\ln a_{\rm p}. \end{equation} \tag{ 10 }$$

  For the same stars, M_⋆ < 0.6 M_☉, we extrapolated this power law to cover the extended mass range 0.1–15M_J. Using the Montet et al. (2014) distribution for M_⋆ < 0.6 M_☉ in place of the Johnson–Cumming distributions substantially reduces the expected Gaia planet yield for two reasons: (1) rather than assuming that stars with M_⋆ < 0.5 M_☉ have the same planet occurrence fixed to the Johnson et al. rate at M_⋆ = 0.5 M_☉, the Montet et al. relation predicts that the planet occurrence rate continues to drop for lower mass stars; (2) although the occurrence rate of small planets around M dwarfs is higher in the Montet et al. distribution, the larger planets and longer period planets which Gaia would be sensitive to are less common.
• 5.  
  Eccentricities: we assumed that eccentricities of the 3D orbits follow a Beta distribution

  $$\begin{equation} P\!_{\beta }\,(e;a,b) \ \propto \ e^{a-1}\, (1-e)^{b-1}, \end{equation} \tag{ 11 }$$

  with a = 0.867 and b = 3.03, as established from fits to radial velocity observations by Kipping (2013). We ignore any possible dependency on period (in practice, the planets detectable with Gaia astrometry are on wide orbits where this assumption, in the absence of tidally circularized orbits, appears justified). The eccentricity was not used in deciding whether an astrometric signal would be detected (Section 4), but it is used in determining the enhanced transit probability for elliptical orbits (Section 6).
• 6.  
  Distribution of low-mass planets: for M_p < 0.3M_J, we used the joint period–radius distribution determined from the Kepler data by Fressin et al. (2013). This extends only out to P = 428 days for the largest planets, and since a power law does not appear to provide a particularly good fit, we did not extend the distributions to longer periods. Radii were converted to masses using the following empirical relation, chosen to match the Fressin et al. (2013) occurrence rate (as a function of R_p) to the Howard et al. (2010) occurrence rate (as a function of M_p) for P < 50 days, and with M_p = 0.3M_J at R_p = 1R_J (the same approach as used by Howard et al. 2012 to compare their radius distribution from Kepler to their earlier mass distribution from Doppler observations; Howard et al. 2010)

  $$\begin{eqnarray} M_{\rm p}&=&1.08\,R_{\,\rm p}^{\,3.45}\qquad\qquad R_{\rm p}\le 1.5\nonumber\\ &=&3.17\,R_{\,\rm p}^{\,0.87}\qquad\qquad 1.5 < R_{\rm p} \le 4.0 \nonumber\\ &=&10.59\,(R_{\,\rm p}/4)^{\,2.07} \qquad 4.0 < R_{\rm p}\le 6.0\nonumber\\ &=&24.51\,(R_{\,\rm p}/6)^{\,2.17}\qquad 6.0 < R_{\rm p} \end{eqnarray} \tag{ 12 }$$

  where M_p and R_p are expressed in Earth units.
• 7.  
  Multi-planet systems: the Cumming et al. (2008) distribution uses only the most significant Doppler signal for stars with multiple planets. Carrying this over into our planet occurrence simulations, we do not account for multiple gas giant systems. The Fressin et al. (2013) distributions based on Kepler, on the other hand, include all planets in the multi-planet systems. In this case, we have partitioned the distribution into their designated mass bins (Table 3), and generate at most one planet per mass bin for each star. In practice, no planets from the parameter space covered by the Fressin et al. (2013) distribution are recovered with α > 2σ_fov, and therefore we do not recover any multi-planet systems.

Table 3 summarizes the fraction of stars with planets from the combination of Kepler transit and radial velocity data, as a function of mass and period. The end result is that for our 260 million simulated stars to r = 17.5, we simulated 254 million planets (with 79 million representing additional planets within a multiple planet system), accompanying 175 million stars.

We stress that our treatment of multi-planet systems is incomplete. First, we do not allow for more than one (low-mass) planet in one mass bin. However, in practice, the entire region of parameter space covered by the Fressin et al. (2013) distribution yields just 1 marginally detectable Gaia planet (with M_p = 0.29M_J, P = 389 d, and α = 1.15σ_fov). To this extent, our simplistic assumptions do not affect the question of astrometric detectability. Second, we have assumed that the astrometric motion of the host star is dominated by a single massive planet. To go further is beyond the scope of this paper: current knowledge of orbit statistics for multiple gas giants at large a_p is limited, while the astrometric motion of the star with respect to the system barycenter for multiple massive planets rapidly becomes more complex (see, e.g., Perryman & Schulze-Hartung 2011, Figures 2–3). We refer to Casertano et al. (2008) for further insight into the detectability (although not the occurrence) of multiple massive systems.

We generate simulated observations of large numbers of exoplanet systems using tools provided by the Gaia AGISLab project (Holl et al. 2012, their Appendix B). Developed as part of the astrometric global iterative solution (AGIS; Lindegren et al. 2012), AGISLab allows the simulation of millions of sources at the level of individual CCD transits, using a comprehensive instrument model, including the scanning law. Using this, we can generate, for any target star based on its sky coordinates, a listing of all field crossings over the mission, giving the time, position angle of the scan, and the parallax factor for each observation. This, together with the assumed along-scan standard error per field crossing (as a function of G), and the seven specified orbital elements of the star's reflex motion, allows us to simulate a full set of representative ("intermediate astrometry") observations.

We then subject these simulated observations for each system to a least-squares orbit fitting algorithm. The objective is to recover the 12 parameters (5 astrometric and 7 Keplerian) describing the star position at each epoch of observation (Figure 4), making the (reasonable) assumption that other deterministic effects incorporated into AGIS (aberration, relativistic light bending, and perspective acceleration) are fully accounted for.

Zoom In Zoom Out Reset image size

For this investigation, we recast the classical orbit elements (a, e, P, t_p, i, Ω, ω) into an equivalent set consisting of the four Thiele–Innes constants (A, B, F, G), together with the frequency f = 1/P, eccentricity e, and mean anomaly at the reference epoch, M₀. While certain aspects of this orbit fitting are rather standard, a number of specific considerations are relevant for Gaia, such as the range of the period search, and the use of a prior density of the orbit eccentricity. Further details are given in Appendix A.

Formal uncertainties of the fitted orbit parameters can be quite misleading, due to the strongly non-linear nature of the fitting procedure, and we use instead a Monte Carlo approach. For a given system (with fixed orbit and observation geometry), we generate N observation sets with independent noise realizations, leading to N different sets of estimated orbit parameters, from which their precision can be estimated. Depending on the kind of investigation, N could range from 1 (when generating the statistics for a large sample of different systems) to 100 or more (when assessing the precision of a specific system).

For illustration, we show results from the Monte Carlo simulations for two transiting systems, i.e., for two systems with "true" i = 90°, and randomly assigned elements ω, Ω, and t_transit. Remaining parameters were simulated as described in Section 3.2.

Our first example has G = 7.8 mag, P = 0.65 yr, e = 0.32, and α = 83 μas, corresponding to a detection at α = 2.4 σ_fov. The number of field crossings (including 20% dead time) is N_fov = 59. Figure 5 shows scatter plots of a_p, P, e, and cos i, and the predicted transit times, for 100 different noise realizations (in all diagrams the long-dashed lines show the true values). The predicted transit times are for ω + ν = 90° or 270°, where ν is the true anomaly. Because astrometry alone cannot determine ω unambiguously, we assume 0 ⩽ ω < 180°, and consequently have to predict two possible transit times per estimated orbit period, but only one true transit time (long-dashed line) per true period.

Zoom In Zoom Out Reset image size

Our second example (Figure 6) has G = 15.5 mag, P = 4.12 yr, e = 0.015, α = 1196 μas (11.5σ_fov), and N_fov = 40. Despite its unusually small number of field crossings (even for its ecliptic latitude of β = −5°; see Table 1) and fainter magnitude (G = 15.5), this orbit is even better determined due to the much higher S/N, and the larger ratio of M_p (=13.8M_J) to M_⋆ (=0.2 M_☉). Interestingly, cos i is about equally well determined in all solutions, even when a_p, P, or e is significantly wrong. We also note that the predictions for t_transit are only good for a few times P around the Gaia observing epoch because of the relatively large uncertainty in P. This will always be a problem for periods larger than a few years, but one that can be improved by radial velocity observations to constrain the period.

Zoom In Zoom Out Reset image size

As discussed in Section 4, the astrometric S/N ratio (α/σ_fov) is a useful zero-order indicator of detectability, but one that in reality depends on the number and distribution of observations, the inclination and eccentricity of the system, and the orbital period in relation to the total length of the observations. We will show that a more precise criterion is given by the likelihood ratio, or, equivalently, the reduction in the minimum χ², when going from the 5 parameter solution to the 12 parameter solution.

Let $\chi ^2_{\rm min}({\rm 12\, parameter})$ be the minimum χ² obtained when adjusting all 12 parameters (Equation (A5)). Omitting the orbit parameters (i.e., setting $\boldsymbol{y}=\boldsymbol{0}$ in Equation (A5)) and fitting only the astrometric parameters ($\boldsymbol{x}$) results in a fit with $\chi ^2\,({\rm 5\, parameter}) \ge \chi ^2_{\rm min}({\rm 12\, parameter})$. The increase in χ² when omitting the orbit parameters is

$$\begin{equation} \Delta \chi ^2 = \chi ^2_{\rm min}({\rm 5\, parameter}) - \chi ^2_{\rm min}({\rm 12\, parameter}), \end{equation} \tag{ 14 }$$

which can therefore be used as a test statistic for the significance of the orbit. In contrast to the S/N, this quantity can be calculated without knowledge of the orbit and is therefore applicable to real data. We note that Δχ² can be small even when α/σ_fov is very large, e.g., if the observations cover only a small part of the orbit. On the other hand, for a fixed S/N per field crossing, Δχ² increases with the number of observations.

Since exp (− χ²/2) is proportional to the likelihood of the model (assuming Gaussian observation noise with the stated standard deviation), this is effectively a likelihood ratio test and therefore close to optimal in terms of its power to detect orbital motion. Based on Wilks' theorem (e.g., Kendall et al. 1983), the distribution of Δχ² in the absence of a companion is expected to follow the chi-squared distribution with seven degrees of freedom. This provides a useful guide for setting the detection threshold. A possible criterion could be Δχ² > 30, with a reasonably low probability of false detection (theoretically ∼10⁻⁴). Higher values of Δχ² give both a more reliable detection and a higher precision of the estimated orbit. In practice, solutions with Δχ² ≃ 30 may be considered marginal, while those with Δχ² > 50 are generally found to be reliable and Δχ² > 100 typically gives orbital parameters determined to 10% or better. The two examples shown in Figures 5 and 6 have Δχ² ≃ 42 and 970, respectively.

For our present primary purpose of assessing detectable numbers, we proceed by provisionally selecting systems above a certain S/N per field crossing (e.g., α > 2σ_fov; see Table 4). We then confirm these provisional detections by also requiring that the fit improves significantly, as measured by Δχ², when proceeding from a 5 parameter to a 12 parameter solution.

The estimation of the number of detectable systems is vastly sped up by using λ + 7 as a proxy for Δχ². Here, λ is the noncentrality parameter (Appendix B) which can be calculated, for simulated data, without actually fitting an orbit. All statistics reported below for Δχ² > X were in fact derived using the criterion λ + 7 > X.

This additional Δχ² criterion is reflected in the bottom part of Table 4. We can see, for example, that the 16,668 planets provisionally detected according to α > 2σ_fov result in 12,893 secure detections according to Δχ² > 30, or to 10,297 detections according to the more restrictive Δχ² > 50. This substantiates our claim (Section 4) that α > 2σ_fov rather than α > 3σ_fov provides a reasonable zero-order estimate of the numbers detectable. In practice, final detection results are then somewhat insensitive to the actual choice of α/σ_fov, in the sense that additional provisional candidates revealed by lowering the S/N threshold are in any case subject to confirmation by the Δχ² criterion. Below α ≲ 0.5 σ_fov more candidates naturally continue to be selected, but the vast majority fail to result in additional confirmed detections.

Our best estimates, pre-selected with α > 0.5 σ_fov (lower thresholds are evidently required when assessing results for a 10 yr mission duration), are indicated in bold in the lower part of Table 4. At Δχ² > 50–30 we find 14,806–27,505 (∼21, 000 ± 6000) systems discoverable over the nominal 5 yr mission duration. As we will quantify in Section 6, some 25–42 of these are likely to be transiting.

Gaia has the potential to observe for considerably longer than the nominal 5 yr (currently limited by its cold gas attitude control mass), and a longer mission would bring very substantial improvements for exoplanet detection. Specifically, (1) it will be possible to detect and reconstruct orbits with periods P ≳ 5 yr; (2) a larger number of systems with P < 5 yr will be detected; (3) the number of false detections will decrease; and (4) the accuracy of the orbit solutions will be greatly improved. To quantify this we have repeated the orbit determinations for the 3500 simulated systems described in Section 6 (viz., 100 realizations of 35 predicted transiting astrometric detections with S/N >2) and an assumed mission length of 10 yr. Resulting Δχ² are typically a factor of 3–4 larger.

Table 4 contains our estimated detection numbers at various levels of S/N and Δχ² (to obtain complete statistics S/N down to 0.5 per field crossing must be considered). We find that about three to four times as many systems could be detected with a comparable Δχ², and that the subset of transiting systems would increase by the same factor. We find 90,751, 53,015, and 25,958 detections at Δχ² > 30, 50, and 100, respectively, with the corresponding numbers of transiting planets being 135, 82, and 31, respectively.

Use of a population synthesis Galaxy model means that we have estimates of the numbers of host stars also as a function of spectral type. In this section we expand on the results for FGK dwarfs and M dwarfs, comparing our findings to previous evaluations.

Rather than attempting detailed predictions in view of largely unknown occurrence rates around rare spectral types more generally, we simply emphasize that the Gaia census should place many new constraints on the wider occurrence of planetary systems across the entire HR diagram. For example, some 1.2% of our recovered planets accompany very low-mass pre-main sequence stars, while some 0.1% (13 out of 16,668 down to r < 17.5 in our simulations at 2σ_fov) are around white dwarfs (identified in the M_⋆ − R_⋆ plane as R_⋆ ∼ R_⊕). The numbers around white dwarfs may be conservative since we have used the current M_⋆ (viz. 0.64–0.74 M_☉) to scale the planet occurrence rate, rather than the mass of the main-sequence progenitor. More significantly, the expected numbers of white dwarfs detectable by Gaia is a strongly increasing function of the limiting magnitude: estimates suggest that some 23,000 disk white dwarfs and some 1100 halo white dwarfs will be detected to G = 20 (Figueras et al. 1999; Carrasco et al. 2014), with the majority of these within 300–400 pc. We have not extended our simulations to r > 17.5 mag since the planet occurrence rate around white dwarfs, and especially over the relevant values of M_p and P, are essentially unknown. The Gaia results will provide an effective exoplanet survey for these nearby systems and will establish, for example, whether the occurrence of gas giants in wide orbits around white dwarfs is consistent with that around main-sequence stars.

7.1.1. FGK Dwarfs

7.1.2. M Dwarfs

M dwarfs provide an interesting subset of host stars for more detailed consideration, given their large numbers at relatively small distances from the Sun, their low stellar mass resulting in large astrometric displacements for a given planet mass, and their topical interest as habitable zone host stars (e.g., Bonfils et al. 2013). Of the 16,668 exoplanets recovered at 2σ_fov for our TRILEGAL magnitude limit of r < 17.5 mag (Table 4), 2097 are around M dwarfs, identified from our simulation as stars with M_⋆ < 0.6 M_☉.

Sozzetti et al. (2014) made a detailed evaluation of the number of planets detectable around M dwarfs with Gaia, focusing on the 33 pc distance limit of the volume-limited LSPM sample (Lépine 2005), and on an extrapolated volume-limited sample out to 100 pc. They adopted magnitude-dependent single measurement errors from recent Gaia models with a typical error of σ_fov = 100 μas for G ≲ 16 mag. They used similar criteria for astrometric detection and orbit characterization as derived by Casertano et al. (2008). They assumed a distribution of planet occurrences defined by M_p = 1M_J, P uniformly distributed over 0.01–15 yr, and e uniformly distributed over 0–0.6, with an overall frequency of f_p = 0.03 (as given by Johnson et al. 2010).

As summarized in Table 5, they predicted some 100 planets detected at 3σ within the 33 pc horizon of the volume-limited LSPM sample (noting that LSPM is neither complete nor entirely comprised of M dwarfs), and 2600 detected planets (and 500 accurate orbits) out of an estimated 4 × 10⁵ objects within 100 pc based on the Besançon Galaxy model. They identified a statistical subset of transiting systems based on a random distribution of orbit inclinations, estimating that ∼40 transiting objects would be sampled out to 100 pc, with accurate orbits for ∼10 (with uncertainties on i of ≈2–10° for i ∼ 90°).

Table 5. Planets Around M Dwarfs Detectable by Gaia, for Two Distance Limits from the Sun

Δd		$N_\star ^{\rm S14}$	$N_{\rm det}^{\rm S14}$		N⋆		Ndet	Ntran		Ndet	Ntran		Ndet	Ntran		Ndet	Ntran		Ndet	Ntran
(pc)								α > 0.5 σfov		α > 1 σfov			α > 2 σfov			α > 3 σfov			α > 6 σfov
0–33(L)		3 150	100		3 150		257	1.2		167	0.6		93	0.3		64	0.2		32	0.1
0–33(C)		⋅⋅⋅	⋅⋅⋅		9 944		344	1.2		225	0.6		122	0.3		91	0.2		42	0.1
0–100(C)		415 000	2600		191 567		4149	11.1		2090	4.7		986	1.8		635	1.1		267	0.4
					Δχ2		Ndet	Ntran		Ndet	Ntran		Ndet	Ntran		Ndet	Ntran		Ndet	Ntran
							α > 0.5 σfov			α > 1 σfov			α > 2 σfov			α > 3 σfov			α > 6 σfov
0–33 pc (C), 5 yr mission					>30		157	0.3		154	0.3		114	0.3		90	0.2		42	0.1
''					>50		121	0.2		121	0.2		107	0.2		87	0.2		42	0.1
''					>100		87	0.2		87	0.2		85	0.2		78	0.2		41	0.1
0–100 pc (C), 5 yr mission					>30		1451	2.4		1415	2.4		912	1.7		625	1.1		267	0.4
''					>50		1047	1.7		1047	1.7		840	1.5		595	1.1		263	0.4
''					>100		644	1.1		644	1.1		625	1.1		517	1.0		255	0.4
0–33 pc (C), 10 yr mission					>30		259	0.5		219	0.5		122	0.3		91	0.2		42	0.1
''					>50		193	0.4		187	0.4		122	0.3		91	0.2		42	0.1
''					>100		140	0.3		140	0.3		118	0.3		91	0.2		42	0.1
0–100 pc (C), 10 yr mission					>30		2609	4.3		2033	3.9		986	1.8		635	1.1		267	0.4
''					>50		1862	3.1		1782	3.1		986	1.8		635	1.1		267	0.4
''					>100		1246	1.9		1246	1.9		956	1.7		635	1.1		267	0.4

Notes. $N_\star ^{\rm S14}$ and $N_{\rm det}^{\rm S14}$ are from Sozzetti et al. (2014), for which their predicted number of astrometrically detected transiting systems are 0 (to 33 pc) and 40 (to 100 pc), respectively. For the 0–33 pc volume, results are given for both the LSPM sample (L; Lépine 2005, with N_⋆ = 3150) and for the predicted complete volume (C), for which N_⋆ is estimated from Galaxy models (TRILEGAL here, and Besançon for S14). Other details are as Table 4. Again, our best estimates of the range of N_det and N_tran are given in bold.

Download table as:  ASCIITypeset image

To compare these results with those from our own simulations, we proceeded as follows. Following Sozzetti et al. (2014), we selected the 3150 reddest (in V−J) non-subgiant sources from the LSPM sample (not all are M dwarfs, but Lépine 2005 does mention that some K dwarfs are in this sample as well). Our own estimates of planet occurrences rates around M dwarfs are dependent on both stellar mass and metallicity (Equation (9)). For the former, we followed Sozzetti et al. (2014) in using the mass–luminosity (in M_J) relation. We estimated G magnitudes based on the LSPM V and J values, adopting V − I = 0.8(V − J) − 0.331 as a fit to the M dwarfs in Lépine & Gaidos (2011, Table 2), then using the G(V − I) relation by Jordi et al. (2010, Table 6). Since metallicities are not available for all LSPM stars, we drew metallicities from the distribution of [Fe/H] values in the solar neighborhood given by Haywood (2001, Table 3). For R_⋆, required to determine transit probabilities, we used an analytic fit to the M_⋆–R_⋆ relation at 4.5 Gyr, and solar metallicity, from the Dartmouth stellar models (Dotter et al. 2008).

Finally, as detailed in Section 3.2, we used the Montet et al. (2014) planet occurrence distribution for M < 0.6 M_☉, and the Johnson–Cumming distribution for M > 0.6 M_☉ (79 of our simulated stars have M_⋆ > 0.6 M_☉, with the highest mass star being 0.85 M_☉). We then simulated planets as before, estimating the number of detections, and the subset of transiting systems. For the corresponding numbers for complete samples of M dwarfs within 33 pc and 100 pc, we simply select stars from our TRILEGAL simulations with M_⋆ < 0.6 M_☉ within these distance limits, and run our detection and transit simulations as before.

The various results, including those based on our Δχ² detection metric, are collected in Table 5, in the same form as for Table 4, for two distance limits (d < 33 pc and d < 100 pc), and for two mission durations (the nominal 5 yr, and a hypothetical 10 yr).

For the LSPM sample our estimate of the planets detected at 3σ is comparable to, if a little lower than, the estimates by Sozzetti et al. (2014), viz. 64 and 100, respectively.

For our complete sample within 33 pc, our number of M stars increases by a factor of 3, while the detected planet numbers increase by only 25%: this is because most of the additional stars have very low masses (63% have M < 0.2 M_☉, compared with 27% in LSPM), while the Montet et al. (2014) relation identifies a decreasing planet occurrence for decreasing stellar mass. Selecting detections through the combination of a lower S/N threshold (α/σ_fov > 0.5) and a Δχ² exceeding 50 and 30, respectively, we estimate 121–157 detections in the complete sample out to 33 pc, with formally only 0.2–0.3 transiting.

For the complete sample within 100 pc, our TRILEGAL-based estimates of the total numbers of M dwarfs (∼190,000) are lower than those of both Sozzetti et al. (2014; ∼415,000), and those implied by the northern 100 pc M dwarf census of Lépine & Gaidos (2013), from which we infer a full-sky estimate of ∼270,000 (including a contribution of 30% due to incompleteness). Our corresponding numbers of planet detections at >3σ_fov are also somewhat less than those estimated by Sozzetti et al. (635 compared to 2600): while Sozzetti et al. (2014) predicted ∼40 transiting systems, we find essentially none (1.1 ± 1.1). Including our Δχ² metric exceeding 50 and 30, respectively, we find 1047–1451 detections out to 100 pc, with just 1.7–2.4 transiting. Again, it can be seen from the lower part of the table that performances improve significantly for an extended 10 yr mission.

These differences in the number of M dwarfs, and therefore in the number of likely planets out to 100 pc, arises in part from our magnitude limit of r = 17.5 (compared to G = 20 adopted by Sozzetti et al. 2014), such that our host star sample is incomplete for the lowest mass M dwarfs (<0.17 M_☉). Other differences originate from the Galaxy models, and the uncertainty on the known planet occurrences. Since Gaia will quantify all of these, we take this analysis no further, concluding that the Sozzetti et al. (2014) estimates are likely to be more complete, such that our estimates of the detectable planets around M dwarfs out to 100 pc (1047–1451) is probably underestimated, perhaps by a factor of 2.

The uncertainties in our numbers remain rather large (with one of the goals of Gaia being to resolve these), perhaps by a factor of 2 or more, due to a number of reasons.

7.2.1. Galaxy Model

7.2.2. Planet Occurrence

7.2.3. Instrument Performance

7.2.4. Bright Star Limit

7.2.5. Mission Accuracy