Toward an initial mass function for giant planets

TL;DR

This paper reconstructs the initial mass distribution of exoplanets by analyzing current eccentricities and using Bayesian methods to account for planets lost through dynamical interactions, revealing that ejected planets are mostly sub-Saturns.

Contribution

It introduces a Bayesian framework combining simulations and observations to infer the primordial exoplanet mass function, accounting for dynamical ejections and collisions.

Findings

Ejected planets are mainly sub-Saturns.

Present-day bimodal distribution is not primordial.

Intermediate-mass planets are preferentially removed by dynamical instability.

Abstract

The distribution of exoplanet masses is not primordial. After the initial stage of planet formation is complete, the gravitational interactions between planets can lead to the physical collision of two planets, or the ejection of one or more planets from the system. When this occurs, the remaining planets are typically left in more eccentric orbits. Here we use present-day eccentricities of the observed exoplanet population to reconstruct the initial mass function of exoplanets before the onset of dynamical instability. We developed a Bayesian framework that combines data from N-body simulations with present-day observations to compute a probability distribution for the planets that were ejected or collided in the past. Integrating across the exoplanet population, we obtained an estimate of the initial mass function of exoplanets. We find that the ejected planets are primarily sub-Saturn type planets. While the present-day distribution appears to be bimodal, with peaks around $\sim 1 M_{\rm J}$ and $\sim 20 M_\oplus$, this bimodality does not seem to be primordial. Instead, planets around $\sim 60 M_\oplus$ appear to be preferentially removed by dynamical instabilities. Attempts to reproduce exoplanet populations using population synthesis codes should be mindful of the fact that the present population has been depleted of intermediate-mass planets. Future work should explore how the system architecture and multiplicity might alter our results.

Figures

Tables

Table 1. Table 1: An overview of the exoplanet dataset used in this work. In this work we assume that planets with μ < 0.05 𝜇 0.05 \mu<0.05 or e < 0.05 𝑒 0.05 e<0.05 acquired their eccentricity without scattering another planet. Our data selection is described in section 4.1 .

		$\mu >0.05$	$\mu \le 0.05$
Total number of exoplanets		506	47
Eccentric ($e\ge 0.05$)		427	15
Discovered by radial velocity		461	23
Discovered by transit		38	23
Discovered by other methods		7	1
Members of multiple-giant systems		171	37
Around O stars	($M_{\star }/M_{\odot }>16$)	0	0
Around AB stars	()	92	0
Around FGK stars	($0.5\le M_{\star }/M_{\odot }\le 1.5$)	388	39
Around M stars	()	26	8

Table 2. Table 2: List of planets included in our sample that were not discovered by either radial velocity or the transit method. The first planet has μ < 0.05 𝜇 0.05 \mu<0.05 .

Planet	mass ($M_J$)	period (day)	semimajor axis (AU)	eccentricity	star mass ($M_{\odot }$)	Discovery method
KOI-620.02	0.024	130.19	0.51	0.008	1.0	TTV
51 Eri b	7.0	14965.0	14.0	0.21	1.75	Imaging
HD 176051 b	1.5	1016.0	1.76	0.0	0.9	Astrometry
KIC 7917485 b	11.8	840.0	2.03	0.15	1.63	Pulsation phase modulation
Kepler-419 c	7.3	675.47	1.68	0.184	1.22	TTV
Kepler-46 c	0.376	57.004	0.2799	0.0145	0.902	TTV
OGLE-2006-109L c	0.271	4931.0	4.5	0.15	0.51	Microlensing
$\beta$ Pic b	7.0	13288.0	13.18	0.323	1.73	Imaging

Table 3. Table 3: Comparison of Θ Θ \Theta and θ 𝜃 \theta (equations 15 and 16 ) for the giant planets in the solar system. The advantage of θ 𝜃 \theta is that it can be calculated for exoplanets without measured radii.

Planet	$\Theta$	$\theta$
Jupiter	5.37	5.20
Saturn	3.57	4.29
Uranus	2.48	2.45
Neptune	4.72	4.28

