Stellar Flybys Interrupting Planet-Planet Scattering Generates Oort Planets

TL;DR

This paper proposes a new mechanism involving stellar flybys that can produce distant 'Oort' planets through planet-planet scattering interruption, supported by numerical simulations and efficiency analysis.

Contribution

It introduces a novel process where stellar flybys halt planet scattering, leading to the formation of distant Oort-like planets, expanding understanding of wide-orbit exoplanets.

Findings

Less than 1% of gas giants become Oort planets via this mechanism.

Up to a few percent of Neptunes and super-Earths may become Oort planets.

Numerical simulations support the efficiency of stellar flybys in creating Oort planets.

Abstract

Wide-orbit exoplanets are starting to be detected, and planetary formation models are under development to understand their properties. We propose a population of "Oort" planets around other stars, forming by a mechanism analogous to how the Solar System's Oort cloud of comets was populated. Gravitational scattering among planets is inferred from the eccentricity distribution of gas-giant exoplanets measured by the Doppler technique. This scattering is thought to commence while the protoplanetary disk is dissipating, $10^6-10^7$ yr after formation of the star, or perhaps soon thereafter, when the majority of stars are expected to be part of a natal cluster. Previous calculations of planet-planet scattering around isolated stars have one or more planets spending $10^4-10^7$ yr at distances >100 AU before ultimately being ejected. During that time, a close flyby of another star in the cluster may dynamically lift the periastron of the planet, ending further scattering with the inner planets. We present numerical simulations demonstrating this mechanism as well as an analysis of the efficiency. We estimate an occurrence of planets between 100 and 5000 AU by this mechanism to be <1% for gas giants and up to a few percent for Neptunes and super-Earths.

Tables

Table 1. Table 1: Initial Conditions for Sample Simulation (run9)

Planet Name	Semi-Major Axis	Mass	Radius	Eccentricity	Inclination	Ascending Node	Argument of Pericenter	Mean Anomaly
	$AU$	$M_{Jupiter}$	$R_{Jupiter}$		$deg$	$deg$	$deg$	$deg$
9A	1.2563	0.6544	1.2229	.0045	5.2181	9.1809	280.4906	171.6866
9B	1.7460	0.1668	0.9838	.0054	7.9660	229.1612	43.2556	70.8627
9C	2.1131	0.1298	0.4769	.0083	9.7999	220.6738	310.6611	102.6931
9D	3.3368	2.1931	1.0545	.0145	0.8957	239.1286	39.4472	176.9094
9E	5.3686	0.1458	1.3875	.0153	5.6120	238.2865	186.4250	32.6289

Table 2. Table 2: Rate of Outcomes for Strong Encounters

Outcome	Occurrence
Lifts Periastron	20%
Rebinds	7%
Unbound	43%
Reduces Periastron	17%
No Flyby Effect	13%

Equations

$$f_{saved,tot}(a)=N(b_{min}<10a)p(a)f_{strong}{f}_{save},$$

$$\Gamma =n\pi b_{max}^{2}v,$$

$$p(b_{min})=\frac{2N{b}_{min}}{b_{max}^2}(1-\frac{b_{min}^2}{b_{max}^2}{)}^{N-1}.$$

$$\frac{M_P}{M_{\oplus }}\gtrsim {\left(\frac{M_{\ast }}{M_{\odot }}\right)}^{3/4}{\left(\frac{t_{\ast }}{10^{9}yr}\right)}^{-1/2}{\left(\frac{a_P}{1AU}\right)}^{3/4}.$$

$$n_o{\tau }_o=\tau (n_o+\frac{1}{2}\dot{n}\tau )$$

$$\tau =\frac{-n_o\pm \sqrt{n_o^2+2n_o\dot{n}{\tau }_o}}{\dot{n}},$$

$$p(x_o>x)=\prod _{i=1}^{N}p(x_i>x)=(1-x{)}^N.$$

$$p(x_o<x)=1-(1-x{)}^N.$$

$$p(x)=\left|\frac{dp(x_o<x)}{dx}\right|=N(1-x{)}^{N-1}$$

$$p(x_i>x)={\int }_x^{1}f(x^{\prime })d{x}^{\prime }=1-x^2$$

$$p(x_o<x)=1-\prod _{i=1}^{N}p(x_i>x)=1-(1-x^2{)}^N$$

$$p(x)=2xN(1-x^2{)}^{N-1}.$$

$${\int }_0^{b_{max}}A2xN(b_{max}^2-x^2{)}^{N-1}dx=1.$$

$$p(x)=2b_{max}^{-2}xN(1-\frac{x^2}{b_{max}^2}{)}^{N-1}.$$

$$p(b_{min})=\frac{2b_{min}}{d}(1-\frac{b_{min}^2}{dN}{)}^{N-1}$$

$$\underset{n\to \infty }{\lim }(1+\frac{x}{n}{)}^n=e^x,$$

$$p(b_{min})=\frac{2b_{min}}{d}{e}^{-\frac{b_{min}^2}{d}}.$$

$$f(x;\sigma )=\frac{x}{{\sigma }^2}{e}^{-x^2/(2{\sigma }^2)},$$

Full text

celestial mechanics, planetary systems, planets and satellites: dynamical evolution and stability

††journal: AJ††software: REBOUND, Forecaster

1 Introduction

The outer parts of planetary systems are still mysterious. Myth and speculation shroud the outer parts of our own system, both in the popular imagination and in scientific circles. From Nirabu, to a companion triggering comet showers (Whitmire & Jackson, 1984), to Planet X perturbing the outer planets (Seidelmann & Harrington 1988, but see Standish 1993), the unknown is intriguing.