Table 4. Table 4: We divide the exoplanet population into three sets (A, B, C) based on the normalized mass μ 𝜇 \mu Eqn.  12 (see figure 3 ). For each set we conduct N-body simulations with a planet mass and semimajor axis ( m 1 subscript 𝑚 1 m_{1} , a 1 subscript 𝑎 1 a_{1} ) representative of that set – the stellar mass is always 1 ​ M ⊙ 1 subscript 𝑀 direct-product 1M_{\odot} . We conduct at least 1,000 simulations per set and select the runs that had an instability after 0.5 Myr, and made sure that there are more simulations than observed exoplanets.

	Exoplanet observations		Simulations
Set	Mass range	Exoplanets	$m_1/M_J$	$a_1/\mathrm{au}$	Usable runs
A		55	0.07	0.4	274
B		241	1	1	410
C	$\mu \ge 2.65$	131	7	2	286

Equations

$$E$$

$$\Lambda$$

$${\alpha }^{-3}\left({\mu }_1+\frac{{\mu }_2}{{\delta }^2}\right)({\mu }_1{\gamma }_1+{\mu }_2{\gamma }_2\delta {)}^2>1+3^{4/3}\frac{{\mu }_1{\mu }_2}{{\alpha }^{4/3}}$$

$$\delta$$

$${\gamma }_k$$

$${\mu }_k$$

$$\alpha$$

$$\Delta$$

$$R_H$$

$${\log }_{10}\left(\frac{t_{\mathrm{ce}}}{\mathrm{yr}}\right)\sim -5-{\log }_{10}\left(\frac{\mu }{10^{-3}}\right)+1.44\Delta {\left(\frac{\mu }{10^{-3}}\right)}^{\frac{1}{12}},$$

$$P(\theta |D)\propto P(D|\theta )P(\theta ),$$

$$\mu =\frac{m/M_{\star }}{M_{\mathrm{J}}/M_{\odot }},$$

$$q\in \{0.1,0.2,0.3,\cdots ,1.0\}\mathrm{and}{m}_2=q{m}_1$$

$$\sigma =\pi R^2\left(1+\frac{v_{\mathrm{esc}}^2}{v_{\infty }^2}\right)=\pi R^2\left(1+\Theta \right),$$

$$\Theta \sim \frac{1}{2}\left(\frac{m}{M_{\star }}\right)\left(\frac{a}{R}\right).$$

$$\theta ={\left(\frac{m}{M_{\mathrm{J}}}\right)}^{2/3}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{-1}\left(\frac{a}{\mathrm{au}}\right)$$

Full text

††pubyear: 2016††pagerange: Toward an initial mass function for giant planets–Toward an initial mass function for giant planets

1 Introduction

The past twenty years of exoplanet observations have revealed a large diversity of planetary systems (e.g. Mayor et al., 2011; Borucki et al., 2011; Batalha et al., 2013), including a large number of planets with orbital eccentricities much higher than those found in the solar system. These eccentricities challenged our understanding of planet formation occurring in a disk, where eccentricities are dampened by planet-disk interactions. Several authors have suggested that these high eccentricities might have a dynamical origin (e.g. Ford et al., 2001; Jurić & Tremaine, 2008). In this view, the eccentricities arise from planet-planet interactions which can cause a planetary system to become dynamically unstable. During a dynamical instability, close encounters between planets lead to either collisions or planet-planet scatterings. This dynamically active period ends only with a collision, or the ejection of one of the planets from the system. Jurić & Tremaine (2008) showed that this type of process can reproduce the observed distribution of exoplanet eccentricities (at least those above $e=0.2$) for a wide range of possible initial conditions. A similar result was obtained by Raymond et al. (2009b); these authors also showed how the eccentricities can be subsequently dampened by a planetesimal belt for less massive planets ($m_{1}+m_{2}<M_{J}$).

Several authors have noticed a mass-eccentricity correlation among exoplanets: low-mass planets tend to be on less eccentric orbits than high-mass planets (e.g. Ribas & Miralda-Escudé, 2007; Ford & Rasio, 2008; Raymond et al., 2010). Raymond et al. (2010) proposed that planet masses in the same system are correlated. Specifically, systems that form high-mass planets also tend to form equal-mass planets, while low-mass planet systems produce a wide diversity of mass ratios. The present work is an extension of this idea. We use the present-day eccentricities of observed exoplanets to estimate the masses of the planets that were ejected or collided. The result of this process is the exoplanet initial mass function (IMF).

This paper is organized as follows. In Section 2 we summarize the theory behind dynamical stability and discuss previous work. In Section 3 we present an introduction to Approximate Bayesian Computation, and describe how the algorithm can be adapted to the problem we want to solve. We describe our core simulations in Section 4, including our data selection method. The results of these simulations are presented in Section 5. In Section 6 we discuss possible caveats and suggest opportunities for future research. Finally, we summarize and conclude in Section 7.

2 Dynamical instability

In this section we cover some important background and key results regarding the stability of planetary systems. In principle, the eccentricity of the surviving planet is determined by the need to conserve energy and angular momentum. For example, in the simple case of two planets, on coplanar orbits the energy $E$ and angular momentum $\Lambda$ of each planet is,

[TABLE]

where $M$ is the stellar mass, and $m$, $a$, and $e$ are the mass, semimajor axis, and eccentricity of the planet (e.g. Davies et al., 2014). In practice, this has limited utility because the problem is degenerate, even in the simplest case of a two-planet system. After a dynamical instability, the escaping planet carries an uncertain amount of angular momentum (Ford et al., 2001). Even if the angular momentum of the ejected planet was known, one would still have only two equations for five unknowns (the mass of the ejected planet, and the initial mass and semimajor axis of the two planets). In sections 3 and 4 we describe the statistical approach that we use to tackle this degeneracy.

There are two common definitions of dynamical stability. In Hill stability, the ordering of the planets, in terms of their proximity to the central star, is fixed. In other words, the planet orbits never cross, but the outermost planet is allowed to escape the system. The other, more stringent definition, is Lagrange stability. In it, planets retain their ordering (i.e. are Hill stable), remain bound to the star, and variations in semimajor axis and eccentricity remain bounded. While Lagrange stability is the more useful definition, it has proven more difficult to solve mathematically than Hill stability.

There is no analytic solution for the evolution of three gravitating bodies, but in some cases it is possible to derive analytic constraints on the motion of the planets (e.g. Zare, 1977). These constraints can be understood as limitations on angular momentum exchange between planets. For example, Gladman (1993) showed that, to first order, a pair of planets on coplanar orbits are guaranteed to be Hill stable whenever

[TABLE]

where

[TABLE]

and where $M$ is the mass of the star, $m_{k}$ is the mass of a planet, and $a$ and $e$ are the planet semimajor axis and eccentricity in barycentric coordinates. Therefore, for given planet masses and eccentricities, there is a critical value of $\delta$ that guarantees Hill stability (see also the discussion in Barnes & Greenberg, 2006).

Raymond et al. (2009a) refer to the entire left-hand-side of Equation (3) as $\beta$, and the right-hand-side as $\beta_{\rm crit}$. They showed that planet-planet scatterings naturally lead to planetary systems at the edge of dynamical stability, with $\beta/\beta_{\rm crit}$ just above 1, and that this result is in agreement with observation. This is further evidence that dynamical instabilities driven by planet-planet scatterings occur frequently in planetary systems.

For a system of two low-mass planets with near-circular coplanar orbits, the Hill stability criterion (Equation (3)) is approximately $\Delta>2\sqrt{3}$ (Chambers et al., 1996), where $\Delta$ is the semimajor axis separation measured in mutual Hill radii,

[TABLE]

where $m_{1}$, $m_{2}$, $a_{1}$, and $a_{2}$ are the masses and semimajor axes of the two planets. For systems with more than two planets there is no analytic solution, but Chambers et al. (1996) showed that these systems are probably unstable for separations up to $\Delta=10$, at least for equal-mass planets. The time before planets experience their first close encounter ($t_{\rm ce}$) grows exponentially with $\Delta$. For a given value of $\Delta$, $t_{\rm ce}$ seems to depend weakly on the number of planets (for systems with more than two planets) and the planet masses (Chambers et al., 1996; Faber & Quillen, 2007). For a planetary system with ten equal-mass planets, Faber & Quillen (2007) found

[TABLE]

where $\mu=m/M$ is the planet-star mass ratio. This is an extremely steep dependence on $\Delta$, meaning that two systems with similar $\Delta$ values can have very different lifetimes.