Around the Sun, searches by optical light have come up empty (Schwamb et al., 2010; Brown et al., 2015). Gravitational probes have been developed (Gaudi & Bloom, 2005; Zakamska & Tremaine, 2005; Fienga et al., 2016), but have found no convincing evidence. The latest motivation for such searches, the clustering properties of high-perihelion objects (Trujillo & Sheppard, 2014; Batygin & Brown, 2016) or collective inclination of the Kuiper Belt (Volk & Malhotra, 2017), have touched off renewed efforts to find such planets. Theories of migration of the gas giants suggest one or more ice giants could have been ejected (Batygin et al., 2012; Nesvorný, 2018), although maybe a planet was instead placed in the Oort cloud.

In other planetary systems, we have good limits on the frequency of gas giant planets from direct imaging from 30 AU to about 300 AU, with an overall occurrence rate of 5-13 $M_{J}$ companions of $0.6^{+0.7}_{-0.5}\%$. For FGK stars, there is an upper limit occurrence rate of $5.1\%$ for companions between 100 AU and 1000 AU (Bowler, 2016). This range to 1000 AU was also studied via a survey with Spitzer Space Telescope (Durkan et al., 2016) that imaged the space around white dwarfs and found, “assuming a mass distribution of ${dn}/{dm}\propto{m}^{-1.31}$, we constrain (at 95% confidence) a population of 0.5-13 $M_{J}$ planets at separations of 100-1000 AU with an upper frequency limit of 9%.” The Wide-orbit Exoplanet search with InfraRed Direct imaging (WEIRD) survey (Baron et al., 2018) extended this range out to 5000 AU, finding a constraint on “the occurrence of $1-13M_{J}$ planetary-mass companions on orbits with a semi-major axis between 1000 and 5000 AU at less than 0.03, with a 95% confidence level.” While direct imaging surveys provide mostly upper limit constraints on wide-orbit companions, the NASA Exoplanet Archive lists 10 confirmed planets beyond 100 AU with mass $<$13 $M_{J}$ (as of April 29, 2019), implying a non-zero lower limit to wide-orbit planets.

Theoretically, we expect planets to be able to form in situ only out to $\sim$30 AU by core accretion (Murray-Clay, 2010). Formation could perhaps occur beyond 100 AU out to several hundred AU by disk instability (limited mainly by the observed sizes of protoplanetary disks; Andrews et al. 2009). However, such planets would form in high-mass disks and would be quite massive ($\sim$10 MJ; Rafikov 2005) and are likely to migrate rapidly inward (Baruteau et al., 2011). An analysis involving photoevaporation of the disk showed that migration may be directed outward (Veras & Armitage, 2004), but even this mechanism would not populate the region beyond several hundred AU.

To our knowledge, only one mechanism has been proposed that can place planets on stable orbits of $\sim$1000 AU: a free-floating planet may bind to a star as the natal cluster of stars dissolves (Perets & Kouwenhoven, 2012), at an efficiency ranging from 1% to 9%. We wondered whether a new mechanism, in which the planet stays with its parent star, could populate the same region.

1.1 Analogy to the Formation of Cometary Orbits

The formation mechanism we study here is closely analogous to the leading mechanism thought to populate the Solar System’s outer reaches with comets. Oort (1950) realized that bodies could transition to very distant orbits by first scattering off the known planets, which raises their aphelion to thousands of AU while the perihelion remains among the planets, and in such a state the passage of stars in the Galaxy can raise the perihelion up out of the planetary region. Then the objects would be put in “cold storage” and occasionally return later, warming up and sprouting a tail as a comet.

Simulations of the formation of the Solar System’s Oort cloud have shown a relatively inefficient process, beginning with an initial rapid population by Jupiter and Saturn and followed by a more gradual period of growth mostly due to Uranus and Neptune, where the fraction of bodies making it to the Oort cloud ranges between 5% and 7.6% (Dones et al., 2004).

Tremaine (1993) extended this idea to planetary perturbers other than the Solar System’s giant planets, inspired by the first discoveries of exoplanets around a pulsar. He leveraged a different external torquing mechanism, that of the Galactic tidal field (Heisler & Tremaine, 1986), which makes the problem especially well-posed. He found that the slower scattering of planetesimals off Neptune-mass planets most efficiently populates the clouds around other stars (called exo-Oort clouds by Veras et al. 2014; Stone et al. 2015).

The most comprehensive model of the Solar System Oort cloud, including the effects of an early cluster environment, has been made by Brasser et al. (2006). They calculated how objects that are being scattered out by planets can be saved from ejection by passing stars. This Oort mechanism can produce objects on orbits similar to Solar System objects such as Sedna (e.g. Morbidelli & Levison 2004, mechanism 4). In this paper we follow some of the assumptions and methods of the simpler and earlier model of Fernández & Brunini (2000).

A population of high-apastron planets would be needed in order to emplace a wide-orbit planet via this mechanism. Planet-planet scattering has been shown to explain most of the eccentric gas giant planets (e.g., Jurić & Tremaine 2008; Chatterjee et al. 2008). This would give a source of planets on their way to be ejected, that can become stranded in the Oort cloud by passing stars. Veras et al. (2009) have computed the distribution of long-period planets as a function of time using this concept, and they find $>$10% of systems undergoing scattering have an outermost planet beyond 100 AU at a time of 5 Myr. This population of planets then decreases over time as the planets are ejected from the system; we propose instead a mechanism to save these planets from ejection and maintain a wide orbit for Gyr timescales. Veras et al. (2009) mention the possibility of giant planets being perturbed onto an Oort-like orbit, but due to an expected low efficiency as well as the difficulties in observing such old and cold giant planets via direct imaging, their work focuses on young giant planets in the process of being scattered and eventually ejected.