3 Approximate Bayesian Computation

We saw in section 2 that it is not possible to exactly determine the mass of an ejected exoplanet given the present-day observables. In this section we explain how one can use modern statistical methods together with N-body simulations to obtain a probability distribution for the ejected planet mass. Our tool of choice is Approximate Bayesian Computation, or ABC (see the introduction by Marin et al., 2011). Like all Bayesian algorithms, ABC is a way to solve the problem

[TABLE]

where $D$ is some data set, $\theta$ is a set of model parameters, and $P(\theta)$ is the Bayesian prior. In its simplest form, the ABC algorithm looks like this

for $i=1:N$ do

     repeat

         $\theta^{\prime}\leftarrow P(\theta)$

         $D^{\prime}\leftarrow{\rm model}(\theta^{\prime})$

     until $\rho(D,D^{\prime})<\epsilon$

     $\theta_{i}\leftarrow\theta^{\prime}$

end for

The function $\rho$ is a distance measure that compares the synthetic dataset $D^{\prime}$ and the observed data $D$, and $\epsilon$ is some critical value. The resulting collection of values $\{\theta_{i}\}$ approximately follows the distribution $P(\theta|D)$.

3.1 Application: two-planet instability

The simplest problem is a two-planet instability that always results in an ejection (i.e. planets are not allowed to collide). In this case, the data $D=(m,e)$ is the mass and final eccentricity of the observed planet, and the model parameter $\theta=q$ is the mass ratio of the ejected and remaining planets ($q=m_{\rm ej}/m_{\rm rem}$). We choose a uniform prior for $P(q)$ (we justify this choice in 4.2). The ABC algorithm becomes

for $i=1:N$ do

     repeat

         $q^{\prime}\leftarrow P(q)$

         $e^{\prime}\leftarrow{\rm simulation}(q^{\prime})$

     until $|e-e^{\prime}|<\epsilon$

     $q_{i}\leftarrow q^{\prime}$

end for

In this investigation we choose $N=1000$. Here we have overlooked the fact that the result of the N-body simulation depends on the initial separation between the two planets. Our approach is to choose a range of orbital separations that are consistent with planet formation: we select only the runs that remain stable for at least 0.5 Myr (see section 4.4). Algorithmically, we can write this as

for $i=1:N$ do

     repeat

         $q^{\prime}\leftarrow P(q)$

         $\Delta^{\prime}\leftarrow P(\Delta)$

         $t^{\prime},e^{\prime}\leftarrow{\rm simulation}(q^{\prime},\Delta^{\prime})$

     until $(|e-e^{\prime}|<\epsilon)\;\;\&\;\;(t^{\prime}>0.5\;{\rm Myr})$

     $q_{i}\leftarrow q^{\prime}$

end for

where $t^{\prime}$ is the time to the ejection. The final step is to allow collisions between planets and record whether or not a collision occurred. Let $C$ be a boolean value that is true when the simulation ends in a collision. The final algorithm becomes

for $i=1:N$ do

     repeat

         $q^{\prime}\leftarrow P(q)$

         $\Delta^{\prime}\leftarrow P(\Delta)$

         $C^{\prime},t^{\prime},e^{\prime}\leftarrow{\rm simulation}(q^{\prime},\Delta^{\prime})$

     until $(|e-e^{\prime}|<\epsilon)\;\;\&\;\;(t^{\prime}>0.5\;{\rm Myr})$

     $q_{i}\leftarrow q^{\prime}$

     $C_{i}\leftarrow C^{\prime}$

end for

Given the list of $\{(q_{i},C_{i})\}$ one can construct a probability distribution of the masses of the two planets:

for $i=1:N$ do

     if $C_{i}^{\prime}$ then

         $m_{1,i}\leftarrow m/(1+q_{i}^{\prime})$

     else

         $m_{1,i}\leftarrow m$

     end if

     $m_{2,i}\leftarrow q_{i}*m_{1,i}$

end for

The set of values $\{m_{1,i},m_{2,i}\}$ approximately follow the posterior distribution $P(m_{1},m_{2}|m,e)$. Finally, we repeat this process for every observed exoplanet that has a measured mass, and for each exoplanet we obtain a set $\{m_{1,i},m_{2,i}\}$. Because all the sets have the same size ($N=1000$), it is straight forward to combine them together. The result is an estimate of the exoplanet IMF.