1.2 Cluster environment

By comparing the counts of young clusters and the rate at which stars in the field of the galaxy were formed, Lada & Lada (2003) concluded “embedded clusters may account for a large fraction of all star formation occurring locally.” If these clusters survive the embedded phase, close encounters with other stars are required to cause individual stars to evaporate from the cluster and reach the field of the galaxy. Usually, it is around stars in the field of the Galaxy, not in star clusters, that exoplanets are found. However, some clusters may dynamically disrupt via prompt evaporation of the gas owing to photoevaporation by the most massive stars, which would make the cluster phase brief and cut down the number of encounters. Some investigators recently found that Gaia parallaxes show the Scorpius-Centaurus OB association was never highly clustered (Wright & Mamajek, 2018; Ward & Kruijssen, 2018). Indeed, statistical clustering of star-forming regions need not imply gravitational binding (Elmegreen, 2018). Nevertheless, the main paradigm remains that most field stars, and hence most exoplanet hosts, formed in clusters and therefore endured flybys of other stars early in the system history (Adams, 2010).

1.3 Plan for this paper

We describe numerical methods by which we explore the mechanism in section 2. In section 3 we show examples of the mechanism, discuss its dependence on planetary and cluster parameters, and estimate its efficiency. Finally, in section 4, we discuss the relation to other mechanisms and propose observational tests of the theory.

2 Methods

Our calculations follow those in the planet-planet scattering literature – a distribution of planets will begin scattering and place some on long-period orbits. Our treatment of what happens next diverges from the literature. Instead of placing an outer boundary (500 AU is common), we will allow planets to remain in the simulation and simulate passing stars. The effect of these flybys could raise the scattering planet’s periastron through torques, saving it from further scattering, an effect which would have been missed in previous calculations. Contrariwise, the flyby could tidally pull the planet away, in which case it would be a true ejection, taken away from the planetary system.

2.1 Initial Conditions

The distributions of initial conditions were chosen in a somewhat arbitrary fashion to provide a wide array of potential planet-planet scattering encounters. The goal is not to have a random sampling of systems as close to the distribution in nature as possible; the goal is to produce a few examples of models in which the planet-saving mechanism occurs, to learn how it operates. The starting point for the simulations was assumed to be after the dissipation of the gaseous protoplanetary disk, once gas giant planets have accreted their massive envelopes.

Planetary masses were chosen uniformly in log $M$ between $M$ = 0.1 and 10 $M_{J}$. The semi-major axis for the innermost planet was chosen uniformly in log $a$ between log $a$ [AU] = 0.0 and 0.1. The remaining semi-major axes were assigned at approximately 5 mutual Hill radii separation (Chambers et al., 1996). If the outermost semi-major axis was beyond 10 AU, a new distribution of masses was drawn and the semi-major axis assignment repeated. A distribution of initial semi-major axes is shown in Figure 1.

The inclination and eccentricity for each planet were drawn from a Rayleigh distribution with scale parameters of 0.1 radians and 0.01, respectively. The ascending node, argument of pericenter, and mean anomaly were all chosen uniformly between 0 and $2\pi$ radians.

The radius of each planet was assigned using the mass-radius prediction function from the Forecaster package (Chen & Kipping, 2017), which forecasts a radius using a random draw from a probabilistic mass-radius relation. Both the central star and flyby star were taken to have solar mass, and the central star solar radius. Planet-planet and planet-central star collisions were allowed during the simulation.

2.2 Simulation Methods

The simulations were computed using the REBOUND package (Rein & Liu, 2012) with initial conditions as described in section 2.1. Computations were performed using the IAS15 integrator (Rein & Spiegel, 2015), a 15th-order adaptive-step integrator that can handle close encounters and highly eccentric orbits, checking for ejections (distance $>$2 x 105 AU) every 0.5 years. Positions, velocities, and astrocentric orbital parameters were recorded every 5 years through the duration of the simulation. Initial runs were computed over 1 Myr while allowing for flyby stars.

Flyby stars were inserted into the simulation in the event that the planets were scattering, as determined by having a planet with a semi-major axis of greater than 20 AU. Each star’s initial location was randomly placed on a 1000 AU sphere centered on the central star with an inward-pointed velocity vector randomly drawn with a one-dimensional dispersion of 1 km s*-1* and probabilistically weighted with the cosine of the incidence angle via rejection sampling (Henon, 1972). The flyby star was removed when it again passed 1000 AU, and another star was sent flying by (if the scattering criterion, any planet’s $a>20$ AU, was met) after a delay of 10 kyr, a timescale that allows any violent instabilities in the system to resolve.

There were ten separate initial runs, which were the parent simulations. The subset of each parent simulation during the passage of a flyby star is a child simulation, which can be considered as independent events due to the time delay between insertions. There were 230 total child simulations. The planetary properties of mass and radius were the same for all the child simulations of a given parent simulation; however, the orbital properties of the planets changed for each child simulation due to the interactions occurring between flyby star passages.

After completion of the initial parent simulation, the child simulation results were examined for potential stranded wide-orbit planets due to a stellar passage. If such a result was seen, a new simulation was initiated from the point of the star’s removal and computed without further flybys over longer timescales (on the order of 10-100 Myr), recording every 50 years, to determine the stability of the post-flyby configuration of the system, particularly looking at the variability of the semi-major axis of the Oort planet.