In section 4 we explain our simulations in greater detail, and we give a rationale for choosing a uniform prior for $P(q)$. We also explain how the stability requirement is implemented in practice.

4 Our simulations

In this section we describe how we selected our list of exoplanets, and we show that exoplanet observations support the use of a uniform prior $P(q)$ as described in section 3. We also explain our simulations in detail, and explain how exoplanet observations inform our choices of initial conditions. At the end of the section we explain how the stability criterion is implemented in detail.

4.1 Data selection

We took the exoplanet catalogue from exoplanets.eu on 16 August 2016. We selected the planets that have mass measurements (i.e. we did not use any mass-radius relation). This means that our sample is dominated by planets discovered by the Radial Velocity method, meaning that we only have the planet’s minimum mass. We then exclude planets with orbital periods shorter than 10 days as many of those may have experienced high-eccentricity migration and tidal circularization. Our resulting dataset included six pulsar planets, so we removed those as well. This left us with 553 exoplanets. Figure 1 shows the semimajor axes and eccentricities of these planets; Table 1 also gives an overview of the dataset.

In this investigation we limit the dynamical simulations to giant planets. To define “giant exoplanet” we use the planet-star mass ratio $m/M_{\star}$, which we normalize with the Jupiter-Sun mass ratio for convenience

[TABLE]

where $M_{\rm J}$ is the mass of Jupiter. We define “giant exoplanet” as those with $0.05<\mu<10$. As a point of reference, Saturn has $\mu=0.30$, Neptune has $\mu=0.054$, and a super-Earth with $m=10M_{\oplus}$ orbiting the Sun would have $\mu=0.031$. Table 2 lists the seven planets that were discovered by a method other than radial velocity or transit.

4.2 Prior distribution of $q$ / mass of planet 2

Our goal is to use an informed prior $P(q)$. Our dataset includes 171 planets with giant planet companions. Figure 2 shows the cumulative distribution of the mass ratios of neighbouring planets in our dataset. The Kolmogorov-Smirnov test cannot distinguish this distribution from the uniform distribution (p-value of 0.96). Therefore, a uniform Bayesian prior is a reasonable choice. We choose 10 values for $q$ spread uniformly between 0 and 1,

[TABLE]

4.3 Mass and semimajor axis of planet 1

The collision cross section between two planets is

[TABLE]

where $R$ is the radius of the planet, $v_{\rm esc}$ is the escape speed at the planet surface, and $v_{\infty}$ is the relative speed of the two planets before gravitational focusing. The focusing factor $\Theta=v_{\rm esc}^{2}/v_{\infty}^{2}$ is commonly known as the Safronov number. In the case of a dynamical instability, where planets approach each other at speeds comparable to the orbital speed, the Safronov number becomes

[TABLE]

When $\Theta\ll 1$, close encounters between planets are are likely to lead to collisions, and when $\Theta\gg 1$ ejections are more common (e.g. Safronov, 1972; Binney & Tremaine, 2008; Ford & Rasio, 2008). The difference is important because collisions are more likely to leave planets in less eccentric orbits than ejections, since they conserve the total angular momentum in the system while ejected planets carry angular (e.g. Ford & Rasio, 2008; Matsumura et al., 2013). We would like our simulations to have planet masses (in terms of $\mu$) and Safronov numbers comparable to the observed population, so that our runs have the correct proportion of collision and ejection events. Unfortunately, most planets in our sample do not have measured radii. Therefore, we make the simplifying assumption that most planets have a density similar to Jupiter, and introduce the quantity

[TABLE]

With this definition, a Jupiter-mass planet orbiting a Sun-like star at 1 AU has $\theta=1$. As a point of reference, table 3 compares $\Theta$ and $\theta$ for the giant planets in the solar system. The two quantities are broadly comparable, but $\theta$ can be calculated for the planets in our sample.