3 Results

3.1 Example Simulation

Initial conditions for one example of the simulations described in section 2.2 are shown in Table 1 and the results shown in Figure 2. The left panels (2a-c) show the first 30 kyr of the initial run, where the inserted stars are shown in black in 2a. While the first flyby star (undeflected $b=263$ AU) has a deep, strongly gravitationally focused encounter that ejects planet 9C (green) from the system entirely, the second flyby star (undeflected $b=827$ AU) appears to lift the periastron of the scattering planet 9B (orange) sufficiently to prevent further scattering and leave the planet on a stable orbit at a semi-major axis of approximately 300 AU. The right panels (2d-f) show the addition of the second simulation over 10 Myr of the system. Planet 9B appears to be stable on its wide orbit resulting from the stellar passage.

After the stellar passages, an additional planet, 9E (purple) undergoes significant orbital excitation from the inner planets. It would be ejected on a shorter timescale than observed, were it not for the intervention of the now-Oort planet, 9B, which lifts its periastron at around 2 Myr in the simulation. The two planets apparently enter a stage of quasi-stability. Planet 9E has an orbital energy closer to zero, so any strong encounters between it and 9B will tend to eject 9E, leaving 9B as an Oort planet.

This simulation run is included as an illustration of the mechanism described in section 1.1. Similar results were seen in other simulations.

3.2 Efficiency

The expected efficiency of saving a planet on a wide orbit of $a$ is given by

[TABLE]

where $N(b_{min}<10a)$ is the expected number of relevant stellar passages with an impact parameter within $10a$, $p(a)$ is the probability of a planet being at a given semi-major axis, $f_{strong}$ is the fraction of strong encounters for $b/a<10$, and $f_{save}$ is the fraction of strong encounters that result in a save. Each of these components will be defined and discussed in detail in the following sections.

3.2.1 Stellar Passage Frequency

In order to understand the feasibility and efficiency of this method of producing stable, wide-orbit planets, the frequency of close stellar passages must be evaluated.

Following the lead of Fernández & Brunini (2000), the natal star cluster is taken to have an initial number density of 25 stars pc*-3* for stars of approximately 1 $M_{\odot}$ and characteristic relative velocities of 1 km s*-1*. The cluster is assumed to dissipate over a timescale of 108 years with the density decreasing linearly to zero over the lifetime of the cluster. This assumption neglects the effects of gas in the cluster. Having gas in the cluster for longer reduces the dissipation rate of the cluster and increases the number of close stellar encounters (Proszkow & Adams, 2009), making this a conservative assumption. Further, the average stellar mass in nearby clusters is 0.5 $M_{\odot}$ (Proszkow & Adams, 2009), so our assumption that all the stars are 1 $M_{\odot}$ reduces the number of stars in the cluster and thereby reduces the interaction rate.

In addition to the relatively “dense” cluster (actually less dense than many stellar clusters), a “loose” cluster with initial density 10 stars pc*-3* was also considered, with the same dissipation profile.

The rate of flyby stars within a specified impact parameter, $b_{max}$, is given by

[TABLE]

where $n$ is the number density of stars in the cluster and $v$ is the incoming relative velocity. The effect of gravitational focusing is neglected. For $v_{\infty}$ = 1 km s*-1*, this approximation is appropriate beyond $\sim$1000 AU.

For statistical modeling, we calculate the passages with an impact parameter of less than 5 x 104 AU over the life of the cluster. A sample draw is a choice of $(b/b_{max})^{2}$ from a uniform distribution between 0 and 1. The time between flybys was modeled as a Poisson process at the rate specified in equation 2. Because the density is not constant over time, the wait time between passages was corrected to account for the change in density over the uncorrected wait time (see Appendix for details). The results of one example cluster model simulation are shown in Figure 3. Close passages can disrupt the planetary system (as demonstrated in Figure 2 with the ejection of planet 9C) in addition to “saving” the planets, so the relevant passages that we consider in our analysis are the recent closest passages — that is, the closest passages without any future closer passages.

A distribution of the total number of passages and the single closest impact parameter over 106 simulations is shown for both a dense cluster (25 stars pc*-3*) and a loose cluster (10 stars pc*-3*) in Figure 4. The probability distribution for the total number of flybys is Poisson, which approaches a normal distribution with mean $N=\int\Gamma dt$ and standard deviation $\sqrt{N}$. For the closest impact parameter, $b_{min}$, the probability distribution function is given by

[TABLE]

In the limit as $N\rightarrow\infty$ and $b_{max}\rightarrow\infty$, $b_{min}$ is distributed as a Rayleigh distribution with scale parameter $\sigma=\sqrt{\frac{b_{max}^{2}}{2N}}$, which is plotted in Figure 4, bottom panels. See Appendix for details of derivation.

Both distributions are dependent on the density of the cluster. Analytically, the total number of passages within the specified $b_{max}$ is expected to increase linearly with the initial density of the cluster. The median closest passage is expected to decrease as the square root of the initial density of the cluster. This was consistent with our simulations at varying densities.

The $b_{min}$ distribution in Figure 4 is for the single closest passage during the cluster’s lifetime. Including all recent closest passages, as illustrated in Figure 3, increases the number of flybys. The average cluster, with an initial density of 25 stars pc*-3*, has 5 of these recent closest flybys. Figure 5 shows how the $b_{min}$ distribution differs for the nominal cluster by including all recent closest flybys. While the single closest passage distribution is Rayleigh-like, the distribution for all recent closest passages has a fatter tail and increased occurrences of more distant flybys.

While our statistical approximation is a simplified approach, it is consistent with more detailed stellar cluster simulations. Using dynamical N-body simulations, Proszkow & Adams (2009) find that, depending on cluster parameters, the total rate of close stellar encounters less than 5 x 104 AU ranges between 0.5 and 13.5 Myr*-1* in the first 10 Myr (assuming our 25 stars pc*-3* are in a sphere of radius $\sim$1 pc and therefore N $\approx$ 100). Our approximation gives an average of 4.5 Myr*-1* encounters in the first 10 Myr.

3.2.2 Stellar Passage Impact

A flyby star is most likely to have a strong impact on planets that have a semi-major axis of the same order of magnitude as the flyby’s impact parameter. As can be seen in Figure 6, with only a couple exceptions, the vast majority of encounters resulting in a significant change in periastron of the planet happened for $b/a\lesssim 10$.

Examining these encounters, we can determine the fraction of encounters with $b/a<10$ that result in a strong encounter. A strong encounter is defined to be one where the planet’s periastron changes by more than 20 AU. We find that there were 51 encounters with $b/a<10$ and 20 of these were strong encounters, giving an $f_{strong}$ of 0.39 $\pm$ 0.07.

For strong encounters seen during our simulations, we find that there are five likely outcomes. The planet may have its periastron lifted (generally resulting in being “saved”); an already unbound planet from scattering may become rebound to the host star (also a save, although more dependent on specific star-planet interaction); the planet may become unbound, either due to planet-planet scattering or from an interaction with the flyby star; the planet may have its periastron reduced (potentially putting a previously safe planet into danger of scattering); or the change may be due to planetary interactions where the flyby star appears to have no effect.

Across the simulations performed, there were 230 flyby stars. Of these 230, 30 resulted in a strong encounter. These 30 encounters are categorized as described above, resulting in the branching probabilities shown in Table 2. Combining the two outcomes that save the planet, the simulated $f_{save}$ is 0.27 $\pm$ 0.08.

The uncertainty for $f_{strong}$ and $f_{save}$ is calculated via the normal approximation of the distribution of error for a binomially-distributed observation.

3.2.3 Distribution of Planets During Scattering

In order for a planet to be saved on a wide orbit by a flyby star, there must first be a planet that has been scattered out to a large semi-major axis at the time that the flyby star passes. The expected distribution of semi-major axes over time depends on the properties of the scattering planet; it must be massive enough to scatter another object to a certain distance in the given system age. Tremaine (1993) found that a necessary condition on the interior planet doing the scattering for the formation of an Oort-type comet cloud is

[TABLE]

For our range of planet parameters (0.1 $\leq$ $M_{P}/M_{J}$ $\leq$ 10; 1 $\leq$ $a_{P}$/1 AU $\leq$ 10), this condition is met within the cluster’s lifetime ($t_{*}$ = $10^{8}$ yrs). Considering in our case the scattered bodies are planets and not comets, the scattering planet must also be at least as massive as the scattered planets (Ford & Rasio, 2008).

To derive the distribution of semi-major axes of scattered planets, an N-body simulation was run having one planet with parameters chosen in a similar manner as described in section 2.1 (with N = 1). The semi-major axis of the planet was chosen uniformly in $a$ between 1 and 10 AU. The randomly generated planet had $M$ = 0.59 $M_{J}$ and $a$ = 9.49 AU. 100 test bodies were added with initial semi-major axes drawn uniformly within 10% of the planet’s semi-major axis and other orbital parameters drawn in the same manner as the planet’s. The test bodies were massless in order to decrease computation time. The effect of including the mass of the planets being scattered would change the outcome of any given energy exchange interaction by a factor of $\sim$2, which would not appreciably alter the resulting distribution. The system was integrated for $10^{8}$ years.

In order to examine the dependence of the distribution on the mass of the scattering planet, this simulation was repeated with planet masses at the limits of the range, $M$ = 0.1 $M_{J}$ and $M$ = 10 $M_{J}$. The same semi-major axis ($a$ = 9.49 AU) was used for consistency.

Taking the amount of time spent at any given semi-major axis (at a resolution of 20 years) and averaging over the lifetime of the cluster (assumed to be $10^{8}$ years) and the number of bodies gives a probability for any one body to be at a given semi-major axis at any point in cluster’s life. The resulting $p(a)$ distribution is shown in Figure 7. At the range of semi-major axes that we are interested in, the probability is relatively flat, dropping less than 1 dex from 100 to 5000 AU for the $M$ = 0.59 $M_{J}$ scatterer. This is due to the fact that bodies that are scattered out to a wide semi-major axis will spend a longer time there before being scattered again, as their period becomes very large and therefore the number of interactions with the scattering planet decreases. For the larger mass scattering planet, the distribution becomes flatter but also reduces in amplitude, as bodies are ejected much more quickly from the system. In addition, the maximum semi-major axis reached before ejection moves inward, seen in the sharp decline in the $p(a)$ distribution. The lower mass scatterer keeps bodies in the system for much longer, giving a higher overall amplitude, but is slower to move the bodies to high semi-major axes, resulting in a steeper probability distribution and a decline in the distribution at larger semi-major axis.

3.2.4 Time Dependence

The passage of time in the cluster affects the circumstances related to this mechanism in two ways. First, planets are likely to begin scattering early in the lifetime of the cluster. As time increases, more planets are ejected until the system reaches a stable configuration, so the probability of having a planet at a large semi-major axis decreases. Secondly, as the cluster dissipates and the stellar density decreases, the single closest flybys become less frequent and the recent closest flybys become more frequent but with a much larger impact parameter.