Figure 3 shows the distribution of $\mu$ and $\theta$ for the exoplanets in our dataset. Clearly there is a wide range of values for both $\mu$ and $\theta$. For this reason, we decided to do three different sets of runs (calling the sets A, B, and C) with the mass and semimajor axis of planet 1 spread along the best-fit line on figure 3. Planet 1 is the more massive planet in each run, and is the one more likely to survive. Our choices for $m_{1}$, and $a_{1}$ — shown in table 4 — are a compromise between covering the full range of observed $\mu$ and $\theta$, while also having our runs concentrated in the regions of the parameter space that have more planets. As it turns out, 56% of the planets in our sample are closer to set B (on a log scale). Therefore, we perform most of our runs in set B (see table 4).

4.4 Other orbital elements

We have already covered how we select $m_{1}$, $a_{1}$, and $m_{2}$. The remaining orbital elements are chosen to make our simulations consistent with planet formation. Giant planets form inside a protoplanetary disk, with typical lifetimes in the order of $\sim$ 3-6 Myr (e.g. Hartmann et al., 1998; Haisch et al., 2001; Mamajek, 2009). Orbital configurations that lead to orbit crossing on a timescale significantly shorter than 3 Myr would lead to ejections and collisions during the disk phase, with enough time for the disk to subsequently dampen the eccentricity of the remaining planet. For this reason, we require that the system experience no collisions or ejections for at least 0.5 Myr.

To conduct our simulations, we select a range of semimajor axes for planet 2 near the Hill stability limit. The exact range of $a_{2}$ varies depending on the stability criterion. Figure 4 shows the values of $\Delta$ (equation (8)) for set B. The $\Delta$ values are equally spaced, in steps of 0.1. For each $\Delta$, we solve for $a_{2}$ and we run ten N-body simulations using the hybrid integrator in mercury (Chambers, 1999). Both planets are initially in circular orbits, but we give them mutual inclinations of a few degrees (similar to the solar system). We assign each planet a random inclination $I\in[0^{\circ},5^{\circ}]$, random longitude of ascending node $\Omega\in[0^{\circ},360^{\circ}]$, and random mean longitude $\lambda\in[0^{\circ},360^{\circ}]$. Mutual inclinations are important because an overly flat system will have an unrealistic number of collisions.

The simulations initially run for 10 Myr. If there is an ejection or collision in less than 0.5 Myr, we do not use the run. If there are no collisions or ejections, we extend the runs in 20 Myr increments up to 70 Myr. Figure 4 shows a bar plot for each value of $(\Delta,m_{2})$ that indicates the fraction of runs that became unstable in less than 0.5 Myr (red), or 0.5-70 Myr (green), or had no instability at the end of the 70 Myr run. Notice that smaller $m_{2}$ requires larger $\Delta$ to be stable for 0.5 Myr. This is the reason why the range of $a_{2}$ varies with $m_{2}$. In a similar way, changing $m_{1}$ also affects the stability region, so the range of $a_{2}$ that are consistent with planet formation change as well.

5 Result of planet-planet scattering

In this section we present the results of our N-body simulations. In section 5.1 we apply the Bayesian algorithm described in sections 3 and 4 to our two-planet simulations and we construct an estimate for the exoplanet IMF. Then in section 5.2 we present a more limited set of simulations with three giant planets and discuss how a three-planet instability differs from the two-planet case.

5.1 Two-planet systems

Here we present our main results. Figure 5 shows the final eccentricity of the surviving planet for every single run in set B that had an instability after 0.5 Myr. By far the most common result was for the lower mass planet to become ejected from the system. In some cases the two planets collided, and in others it was the inner, more massive planet that was ejected. As a general rule, collisions are associated with low-eccentricity outcomes, even for equal-mass planets. The connection between collisions and dynamically quiet histories has already been noted by other authors in the context of three-planet instabilities (Matsumura et al., 2013; Carrera et al., 2016).

Perhaps the most important feature of Figure 5 is the correlation between the eccentricity of the remaining planet $e$ and the mass ratio of the planet that was ejected or destroyed over the mass of the surviving planet. In most cases $q=m_{2}/m_{1}$ where planet 2 is the outer planet; but in the runs where the inner (more massive) planet is ejected, we write $q=m_{1}/m_{2}$. The correlation between $e$ and $q$ is not surprising, but it is important — it is what allows us to produce a probability distribution of $q$.