By splitting the time of the cluster into two time periods, an early period from 0 to $10^{7}$ years and a late period from $10^{7}$ to $10^{8}$ years, we can capture some of this time dependence. The effects of the time period splitting on both the number and impact parameter of the flyby stars and the probability of having a planet at a given semi-major axis are shown in Figure 8 for the $M$ = 0.59 $M_{J}$ scatterer. The $p(a)$ is the same as the black curve in Figure 7, except here the spacing is 1 AU and shown linearly in semi-major axis and split between the two time periods. As expected, the probability of having a planet at a large semi-major axis decreases significantly in the late period. The $f_{strong}$ and $f_{save}$ factors are unaffected by the timing of the flyby event.

The single closest flyby, shown in the bottom panels of Figure 4, has a probability per unit time which is almost twice as much in the early period compared with the average of the late period. As seen in Figure 5, these closest flybys usually have a much lower impact parameter. This explains why the difference between $N(b_{min}<10a)$ in the early and late periods is much less at the lower end of the semi-major axis space. However, because our area of interest is at large semi-major axes, even the higher impact parameters are within a range to have a potential effect ($b/a<10$), and so the fact that the recent closest flybys are more frequent in the late period is the dominant effect. This can be seen in Figure 8 most strongly for planets with $a$ = 5000 AU, where on average there are fewer than 1 recent closest flybys within 50,000 AU ($b=10a$) in the first $10^{7}$ years, but there are $\sim$4 within that value between $10^{7}$ and $10^{8}$ years.

3.2.5 Total Efficiency

Returning to equation 1, we can now calculate an expression for the total efficiency as a function of a planet’s semi-major axis. In each time period, summing the simulated distribution of all $b_{min}$ from 0 to $10a$ gives $N(b_{min}<10a)$. The distribution obtained from the scattering simulation with $M_{scattering}$ = 0.59 $M_{J}$, split into two time periods as shown in Figure 8, gives $p(a)$. Both of these distributions are calculated with a resolution of 1 AU. Combined with $f_{strong}$ and $f_{save}$ from section 3.2.2, we find an efficiency as a function of $a$ as shown in Figure 9.

In order to consider the total expected efficiency of wide-orbit planets resulting from stellar flybys, the distribution shown in Figure 9 is integrated from 100 AU to 5000 AU ($b_{max}$/10). For the nominal cluster, this total efficiency is 0.282$\pm$0.007%. For the loose cluster, this total efficiency is 0.226$\pm$0.006%. The uncertainty cited here is solely from statistical effects and does not incorporate uncertainties related to the simplifying assumptions used in our models.

The efficiency calculation is strongly dependent on the distribution of semi-major axis used to calculate $p(a)$, as can be seen in Figure 10, which shows the efficiency calculated with the $p(a)$ term set equal to 1. The $p(a)$ distribution itself is sensitive to the properties of the scattering planet, the host star, and the population of planets formed within a system, as well as the stochastic history of planet-planet interactions. This can be seen clearly in Figure 7, where the variation of mass of the scattering planet has a large effect on the shape and magnitude of the $p(a)$ distribution.

Therefore the preceding efficiency calculation should be considered a sample calculation given one particular scattering scenario. Repeating the calculation for the other simulated $p(a)$ distributions gives a total efficiency of 0.007$\pm$0.001% ($M_{scattering}$ = 10 $M_{J}$) and 1.64$\pm$0.02% ($M_{scattering}$ = 0.1 $M_{J}$) for the nominal cluster and 0.006$\pm$0.001% ($M_{scattering}$ = 10 $M_{J}$) and 1.07$\pm$0.02% ($M_{scattering}$ = 0.1 $M_{J}$) for the loose cluster.

4 Conclusions and Discussion

4.1 Summary of Results

In this paper, we have explored how planet-planet scattering can be interrupted by passing stars, creating stable wide-orbit planets. The occurrence of planets between 100 and 5000 AU (Oort planets) by this mechanism is likely to be relatively small; we estimate it is $<$2% relative to the number of planets that are ejected. The mechanism operates in planetary systems that underwent planet-planet scattering while the host star still orbited in a stellar cluster, which may be the majority of systems that have at least one giant planet.

The more massive a planet is, the less likely it is to become an Oort planet. It requires another, more massive planet to scatter it out to large distances, and the time spent in the system prior to ejection is smaller for more massive planets. Spending less time at large distances while bound makes massive planets less likely to interact with a close stellar flyby. Loosely-bound Neptune- or super-Earth-mass planets are more likely to become stranded as Oort planets at a rate up to a few percent.

4.2 Relation to Other Work

A completely distinct way a planet might be added on a long-period orbit is by the dissolution of the cluster while free-floating planets are present. The capture probabilities and statistics of such captures have been worked out by Perets & Kouwenhoven (2012) and were seen in cluster simulations of van Elteren et al. (2019). The semi-major axis distribution of capture-formed planetary systems seen in Perets & Kouwenhoven (2012) (see their Figure 4) peaks at larger semi-major axes than we find for the Oort planets, so these captured wide-orbit planets could form part of a distinct but overlapping population as Oort planets.

Chatterjee et al. (2011) proposed a method for circularizing eccentric, long-period planets resulting from planet-planet scattering. This mechanism operates via dynamical friction with an outer disk (or with a cold planetesimal belt, i.e. Eriksson et al. 2018). The semi-major axes of such planets are limited by the size of the disk or belt, and simulations in Chatterjee et al. (2011) showed additional inward migration after circularization. Such planets would be on wide orbits, but of $\lesssim$100s of AU, mostly below the range we are considering here.