Figure 6 makes this point more concrete. Here we choose four sample eccentricities ($e_{\rm fin}\in\{0.2,0.4,0.6,0.8\}$), and for each one we compute a probability distribution $k(q)$ using the data in set B (Figure 5). Conceptually, for any given value of $e$ one can draw a vertical line in Figure 5 and select all the points near that line. We use kernel density estimation (KDE) to convert these points into a smooth distribution. A KDE is conceptually equivalent to replacing each point with a Gaussian distribution with standard deviation $h$ (called the bandwidth). In Figure 6 we chose a wide bin ($e_{\rm fin}\pm 0.1$) to focus on the overall shape of the curve, but notice that in Figure 5 there is a relatively small number of points for any given $e_{\rm fin}$. As the number of simulations continues to increase, our statistics will gradually improve. Finally, to produce the exoplanet planet IMF, we generalize the procedure:

• •

For planets with $\mu\geq 0.05$ and $e\geq 0.05$ we randomly select $N=1000$ runs (see section 3) with $|e_{\rm sim}-e_{\rm obs}|<0.02$ and compute the initial planet masses $m_{1}$ and $m_{2}$.

• •

For planets with with $\mu<0.05$ or $e<0.05$ we assume that there was no planet-planet scatter, and we simply copy the planet mass $N=1000$ times.

Together, these sets give a synthetic population that approximately follows the initial mass function of the exoplanets in our dataset. As in Figure 6, we compute the kernel density $k(\mu)$. Figure 7 shows the final result for all the planets in our dataset, and for the planets orbiting Sun-like stars. One interesting feature of these plots is that the observed exoplanet population appears to have a “valley” at around $\mu=0.2$ (sub-Saturn planets), but the synthetic IMF suggests that this valley is not primordial, but was carved out later by planet-planet scatterings. If this interpretation is correct, there should be a population of free-floating Saturn-like planets that were ejected from their host system by a more massive planet. However, note that Veras & Raymond (2012) have shown that dynamical instabilities cannot be the primary source of the observed population of free-floating giant planets (Clanton & Gaudi, 2017), as that would require an implausible number of planetary ejections per system.

5.2 Three-planet systems

Our results so far only apply to two-planet systems, and it is important to understand to what extent they generalize to multiple-planet systems. In this section we present a set of N-body simulations with three giant planets, and we discuss the similarities and differences between two-planet and three-planet systems.

All our three-planet simulations have a Jupiter-mass planet at 1 au, and two exterior planets with half the mass of Jupiter. We call this set J+J/2+J/2. The planets have a uniform separation in terms of their mutual Hill radii (equal $\Delta$, see Equation (8)). We compare these runs against the two-planet runs with $q=1$ (we call them J+J) and $q=0.5$ (we call them J+J/2). The three-planet runs are generally less stable, and we need to be more widely spaced ($\Delta>4$) for them to last 0.5 Myr. That said, in terms of orbital energy and angular momentum the difference is not very large. For example, in a J+J/2+J/2 system with $\Delta=4.1$ the two outer planets have 79% of the binding energy and 113% of the angular momentum of the outer planet in a J+J system with $\Delta=2.5$.

Figure 8 shows the final eccentricity of the surviving planets in each set of runs. Broadly speaking, J+J/2+J/2 produces eccentricities similar to those in J+J/2, and somewhat lower than those in J+J. In J+J/2+J/2 there are a few runs where only one planet survived. These planets are more eccentric, and are generally consistent with J+J.

6 Some caveats

In this section we discuss sources of bias or uncertainty that might affect the accuracy of our results. When possible, we discuss ways that these sources of error might be corrected or at least measured.

6.1 Observational errors and biases

The exoplanet sample used in this investigation is plagued with observational biases. There are completeness issues because massive planets are easier to detect than lower-mass planets. In addition, eccentricities from radial velocity surveys are unreliable and the RV signal caused by two planets in circular orbits can be difficult to distinguish from the signal caused by a single planet in an eccentric orbit. De-biasing RV surveys is a difficult problem that lies beyond the scope of the present investigation.