Fregeau et al. (2006) studied the cross-sections for interactions between a circular planet orbiting a host star and an intruder star of equal mass. They found a novel region of parameter space between the hard-soft boundary and the slow-fast boundary where there is about a $1/3$ probability that the planet will be stolen by the incoming star. This process occurs via an elastic (rather than head-on) collision of the two stars that — from the perspective of the planet — stops the intruder star and kicks out the original star. The other $2/3$ probability is ionization, in which the planet is stripped from its host star. This is the regime relevant to $100-1000$ AU planets ($v_{\rm orb}\simeq 1-3$ km/s) embedded in an open cluster ($v_{\sigma}\sim 1$ km/s).

When a planet is swapped into an orbit around a different star (of equal mass), Fregeau et al. (2006) calculates its semi-major axis on average increases by $\sim$25 %, and its new eccentricity is distributed thermally ($f(e)=2e$). What Fregeau et al. (2006) did not calculate is the cross-section for preserving the planet around its own host, yet changing the planet’s eccentricity by a certain amount. Since we are presuming a more highly eccentric planet ($[1-e]\ll 1$, due to planet-planet scattering) initially around the first host star, then the orbit around the new host star will tend to be more circular, and thus whatever inner planetary system the new star has might be undisturbed by it (depending, of course, on the outer edge of that system). In our simulation, the fate of the planets was only determined with respect to the initial host star; it is possible but unknown if some of the planets that became unbound from their host star due to flybys in fact became bound to the flyby star. Finally, Fregeau’s model considered equal-mass stars, which may be required for the elastic encounter to trade a planet (Levison et al., 2010). Such a coincidence in flyby stars is rare, given the broadness of the stellar initial mass function.

The similar model of Reipurth & Mikkola (2012) tracks the “dynamical unfolding” of triple-star orbits, as the gaseous core that gives them birth evaporates. There, passing stars are considered as somewhat of a nuisance, whereas we treat them a vital mechanism for raising the periastron. We found this mechanism rather ineffective for 10 $M_{J}$ planets; it is likely even less effective for stars. Recently, Rosa & Kalas (2019) proposed an interior binary HD 106906 kicked out a planet, but that planet was saved by a flying-by star (candidates were identified with Gaia). Since we found that even massive planets would not efficiently create exo-Oort planets, an interior binary star would be even less likely as the initial scatterer in our mechanism.

Gladman & Chan (2006) say rogue planets may account for high perihelion scattered-disk objects such as Sedna. Maybe scattering among planets, with a quickly-lost middle planet, can do the same in exoplanetary systems. The idea was made much more specific in a suite of numerical experiments by Silsbee & Tremaine (2018). It deposits planets at a factor of 3-10 beyond the initial distribution of planets, at least for our Solar System’s architecture.

The change of eccentricity due to an outer body has been studied in a number of contexts. Heggie & Rasio (1996) calculated the eccentricity change in a binary star that results from a flyby of a star in a cluster setting. The analytic expressions for eccentricity excitation, starting at zero eccentricity and in the limit that the flyby is slow, was derived by Kobayashi & Ida (2001). A similar calculation, but starting from an eccentricity of near unity, was done by Katz & Dong (2012). The latter’s scientific question was whether the periastron of an inner binary in a triple can be changed significantly in one orbit of that binary. Thus it is very similar to our question, though they wanted to lower the periastron, and we want to raise it. Since they assumed the outer body was bound to the inner binary, it must be moving much slower than the inner binary components, so they were able to approximate the outer body as fixed at a given location. This approximation is the opposite extreme of the usual Oort-cloud approximation (the impulse approximation) in which the intruder star flies by much faster than the orbit of the comet. We are treating flybys in which the passage can be comparable to the orbit of the planet, so the details of the calculation may be different.

Simulations of planetary systems in stellar clusters have been done. Zheng et al. (2015) modeled single planet systems in evolving stellar clusters with varying semi-major axes of the planet. They found a sharp decrease in the survival rate of planets with $a$ $\gtrsim$ 2000 AU. The planets are on circular orbits for the entire lifetime of the cluster, which would produce a more destructive effect than predicted in our mechanism. Hao et al. (2013) modeled multi-planet systems in a stellar cluster. They found that planet-planet interactions are an important factor in understanding the impact of a dense stellar environment. Their planets showed eccentricity and semi-major axis growth over 100 Myr in a cluster, consistent with our expectations. Their cluster is of constant density over 100 Myr, which increases the number of close stellar interactions and reduces the survival rate of wide-orbit planets.

4.3 Current and Potential Observations

The obvious candidate for observing this population of wide-orbit planets directly requires imaging. The giant planets involved in planet-planet scattering will emit in the infrared and are expected to be observable around young stars. Recent observational surveys for wide-orbit planets have determined an upper frequency limit of 9% for giant planets (0.5-13 $M_{J}$) at 100-1000 AU (Durkan et al., 2016) and 3% for giant planets (1-13 $M_{J}$) at 1000-5000 AU (Baron et al., 2018), rather larger than our predictions but not contradictory.

The current population of known wide-orbit planets is quite small, with only 10 confirmed planets with mass $<$13 $M_{J}$. The lowest mass of these is $\sim$6 $M_{J}$ (Ross 458 c; Burgasser et al. 2010). Planets of such high masses are very inefficiently formed via the Oort mechanism in our model, and we do not expect that any of the currently-known wide-orbit giant planets are formed via this pathway. Improved direct imaging detection limits, with sensitivities to lower masses and older star systems, could lead to detections of an Oort planet population.