We did explore the completeness issue with a reduced subset of observed planets. We selected planets in the range of periods 10 to 300 days and masses between 0.1 and 10 $M_{\rm J}$, and repeated our procedure. This subset suffers from small number statistics, as there are only 27 planets in the set; nonetheless, our key conclusion continues to hold for this data set — that the observed exoplanet population has a deficiency in sub-Saturn-mass planets, and that that this gap is not primordial, but is the result of ejections.

6.2 System architecture

This investigation is limited to planetary systems with two dynamically dominant planets. In section 5.2 we quantified some of the key differences between 2-planet and 3-planet systems. That said, a full investigation is not feasible because the problem is too degenerate when there are more than two giant planets in the system.

In all our simulations the inner planet is the most massive. It is possible that if the order of the planets was reversed, the simulation results would be different. To test this, we conducted 550 simulations with the same parameters as J+J/2 (i.e. set B with $q=0.5$), but we reversed the order of the planets, so that the inner plant has a mass of 0.5 $M_{\rm J}$ and the outer planet has 1 $M_{\rm J}$. We call this new set J/2+J. The two sets of runs give similar results. The eccentricity distributions are similar, and in both cases it is always the less massive planet that becomes ejected. One difference is that J/2+J is more collisional than J+J/2. In J+J/2 there are 2.5 ejections per collision, while in J/2+J there are only 1.3 ejections per collision. Figure 9 gives an overview of the simulation results. The two sets of runs have the same mean eccentricity, but J/2+J may have a wider eccentricity distribution.

6.3 Eccentricity damping

The next complication is that gravitational interactions with a planetesimal belt can dampen a planet’s eccentricity (Raymond et al., 2010). In fact, some authors suggest that this might have occurred in the solar system at the time of the Late-Heavy Bombardment (Gomes et al., 2005). In this way, a planetesimal belt may erase information about the dynamical history of the planetary system. Fortunately, this is mainly a concern for low-mass (sub-Saturn) planets in wide orbits (Raymond et al., 2010). Hence, it is unlikely that this type of damping has a major impact in our results. Nonetheless, it is desirable to quantify how often this occurs. Raymond et al. (2010) suggested a search for $K\approx 5{\rm m\,s}^{-1}$ planets with long periods ($P\sim 10$ yr). Another possibility is to look for planetary systems with highly depleted debris belts, which is an independent indicator of past dynamical instability (Raymond et al., 2011, 2012) and measure how often these systems have low eccentricities (which could suggest eccentricity damping).

6.4 Secular effects

Finally, in a multiple planet system, planets can exchange angular momentum over secular timescales. For example, if a planetary system with three giant planets ejects one, the two remaining planets will exchange eccentricity periodically. In other words, some of the eccentricities that we observe today may differ from those at the end of the dynamical instability. To tackle this one could model the secular evolution of the remaining planets and note the range of eccentricities that the planets acquire. This, in turn, gives a range of possible solutions for the ejected planet mass. We will revisit this idea in future work.

7 Summary and conclusions

In this work we developed a novel Bayesian framework, supported by N-body simulations, to estimate the probable masses of planets that have suffered collisions or ejections, using the present-day masses and orbits of the surviving planets. When applied to the entire exoplanet population, this technique yields the exoplanet initial mass function. We then demonstrated the use of this technique using present-day observations and several thousand N-body integrations.

The observed exoplanet population has a paucity of sub-Saturn-mass planets. We find that this gap is not primordial, and is mainly the result of dynamical instabilities where more massive giant planets eject less massive ones. This in turn implies that there is a population of free-floating sub-Saturn planets that might be detectable by the upcoming WFIRST telescope through its micro-lensing survey (Clanton & Gaudi, 2017).

Acknowledgements

The authors are supported by the project grant “IMPACT” from the Knut and Alice Wallenberg Foundation (KAW 2014.0017), as well as the Swedish Research Council (grants 2011-3991 and 2014-5775), and the European Research Council Starting Grant 278675-PEBBLE2PLANET. D.C. acknowledges Dimitri Veras for helpful discussions and guidance in calculating the Hill stability limit. Computer simulations were performed using resources provided by the Swedish National Infrastructure for Computing (SNIC) at the Lunarc Center for Scientific and Technical Computing at Lund University. Some simulation hardware was purchased with grants from the Royal Physiographic Society of Lund.