Potentially, microlensing could infer this population and do so to much lower masses, with sensitivity to free-floating or loosely-bound planets down to Earth-mass objects (Mróz et al., 2018). Sumi et al. (2011) reported a large population of free-floating Jupiter-mass objects (1.8 objects per main sequence star) and, combined with direct imaging constraints on nearby potential hosts, determined that more than 75% of these objects are either unbound or bound to stars further than 500 AU. A larger survey by Mróz et al. (2017) did not find evidence of this excess of Jupiter-mass objects but did find indications of Earth-mass and super-Earth-mass free-floating planets.

The systems with a loosely bound planet by our mechanism are the same systems that bear signatures at $<$5 AU of planet-planet scattering, the primary signature being giant planets with large eccentricity (Jurić & Tremaine, 2008; Chatterjee et al., 2008; Cumming, 2010). If a system with far-out planets is found, a concerted search closer to the star is warranted, e.g., by closer-in imaging (Bryan et al., 2016) or by radial velocity.

4.4 Future Work

This work focused on reproducing the Oort mechanism for giant planets and estimating a frequency that such planets might be saved rather than ejected for sun-like systems. More physically accurate simulations across a wider range of expected system properties, done for a more statistically robust number of simulations, could provide better predictions for the efficiency of the process as well as the expected properties of Oort planets. This more in-depth evaluation of Oort planets is deferred to future work.

Oort planets can act as external perturbers to the inner systems. In cases where the semi-major axis of the stranded planet lies beyond $\sim$3000 AU, over the lifetime of the star it may be torqued back into the inner system by Galactic tides and passing stars (Heisler & Tremaine, 1986), stimulating a new epoch of scattering. Hence Oort planets may become exo-Nemesis planets. The outcome of the instability generated by this long-timescale fuse awaits future work.

We acknowledge support of grant NASA-NNX17AB93G through NASA’s Exoplanet Research Program. Simulations in this paper made use of the REBOUND code which can be downloaded freely at

http://github.com/hannorein/rebound.

Appendix A Correcting for Variable Density

In determining the wait time between flyby stars, the fact that density is decreasing over time must be taken into account. An initial uncorrected wait time, $\tau_{o}$, is calculated using equation 2, assuming that the density is constant at some initial density, $n_{o}$. An example considering a flyby star during the last 1% of a cluster’s life is shown as a blue rectangle in Figure 11, where $n_{init}=25\>\text{stars}\,\text{pc}^{-3}$. The additional delay due to the decaying density, assuming $\dot{n}$ is a negative constant over the life of the cluster, is found by calculating corrected wait time, $\tau$, that has the equivalent area underneath the density curve, as shown by the magenta trapezoid in Figure 11.

[TABLE]

Solving for $\tau$ leads to

[TABLE]

where the positive root is the corrected wait time.

This process is iterated over the lifetime of the cluster, with $n_{o}$ for each new flyby taken to be the density at the time of the most recent flyby and a new random $\tau_{o}$ generated until $\tau$ implies a wait time beyond the age of the cluster.

Appendix B Probability Density Function of $b_{min}$

Consider $x_{i}$ drawn from a uniform distribution between 0 and 1. The probability that $x_{i}$ is greater than $x$ is given by $1-x$. After $N$ such draws, reordering the results by index number, the probability that $x_{o}$ (the smallest number drawn) is greater than $x$ is given by

[TABLE]

Therefore the probability that $x_{o}$ is smaller than $x$ is given by

[TABLE]

The probability density of $x$ being the smallest number after N draws is then computed as the derivative of $p(x_{o}<x)$.

[TABLE]

This can be expanded to consider $f(x)=2x$ for $x$ drawn from a uniform distribution between 0 and 1. This is the same as drawing $x^{2}$ uniformly between 0 and 1, which is what we do for $x=\frac{b}{b_{max}}$. Now the probability that $x_{i}$ is greater than $x$ is given by integrating $f(x)$ from $x$ to 1.

[TABLE]

The same logic as before can be applied to find the probability of $x_{o}$ being less than $x$.

[TABLE]

Finally, taking the absolute value of the derivative with respect to x, the probability density of $x$ being the smallest number after N draws is given by

[TABLE]

Expanding further to a number $x^{2}$ drawn uniformly between 0 and $b_{max}^{2}$, we introduce a normalization factor $A$ and integrate over the interval to normalize:

[TABLE]

We find $A=b_{max}^{-2N}$ and therefore the probability density function can be written as

[TABLE]

This is the same as equation 3 where $x=b_{min}$. Now, as we consider the limit where $N\rightarrow\infty$ and $b_{max}\rightarrow\infty$, we can establish a constant $d=\frac{b_{max}^{2}}{N}$. Equation B8 can then be rewritten as

[TABLE]

Given the limit definition of Euler’s formula,

[TABLE]

we can rewrite B9 as

[TABLE]

The probability density function of a Rayleigh distribution is

[TABLE]

so the probability density function of $b_{min}$ is a Rayleigh distribution with $\sigma=\sqrt{\frac{b_{max}^{2}}{2N}}$.

This derivation was inspired by the technique of finding the distribution of the highest peak among many in a periodogram (Scargle, 1982, Section IIIa). A branch of statistics known as record statistics (e.g., Chandler 1952) can be put to service to find the distribution of the $N$th “recent closest passage,” though in practice we resorted to numerical draws to determine the distribution used for the efficiency calculation and illustrated in Figure 5.