In-Situ Scattering of Warm Jupiters and Implications for Dynamical Histories

TL;DR

This study investigates how planet-planet scattering at sub-AU distances influences the eccentricities of warm Jupiters, revealing that a significant fraction may have undergone in-situ scattering, affecting their dynamical histories.

Contribution

It provides a systematic analysis of scattering outcomes for warm Jupiters at sub-AU distances, including effects of relativistic precession, and links these to observed eccentricity distributions.

Findings

Approximately 60% of solitary WJs may have experienced in-situ scattering.

Two-planet systems are formed mainly through collisions with low eccentricities.

One-planet systems result from ejections and collisions, with higher eccentricities.

Abstract

Many warm Jupiters (WJs) have substantial eccentricities, which are linked to their formation and migration histories. This paper explores eccentricity excitation of WJs due to planet-planet scattering, beginning with 3-4 planets in unstable orbits, with the innermost planet placed in the range (0.1 - 1)AU. Such a setup is consistent with either in-situ formation or arrival at sub-AU orbits due to disk migration. Most previous N-body experiments have focused on "cold" Jupiters at several AU, where scattering results in planet ejections, efficiently exciting the eccentricities of surviving planets. In contrast, scattering at sub-AU distances results in a mixture of collisions and ejections, and the final eccentricities of surviving planets are unclear. We conduct scattering experiments for a range of planet masses and initial spacings, including the effect of general relativistic apsidal precession, and systematically catalogue the scattering outcomes and properties of surviving planets. A comparable number of one-planet and two-planet systems are produced. Two-planet systems arise exclusively through planet-planet collisions, and tend to have low eccentricities/mutual inclinations and compact configurations. One-planet systems arise through a combination of ejections and collisions, resulting in higher eccentricities. The observed eccentricity distribution of solitary WJs (lacking detection of a giant planet companion) is consistent with roughly 60% of the systems having undergone in-situ scattering, and the remaining experiencing a quiescent history.

Tables

Table 1. Table 1: Parameters and results for the different sets of simulations. In all simulations, the initial innermost planet semi-major axis is sampled uniformly in a 1 = [ 0.1 − 1 ] subscript 𝑎 1 delimited-[] 0.1 1 a_{1}=[0.1-1] AU, and the initial inclinations are sampled uniformly in the range [ 0.1 ∘ − 2 ∘ ] delimited-[] superscript 0.1 superscript 2 [0.1^{\circ}-2^{\circ}] . For fiducial , near-eq-mass , lognorm-mass , and fiducial-K-3 , the initial eccentricities are sampled in the range [ 0.01 − 0.05 ] delimited-[] 0.01 0.05 [0.01-0.05] . For fiducial-K-5 , we increase the initial eccentricities so that instabilities occur within a practical amount of time (see Section 2.4 ). The lognorm-mass  simulations have all masses sampled from a log-normal distribution, with a mean of log ⁡ ( m p / M J ) = 0 subscript 𝑚 p subscript 𝑀 J 0 \log(m_{\rm p}/M_{\rm J})=0 , width 0.5 0.5 0.5 , and upper and lower limits of m p = 4 ​ M J subscript 𝑚 p 4 subscript 𝑀 J m_{\rm p}=4M_{\rm J} and m p = 0.25 ​ M J subscript 𝑚 p 0.25 subscript 𝑀 J m_{\rm p}=0.25M_{\rm J} . The columns in the table are (from left to right), the simulation name, number of simulations ( N run subscript 𝑁 run N_{\rm run} ), choice of planet masses, initial orbital spacing in units of mutual Hill radii ( K 𝐾 K , see equations [ 1 ]-[ 2 ]), fraction of one and two-planet systems produced at the end of “Phase 2” of the integration ( f 1 ​ p subscript 𝑓 1 p f_{\rm 1p} and f 2 ​ p subscript 𝑓 2 p f_{\rm 2p} ). The columns labeled f pp subscript 𝑓 pp f_{\rm pp} , f ps subscript 𝑓 ps f_{\rm ps} , and f ej subscript 𝑓 ej f_{\rm ej} indicate the fraction of systems in which the first (or only) planet loss resulted from a planet-planet collision, planet-star collision, or planet ejection respectively (see Fig.  2 ).

Name	$N_{\mathrm{run}}$	Initial Masses $\left(M_J\right)$	$K$	$f_{1p}$	$f_{2p}$	$f_{\mathrm{pp}}$	$f_{\mathrm{ps}}$	$f_{\mathrm{ej}}$
fiducial	3313	$0.5,1.0,2.0$	4	0.52	0.48	0.65	0.06	0.29
near-eq-mass	973	$0.9,1.0,1.1$	4	0.34	0.66	0.78	0.10	0.12
lognorm-mass	965	See caption	4	0.50	0.49	0.69	0.00	0.31
fiducial-K-3	986	$0.5,1.0,2.0$	3	0.60	0.41	0.76	0.00	0.24
fiducial-K-5	949	$0.5,1.0,2.0$	5	0.50	0.48	0.63	0.00	0.37

Table 2. Table 2: Scattering outcomes and properties of the two-planet systems at the end of “Phase 2” of the integration. The quantities f pp , 2 ​ p subscript 𝑓 pp 2 p f_{\rm pp,2p} , f ps , 2 ​ p subscript 𝑓 ps 2 p f_{\rm ps,2p} , f ej , 2 ​ p subscript 𝑓 ej 2 p f_{\rm ej,2p} are the branching ratios for planet-planet collisions, planet-star collisions, and planet ejections respectively (see Fig.  2 ). The remaining columns show the inner and outer planet eccentricities ( e in , e out subscript 𝑒 in subscript 𝑒 out e_{\rm in},e_{\rm out} ), the mutual inclination I mut subscript 𝐼 mut I_{\rm mut} , and the semi-major axis ratio α ≡ a out / a in 𝛼 subscript 𝑎 out subscript 𝑎 in \alpha\equiv a_{\rm out}/a_{\rm in} . The superscript “ave” denotes the unweighted average, and the superscripts 50 , 90 50 90 50,90 indicate the 50th (median), and 90th percentiles.

Name	$f_{\mathrm{pp},2p}$	$f_{\mathrm{ps},2p}$	$f_{\mathrm{ej},2p}$	$e_{\mathrm{in}}^{\mathrm{ave}}$	$e_{\mathrm{in}}^{50}$	$e_{\mathrm{in}}^{90}$	$e_{\mathrm{out}}^{\mathrm{ave}}$	$e_{\mathrm{out}}^{50}$	$e_{\mathrm{out}}^{90}$	$I_{\mathrm{mut}}^{\mathrm{ave}}$ $(\circ )$	$I_{\mathrm{mut}}^{50}$ $(\circ )$	$I_{\mathrm{mut}}^{90}$ $(\circ )$	${\alpha }_{\mathrm{ave}}$	${\alpha }_{50}$	${\alpha }_{90}$
fiducial	0.48	0.0	0.0	0.15	0.11	0.34	0.12	0.09	0.27	5.25	2.46	13.26	2.46	2.06	3.02
near-eq-mass	0.66	0.0	0.0	0.18	0.14	0.38	0.15	0.12	0.3	7.16	3.01	18.39	2.62	2.04	3.69
lognorm-mass	0.49	0.0	0.0	0.14	0.1	0.31	0.12	0.08	0.25	5.20	2.28	14.73	2.3	2.01	2.74
fiducial-K-3	0.40	0.0	0.0	0.13	0.08	0.32	0.09	0.06	0.22	4.99	2.34	12.56	1.93	1.64	2.5
fiducial-K-5	0.48	0.0	0.0	0.15	0.12	0.31	0.15	0.11	0.33	5.14	2.16	12.02	2.63	2.37	3.42

Table 3. Table 3: Scattering outcomes and properties of the one-planet systems at the end of “Phase 2” of the integration. The quantities f pp , 1 ​ p ( 1 ) superscript subscript 𝑓 pp 1 p 1 f_{\rm pp,1p}^{(1)} , f ps , 1 ​ p ( 1 ) superscript subscript 𝑓 ps 1 p 1 f_{\rm ps,1p}^{(1)} and f ej , 1 ​ p ( 1 ) superscript subscript 𝑓 ej 1 p 1 f_{\rm ej,1p}^{(1)} denote the branching ratios for the first planet loss (see Fig.  2 ). The quantities f pp , 1 ​ p ( 2 ) superscript subscript 𝑓 pp 1 p 2 f_{\rm pp,1p}^{(2)} , f ps , 1 ​ p ( 2 ) superscript subscript 𝑓 ps 1 p 2 f_{\rm ps,1p}^{(2)} and f ej , 1 ​ p ( 2 ) superscript subscript 𝑓 ej 1 p 2 f_{\rm ej,1p}^{(2)} have identical meanings, but for the second planet loss. The remaining columns indicate the average, median, 10th, and 90th percentile eccentricities.

Name	$f_{\mathrm{pp},1p}^{\left(1\right)}$	$f_{\mathrm{ps},1p}^{\left(1\right)}$	$f_{\mathrm{ej},1p}^{\left(1\right)}$	$f_{\mathrm{pp},1p}^{\left(2\right)}$	$f_{\mathrm{ps},1p}^{\left(2\right)}$	$f_{\mathrm{ej},1p}^{\left(2\right)}$	$e_{\mathrm{ave}}$	$e_{50}$	$e_{10}$	$e_{90}$
fiducial	0.17	0.06	0.29	0.01	0.03	0.48	0.30	0.27	0.10	0.55
near-eq-mass	0.12	0.10	0.11	0.01	0.02	0.31	0.47	0.42	0.23	0.78
lognorm-mass	0.2	0.0	0.3	0.04	0.0	0.46	0.30	0.26	0.10	0.56
fiducial-K-3	0.35	0.0	0.24	0.06	0.0	0.53	0.29	0.26	0.10	0.55
fiducial-K-5	0.15	0.0	0.36	0.03	0.0	0.47	0.30	0.26	0.09	0.55

Table 4. Table 4: Properties of both one and two-planet systems produced from the 4-planets  simulations (see Section 2.7 ). The columns labeled f 1 ​ p subscript 𝑓 1 p f_{\rm 1p} and f 2 ​ p subscript 𝑓 2 p f_{\rm 2p} have the same meanings as in Table LABEL:table1 . The next three columns indicate the average, 50th and 90th percentile eccentricities of the one-planet systems. The remaining columns indicate properties of two-planet systems, with the same notation as in Table LABEL:table_2p .

Name	$f_{1p}$	$f_{2p}$	$e_{\mathrm{ave}}$	$e_{50}$	$e_{90}$	$e_{\mathrm{in}}^{\mathrm{ave}}$	$e_{\mathrm{in}}^{50}$	$e_{\mathrm{in}}^{90}$	$e_{\mathrm{out}}^{\mathrm{ave}}$	$e_{\mathrm{out}}^{50}$	$e_{\mathrm{out}}^{90}$	$I_{\mathrm{mut}}^{\mathrm{ave}}$	$I_{\mathrm{mut}}^{50}$	$I_{\mathrm{mut}}^{90}$	${\alpha }_{\mathrm{out}}^{\mathrm{ave}}$	${\alpha }_{\mathrm{out}}^{50}$	${\alpha }_{\mathrm{out}}^{90}$
4-planets	0.75	0.23	0.35	0.31	0.64	0.2	0.17	0.4	0.2	0.17	0.4	7.21	4.01	15.76	4.13	2.87	4.98

Equations

$$a_i-a_{i-1}=K{R}_{\mathrm{H},\mathrm{mut}},$$

$$R_{\rm H,mut}=\frac{a_{i-1}+a_{i}}{2}\bigg{(}\frac{m_{i-1}+m_{i}}{3M_{\star}}\bigg{)}^{1/3}.$$

$$K_{\mathrm{RV}}=\frac{28.4m/s}{\sqrt{1-e^2}}\left(\frac{m_{\mathrm{p}}\sin i_{\mathrm{sky}}}{M_{\mathrm{J}}}\right){\left(\frac{P}{1\mathrm{yr}}\right)}^{-1/3}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{-2/3}.$$

$$f(e)=\mathcal{F}{f}_{\mathrm{circ}}(e)+(1-\mathcal{F})f_{\mathrm{scat}}(e),$$

$$\mathcal{L}=\displaystyle\prod_{i=1}^{n_{\rm obs}}\bigg{[}\mathcal{F}f_{\rm circ}(e_{i})+(1-\mathcal{F})f_{\rm scat}(e_{i})\bigg{]},$$

$$v_{\infty }^2=v_{\mathrm{coll}}^2-v_0^2,$$

$$b_{\infty }{v}_{\infty }=b_{\mathrm{coll}}{v}_{\mathrm{coll}},$$

Full text

††pagerange: In-Situ Scattering of Warm Jupiters and Implications for Dynamical Histories–LABEL:lastpage††pubyear: 2019

1 Introduction

The origin of warm Jupiters (WJs, giant planets with orbital periods in the range of $\sim 10-300$ days) remains an unsolved problem in exoplanetary science. Whether WJs previously migrated from farther out, are currently in the process of migration, or formed in-situ, is uncertain. In-situ formation of hot Jupiters (HJs) and other close-in planets was considered previously by Lee et al. (2014); Batygin et al. (2016); Boley et al. (2016); Lee & Chiang (2016). Giant planet migration has been extensively studied in the context of HJ formation, and comes in two types. One is disk migration, in which planets are transported inwards due to torques from the protoplanetary disk (e.g. Lin et al., 1996; Tanaka et al., 2002; Kley & Nelson, 2012; Baruteau et al., 2014). The second is high-eccentricity migration, in which the planet’s eccentricity is excited to an extreme value by a stellar or planetary companion(s), so that tides raised on the planet at pericenter passages shrink and circularize the orbit. High-eccentricity migration comes in several distinct flavors, depending on the details of the eccentricity excitation. Possibilities include excitation from an inclined companion due to Lidov-Kozai cycles (Lidov, 1962; Kozai, 1962) or other secular perturbations (Wu & Murray, 2003; Fabrycky & Tremaine, 2007; Correia et al., 2011; Naoz et al., 2012; Petrovich, 2015a, b; Anderson et al., 2016; Muñoz et al., 2016; Hamers et al., 2017; Vick et al., 2019), planet-planet scatterings, possibly combined with secular interactions (Rasio & Ford, 1996; Nagasawa et al., 2008; Nagasawa & Ida, 2011; Beaugé & Nesvorný, 2012), and secular chaos (Lithwick & Wu, 2011, 2014; Teyssandier et al., 2019). See also Dawson & Johnson (2018) for a review.

WJs are observed to have a wide range of eccentricities. The radial velocity WJs from the Exoplanet Orbit Database111exoplanets.org, accessed July 4, 2019. (with minimum planet masses above $0.3M_{\rm J}$ and semi-major axes in the range $0.1-1$ AU) have an average eccentricity $e\simeq 0.24$, and $30\%$ of the planets have $e>0.3$, with a maximum $e\simeq 0.93$. Theories of planet formation and migration must be able to account for these observations.

Several mechanisms have been proposed in exciting WJ eccentricities. One seemingly natural explanation is high-eccentricity migration. In this scenario, WJs are caught in the act of inward migration, eventually to become HJs in circular orbits. In any high-eccentricity migration theory, a small pericenter distance is required, so that tidal dissipation may shrink the orbit within the lifetime of the host star. The majority of observed WJs have pericenter distances too large to allow for efficient tidal dissipation. These WJs may be undergoing secular eccentricity oscillations driven by an external companion, so that they are observed during a low-eccentricity phase of the secular cycle (Dong et al., 2014; Dawson & Chiang, 2014; Petrovich & Tremaine, 2016). This scenario requires a relatively close or massive companion, so that the secular perturbations that lead to eccentricity oscillations are not suppressed by general relativistic apsidal precession.

As a planet is transported inward via high-eccentricity migration, it eventually becomes decoupled from the companion. Unless the perturber is particularly “strong” (i.e. close or massive), this decoupling has usually occurred by the time the migrating planet reaches WJ territory, causing the eccentricity to freeze at the maximum value. With a small (fixed) pericenter distance, the planet migrates rapidly into HJ territory, with little time spent in the WJ phase (see, e.g., Fig. 1 of Anderson et al., 2016)222When dynamical (chaotic) tides are incorporated into the high-eccentricity migration model, highly eccentric ($e\gtrsim 0.9$) WJs form rapidly, and subsequently circularize to become HJs (Vick et al., 2019).. Consequently, forming WJs that are likely to be observed requires a somewhat specific initial architecture. The perturber must be strong enough to drive eccentricity oscillations at sub-AU distances, but weak enough to ensure dynamical stability of the initial configuration.

Indeed, HJ population syntheses that have chosen a broad range of properties for the perturber do not produce any appreciable numbers of WJs (Petrovich, 2015a; Anderson et al., 2016; Hamers et al., 2017). Observationally, the number of WJs exceeds the number of HJs. The study by Petrovich & Tremaine (2016) attempts to replicate the observed relative numbers of HJs and WJs by choosing as their initial condition a population of planets initially located at $1$ AU with strong perturbers in a narrow range of semi-major axis, located at $5-6$ AU. However, these initial conditions still produce too few WJs compared to HJs.

Another difficulty of high-eccentricity migration in the context of WJ formation is the fact that unlike HJs, a large fraction of WJs have one or more close, low-mass neighbors (Huang et al., 2016). Such a configuration is difficult to produce in a violent high-eccentricity migration scenario. Furthermore, Antonini et al. (2016) examined the subset of observed WJs with characterized (in terms of mass, semi-major axis and eccentricity) external giant planet companions, and found that most systems are inconsistent with a traditional high-eccentricity migration origin (where the WJ formed beyond 1 AU), due to the initial configurations being dynamically unstable.

These difficulties of high-eccentricity migration in reproducing properties of WJs indicate that in-situ formation or disk migration may be responsible for forming many, if not most WJs. However, both in-situ formation and disk migration by themselves have difficulties in producing eccentric WJs. As a result, mechanisms for exciting eccentricity are needed. Planet-disk interactions are capable of exciting eccentricities, but are limited to modest values $\lesssim 0.2$ (Goldreich & Sari, 2003; Tsang et al., 2014; Duffell & Chiang, 2015; Ragusa et al., 2018). Recently, Petrovich et al. (2019) studied a mechanism for transferring eccentricity from an outer planet to a WJ during the dispersal of a massive protoplanetary disk.

In Anderson & Lai (2017), we considered the possibility that eccentric WJs arise due to secular perturbations from a distant companion, without requiring that the WJs be undergoing high-eccentricity migration. Eccentricity may be excited by a highly inclined perturber via Lidov-Kozai cycles, or by an eccentric coplanar perturber due to an apsidal precession resonance. Taking the sample of WJs with external planetary companions with characterized orbits, and assuming that the WJ formed in a circular orbit, we found that relatively high mutual inclinations are needed ($\sim 50^{\circ}-60^{\circ}$) to generate the observed eccentricities. This finding is intriguing, and consistent with previous evidence for high inclinations in many of the same systems found by Dawson & Chiang (2014). However, generating such inclinations is non-trivial, and requires an early scattering event in the system’s history.

Planet-planet scattering itself remains another possibility in producing eccentric WJs (Mustill et al., 2017; Frelikh et al., 2019), or scattering in combination with dynamical tides (Marzari & Nagasawa, 2019). A substantial literature of giant planet scattering work exists, ranging from scattering of two planets to ten or more planets (e.g. Chambers et al., 1996; Lin & Ida, 1997; Ford et al., 2001; Adams & Laughlin, 2003; Chatterjee et al., 2008; Ford & Rasio, 2008; Jurić & Tremaine, 2008; Nagasawa & Ida, 2011), but most studies have focused on “cold Jupiters,” giant planets located at several AU. Giant planet scattering close to the host star is less explored.

Scattering outcomes depend on the “Safronov number”, the squared ratio of the escape velocity from the planetary surface to the planet’s orbital velocity. When the Safronov number is much less than unity, close encounters between planets result in planet-planet collisions, with the collision product having a low eccentricity. When the Safronov number is much greater than unity, planet ejections are expected, efficiently raising the eccentricities of the remaining planets. WJs have Safronov numbers of order unity and thus lie in a regime in which a combination of collisions and ejections may occur, so that the degree of eccentricity excitation is unclear. Petrovich et al. (2014) undertook a scattering study of primarily HJs, with close initial spacings, and found inefficient eccentricity excitation due to a preponderance of collisions. Whether this finding holds for WJs with a wider range of initial spacing has yet to be thoroughly explored.

This paper presents a systematic study of planet-planet scattering for systems of closely-spaced WJs starting with initially low eccentricities. This setup is consistent with either in-situ formation or arrival at a sub-AU orbit by disk migration. The goal of this paper is two-fold: (1) On observationally-motivated grounds, we aim to identify to what extent planet scattering may be contributing to eccentric WJs; (2) On the theoretical side, we aim to catalogue the scattering outcomes for planets in the range $0.1-1$ AU, where a rich variety of collisions and ejections are expected. We conduct N-body scattering experiments of three or four giant planets, including the effect of general relativistic apsidal precession. We explore a variety of choices for planet masses and initial spacings, determine the branching ratios for various outcomes, and analyze the properties of the “remnant” systems after scattering. Radial velocity observations have yielded samples of solitary WJs and WJs with an external giant planet companion, with eccentricity measurements for both samples. We compare the results of the scattering experiments with the observed system properties.

This paper is organized as follows. In Section 2, we describe the setup of our N-body calculations, present branching ratios, demonstrate how the scattering results depend on various parameters, and discuss the properties of the surviving planetary systems. In Section 3 we compare the results of Section 2 with observations. We present our conclusions in Section 4.

2 Scattering Experiments

2.1 Setup & Canonical Parameters

We begin with a system of three planets, with masses $m_{1}$, $m_{2}$, $m_{3}$, orbiting a host star with mass $M_{\star}=1M_{\odot}$ and radius $R_{\star}=1R_{\odot}$. For each planet, we sample the initial eccentricities uniformly in the range $[0.01,0.05]$, inclinations (relative to an arbitrary reference frame) uniformly in the range $[0^{\circ},2^{\circ}]$, and argument of pericenter, longitude of ascending node, and mean anomaly uniformly in $[0,2\pi]$. The initial semi-major axes are specified in units of mutual Hill radius $R_{\rm H,mut}$, so that

[TABLE]

where

[TABLE]

We adopt a fiducial spacing of $K=4$. The innermost semi-major axis $a_{1}$ is sampled uniformly in the range $[0.1-1]$AU. The planet masses are chosen to be $0.5$, $1$ and $2$ $M_{\rm J}$, with randomly assigned ordering. We draw a sample of over 3000 systems with these parameters and evolve using N-body integrations. This set of simulations constitutes our fiducial sample, which we will refer to as fiducial (see also Table LABEL:table1). Sections 2.4, 2.5 and 2.6 explore how the results depend on $K$, $a_{1}$, and planet masses, and Section 2.7 considers scattering of four planets.

The choice of these initial conditions is broadly consistent with an initial ensemble of sub-AU giant planets that is stabilized under the presence of a gaseous disk, and later undergoes instabilities following disk dispersal. The ranges for the initial eccentricities and inclinations are chosen to be small, but non-zero, reflecting that small but finite eccentricities and inclinations are expected under some circumstances (e.g. by planet-disk interactions). Furthermore, these ranges, as well as our choice to space the orbits in fixed units of mutual Hill radius, are similar to previous “cold Jupiter” N-body scattering experiments, allowing for comparison between this paper and previous work.

The N-body calculations are performed using rebound (Rein & Liu, 2012). We include the effects of apsidal precession due to general relativity using the gr-potential option in reboundx333https://github.com/dtamayo/reboundx. When the separation between any two bodies becomes less than the sum of their radii, we assume the bodies merge, conserving mass and momentum, as in the built-in rebound collision routine (see also Appendix A). Since we consider young giant planets, we set the radius of each planet to $R_{p}=1.6R_{\rm J}$. Planets are considered ejected if the distance from the host star exceeds $1000$ AU, and are subsequently removed from the simulation.

We integrate systems of three unstable planets using the IAS15 integrator in rebound (Rein & Spiegel, 2015) for a timespan of $10^{6}P_{1}$, with $P_{1}$ the initial orbital period of the innermost planet. We refer to this initial, highly unstable phase (with close encounters eventually resulting in collisions or ejections) as “Phase 1” of the integration. After this phase nearly $100\%$ of the three-planet systems have become destabilized due to collisions or ejections. We continue to integrate the remaining two-planet systems for a timespan of $10^{8}P_{\rm in}$ (or until another collision or ejection has occurred), where $P_{\rm in}$ is the orbital period of the inner planet at the end of “Phase 1”. We refer to this longer-term integration as “Phase 2,” and use the hybrid integrator Mercurius. Mercurius utilizes a symplectic Wisdom-Holman integrator WHFAST for large separations between planets (Rein & Tamayo, 2015; Wisdom & Holman, 1991), switching to IAS15 when the separation between any two bodies becomes less than a critical value. We choose this critical value to be $5$ Hill radii. The timestep for the Wisdom-Holman integrator is chosen to be $0.02P_{\rm in}$. Repeating a subset of the fiducial sample with a timestep $0.01P_{\rm in}$ yielded statistically identical results.

2.2 Scattering Outcomes

Figure 1 shows the fractions of one, two, and three-planet systems as a function of time. The left panel shows “Phase 1” of the integration using IAS15. Since the initial systems are highly unstable, the fraction of three-planet systems quickly decays, eventually reaching a negligible value after $10^{6}$ initial orbital periods of the innermost planet. During the long-term follow-up integration of the two-planet systems (“Phase 2”, right panel of Fig. 1), the fraction of two-planet systems declines, eventually approaching a constant value after $10^{8}$ orbits. The vast majority of the two-planet systems are undergoing secular interactions at this point, with constant semi-major axes and oscillating eccentricities and inclinations. We thus conclude that the majority of the remaining 2-planet systems are stable at this time.

We note that a number of these two-planet systems could undergo instabilities over much longer timescales. Several empirical stability criteria for two-planets in eccentric/inclined orbits have been proposed (e.g. Mardling & Aarseth, 2001; Petrovich, 2015c), but they are only reliable sufficiently far from the “fuzzy” stability boundary. Using the Petrovich (2015c) stability criterion, we find that over $90\%$ of our fiducial two-planet systems are classified as unstable, but fall in the uncertain regime of parameter space near the boundary. Thus, the Petrovich (2015c) stability criterion is overly conservative in this application. In addition, his stability criterion does not consider the potentially stabilizing effects of GR apsidal precession, as well as the possibility of planet-planet collisions. As a result, full numerical integrations over long timescales are needed. To evaluate how likely instabilities may occur over much longer timescales, we randomly choose 30 of the closely-spaced two-planet systems (with $a_{\rm out}/a_{\rm in}<3$) from the fiducial sample, and integrate for an additional $10^{9}$ orbits of the inner planet, i.e. an order of magnitude longer. Of these 30 systems, only 3 became destabilized. We note that $10^{9}$ orbits at typical WJ distances may still be relatively short compared to the lifetimes of observed systems. We thus conclude that at least $\sim 10\%$ of the two-planet systems are expected to eventually go unstable.

The fractions of one and two-planet systems at the end of Phase 2 are $f_{\rm 1p}=52\%$ and $f_{\rm 2p}=48\%$ respectively. Figure 2 shows the chain of events and branching ratios that produced these one and two-planet systems. A close encounter can result in a planet “loss” due to a planet-planet collision, planet-star collision, or planet ejection. Considering the first planet loss only (so that the total number of planets is reduced from three to two), we introduce the variables $f_{\rm pp}$, $f_{\rm ps}$, $f_{\rm ej}$, which quantify the fractions of three-planet systems resulting in planet-planet collisions, planet-star collisions, and planet ejections respectively. Fig. 2 shows that the most common scattering outcome of the original three-planet systems is a planet-planet collision ($f_{\rm pp}\simeq 65\%$), followed by planet ejection ($f_{\rm ej}\simeq 29\%$) and planet-star collision ($f_{\rm ps}\simeq 6\%$). See also Table LABEL:table1. Note that following the first planet loss, some systems may later lose a second planet, which is not included in the calculation of $f_{\rm pp}$, $f_{\rm ps}$, $f_{\rm ej}$.

In our numerical experiments, the two-planet systems arise entirely due to planet-planet collisions. The one-planet systems arise from a combination of ejections and collisions, with a modest preference for ejections for the first planet loss, and an overwhelming preference for ejections for the second planet loss. Branching ratios separated into results for one and two-planet systems are also listed in Tables LABEL:table_2p and LABEL:table_1p.

2.3 Properties of One and Two-Planet Systems

In this section we analyze the orbital characteristics of the one and two-planet systems at the end of Phase 2 of the integration. In most figures we also show properties of observed WJ systems; however, we delay comparison with the observed systems until Section 3.

Figure 3 shows the properties of the two-planet systems. The subscript “in” indicates the inner planet and “out” indicates the outer planet. We comment on the key features here, but see the caption and Table LABEL:table_2p for additional details. In Fig. 3a, we present the relative spacing, $a_{\rm out}$ versus $a_{\rm in}$. The two-planet systems are quite closely-spaced, with $90\%$ of systems satisfying $a_{\rm out}/a_{\rm in}<3$. In Fig. 3c, we show semi-major axis versus eccentricity for the inner planet ($a_{\rm in}$ versus $e_{\rm in}$). Smaller values of $a_{\rm in}$ lead to a wider range of eccentricities. The eccentricity distribution in Fig. 3e is peaked at small values, and has mean and 90th percentile values $e_{\rm in}^{\rm ave}=0.15$, $e_{\rm in}^{90}=0.34$. The mutual inclinations of the two-planet systems (Fig. 3g) have mean and 90th percentiles $I_{\rm mut}^{\rm ave}=5.25^{\circ}$ and $I_{\rm mut}^{90}=13.26^{\circ}$. The typically low eccentricities and inclinations are consistent with the fact that the two-planet systems arose entirely due to planet-planet collisions (Ford et al., 2001).

Figure 4 shows the properties of the one-planet systems. Most of these planets are WJs, but a small fraction have semi-major axes greater than 1 AU. Compared to the two-planet systems, the degree of eccentricity excitation for the one-planet systems is higher, due to the fact that each one-planet system has previously suffered at least one ejection. The distribution of eccentricities shown in Fig. 4 peaks near $0.25$ and has a long tail. The mean and 90th percentile eccentricities of the single-planet systems are $e_{\rm ave}=0.3$ and $e_{90}=0.55$ respectively.

2.4 Dependence on initial spacing $K$

We now examine how the results depend on the initial spacing parameter $K$. Previous work has shown that when planet collisions are infrequent, $K$ primarily affects the time for instabilities (Chambers et al., 1996), rather than the degree of eccentricity excitation itself. However, the situation is unclear when planet collisions are common.

The fiducial simulations have $K=4$. To explore how the results depend on $K$, we perform additional sets of $\sim 1000$ simulations with $K=3$ and $K=5$ (titled fiducial-K-3 and fiducial-K-5 respectively), keeping the other parameters the same as fiducial. However, for fiducial-K-5, the timescale for instabilities to develop is sufficiently long to be impractical, so we slightly increase the initial eccentricities of the planets: We set the initial $e_{1}=0.08$ and assign the $e_{2}$ and $e_{3}$ so that $a_{3}(1-e_{3})-a_{2}(1+e_{2})=3R_{\rm H,mut}$. Depending on the mass ordering, the initial eccentricities of the outer planets are at most $\sim 0.11$ Such eccentricities may conceivably develop due to planet-disk interactions.

Inspection of Tables LABEL:table1 - LABEL:table_1p reveals that the value of $K$ does not drastically affect the results. As expected, the fraction of planet-planet collisions is highest for the smallest $K$, with $f_{\rm pp}=0.76$ for $K=3$, compared to $f_{\rm pp}=0.63$ for $K=5$. The eccentricities listed in Tables LABEL:table_2p and LABEL:table_1p do not differ appreciably. The final spacing of the two-planet systems ($\alpha\equiv a_{\rm out}/a_{\rm in}$) is mildly dependent on $K$, with average $\alpha_{\rm ave}=1.93,2.46,2.63$ for $K=3,4,5$, and 90th percentile spacing $\alpha_{90}=2.5,3.02,3.42$.

2.5 Dependence on initial innermost semi-major axis

Next we discuss how the scattering results depend on the initial innermost planet semi-major axis $a_{1}$. Taking fiducial (with $K=4$) we split the simulations ($3313$ total) into four bins of $a_{1}$, each with width $0.225$ AU and centered at $0.21$, $0.44$, $0.66$, and $0.89$ AU.

Figure 5 shows how the results depend on $a_{1}$. The top panel shows the frequencies of collisions and ejections. In the first bin, centered at $0.21$ AU, planet-planet collisions are very common, with $f_{\rm pp}\simeq 82\%$ and $f_{\rm ej}\simeq 15\%$. As $a_{1}$ increases, the frequency of ejections increases, so that in the last bin, centered at $0.89$ AU, $f_{\rm pp}\simeq 50\%$. The fraction of planet-star collisions remains nearly constant at $f_{\rm ps}\simeq 6\%$ across all bins.

The dependency of $f_{\rm pp}$ and $f_{\rm ej}$ on $a_{1}$ leads to a variation of the relative numbers of one and two-planet systems, as shown in the middle panel of Fig. 5. In the first $a_{1}$ bin, we have $f_{\rm 1p}=37\%$ and $f_{\rm 2p}=63\%$. The larger value of $f_{\rm 2p}$ at small $a_{1}$ is a direct consequence of the fact that collisions are most common close to the host star, leading to systems with lower eccentricities, and hence, more stable two-planet systems. As $a_{1}$ increases, the number of ejections increases, so that $f_{\rm 1p}$ increases, reaching $f_{\rm 1p}=63\%$ and $f_{\rm 2p}=37\%$ in the last bin. This dependency of collisions/ejections with $a_{1}$ leads to a slight increase in the average eccentricity with $a_{1}$, as shown in the bottom panel of Fig. 5. Note that the eccentricities shown in Fig. 5 include both the one and two-planet systems.

Figure 6 illustrates how the initial $a_{1}$ affects the properties of the final two-planet systems. The initial $a_{1}$ primarily determines the final semi-major axis, since scattering typically limits changes in semi-major axis to factors of order unity. The final planet eccentricities, mutual inclination, and relative spacing depend weakly on $a_{1}$. A similar weak dependence on $a_{1}$ for one-planet systems is also found (not shown).

2.6 Dependence on Planet Masses

Finally, we explore how the scattering results depend on planet masses. When planet masses are unequal, scattering leads to fewer close encounters between planets, since a massive planet can easily eject a low-mass planet. For equal-mass planets, ejection becomes more difficult and planets may suffer a higher number of close encounters, leading to a more efficient eccentricity excitation (Ford & Rasio, 2008). Besides fiducial (with masses $0.5,1,2M_{\rm J}$), we conduct two additional sets of simulations, with all parameters (except for the masses) chosen/sampled identically to fiducial. We consider the case of nearly equal mass planets, with $m_{\rm p}=0.9,1,1.1M_{\rm J}$, and refer to this set of simulations as near-eq-mass. We also consider masses sampled from a log-normal distribution, referred to as lognorm-mass, with mean $1M_{\rm J}$, width $\log_{10}(m_{\rm p}/M_{\rm J})=0.5$, and upper and lower limits $0.25M_{\rm J}$ and $4M_{\rm J}$. These choices are somewhat arbitrary, but qualitatively mimic the observed distribution of giant planet masses. In all simulations, the ordering of the three masses is randomly assigned.

Figure 7 compares the eccentricity distributions of fiducial, near-eq-mass, and lognorm-mass, showing the combined eccentricity distribution of one and two-planet systems (left panel), as well as eccentricities split into one-planet and two-planet systems (middle and right panels). We only show results for simulated one-planet systems with $a<1$AU in order to compare with observed WJ systems in Section 3. The results for the two planet systems show the inner planet eccentricity only. As expected, near-eq-mass produces the highest eccentricities (see also Tables LABEL:table_2p - LABEL:table_1p). Only minor differences are observed between fiducial and lognorm-mass.

The choice of planet masses affects the relative spacing of the two-planet systems, but only marginally. Inspecting Table LABEL:table_2p, the average values of $a_{\rm out}/a_{\rm in}$ for fiducial, near-eq-mass, and lognorm-mass are $2.46$, $2.62$, and $2.3$ respectively. As a result, the finding that in-situ scattering produces very compact two-planet systems is robust with respect to the choice of planet masses.

2.7 Scattering of Four Planets

All scattering experiments discussed thus far began with three unstable planets. We also constructed a set of four-planet systems, with masses $m_{\rm p}=(0.5,1,1.25,2)M_{\rm J}$ (placed in random order), and all other parameters chosen the same as fiducial. We term this set of simulations 4-planets, integrate the systems through “Phase 1” and “Phase 2,” and analyze the properties of the surviving planets. All four-planet systems are destabilized, and the fractions of one, two, and three-planet systems at the end of “Phase 2” are $0.75$, $0.23$, and $0.02$ respectively (see also Table LABEL:table_4p). Adding an additional planet allows more possibilities for ejections. As a result, while fiducial (with three planets) produces roughly equal numbers of one and two-planet systems, 4-planets produces over three times more one-planet systems compared to two-planet systems. The average eccentricity for one-planet systems from 4-planets is $e_{\rm ave}=0.35$ (see also Table LABEL:table_4p), compared with $e_{\rm ave}=0.3$ from fiducial.

Figure 8 shows the two-planet systems produced by 4-planets, along with the fiducial results for reference. The top panel of Fig. 8 shows that 4-planets produces more hierarchical systems, with mean value of $\alpha=a_{\rm out}/a_{\rm in}=4.13$. The remaining panels of Fig. 8 show the inner and outer planet eccentricities. Unsurprisingly, the eccentricities of both planets tend to be higher, with an average inner planet eccentricity $e_{\rm in}^{\rm ave}\simeq 0.2$ for 4-planets compared to $e_{\rm in}^{\rm ave}\simeq 0.15$ for fiducial.

3 Comparison with Observations

The Exoplanet Orbit Database lists 106 RV giant planet systems ($m_{\rm p}\sin i>0.3M_{\rm J}$) with semi-major axes in the range $[0.1{\rm AU}-1]{\rm AU}$, and with eccentricity and mass measurements listed. Of these 106 systems, 23 have an external giant companion with a characterized orbit and minimum mass. The majority of these two-planet systems were previously discussed by Antonini et al. (2016), in the context of high-eccentricity migration, and by Anderson & Lai (2017), in the context of secular perturbations. We augment this sample with the well-known Kepler-419 system (Dawson et al., 2014), which is not flagged by our search criteria due to lack of planet mass listings on exoplanets.org. This yields a sample of 24 WJs with an external giant planet companion.

In the following we compare our scattering results to this observed sample of WJs. Unless otherwise stated, we will refer to our three-planet scattering results.

3.1 Eccentricities

As discussed above, in-situ scattering of three planets leads to a range of eccentricities in the final (post-scattering) planetary orbit(s). The left panel of Fig. 7 presents the combined eccentricity distributions of the one and two-planet systems, consisting of $e$, the eccentricities of single WJs, and $e_{\rm in}$, the inner planet eccentricities of the two-planet systems. We see that in-situ scattering of three planets with a variety of mass choices reproduces the observed distribution fairly well.

However, separating the one and two-planet systems tells a somewhat different story. Recall that the two-planet systems tend to have relatively low eccentricities due to the fact that they formed exclusively via planet-planet collisions (see Table LABEL:table_2p). Comparing the observed and simulated eccentricity distributions in Figs. 3e and 7, we see that although scattering is able to produce the full range of observed eccentricities, there is a clear discrepancy in the shapes of the distributions: the observed distribution is much flatter. This suggests that in-situ scattering of three planets cannot adequately account for the full sample of observed two-planet systems. The addition of a fourth planet alleviates, but does not completely resolve this problem (see the middle panel of Fig. 8).

The single-planet systems, which suffered at least one ejection, tend to have much higher eccentricities. Figs. 4 and 7 reveal that the in-situ scattering of three planets can reproduce the large eccentricities of the observed solitary WJs. But the observed peak at low eccentricities is not reproduced.

The tension between the observed and simulated eccentricity distributions of both the one and two planet systems may be relieved if some of the observed solitary WJs have an undetected companion. If some of the observed one-planet systems are actually two-planet systems, the number of circular one-planet systems may decrease, and the number of circular two-planet systems may increase, leading to a better match with the scattering simulations. To evaluate this possibility, we calculate the radial-velocity semi-amplitude of the outer planet in our simulated two-planet systems from fiducial, given by (Cumming et al., 2008)

[TABLE]

In Fig. 9 we show the RV semi-amplitude versus orbital period of the outer planet in the two-planet systems, where we have assumed an isotropic distribution of sky-projected inclinations $i_{\rm sky}$. Scattering produces a population of outer planets with RV semi-amplitudes mostly in the range $\sim 10-100{\rm m}/{\rm s}$. Due to the compactness of the two-planet systems, many systems have large RV amplitudes and short-orbital periods (with $P_{\rm out}\lesssim 1$ yr) which should be readily detected, but are typically not observed. We note, however, the possibility of undetected companions due to low-cadence RV observations.

To test whether undetected companions to solitary WJs may reconcile the disparities in eccentricity distributions, we apply RVs cuts of $10$ and $100$ m/s to the two-planet fiducial simulations, and classify any two-planet systems as one-planet systems if the outer planet has an RV semi-amplitude less than the specified cut. The results are shown in Fig. 10. A 10 m/s RV cut does not appreciably change the original eccentricity distributions. A 100 m/s RV cut does provide a significantly better agreement with observations. However, since planets with tens of m/s RV amplitudes are readily discovered, there is unlikely to be a large population of undetected companions to observed single WJs with the properties obtained from our scattering simulations. Thus, we conclude that in-situ scattering is unable to separately reproduce the full observed eccentricity distributions of one and two-planet systems, in spite of the apparent agreement in the combined eccentricity distributions (left panel of Fig. 7).

Instead, the eccentricity distribution of one-planet systems is consistent with in-situ scattering having occurred in a significant fraction of systems, alongside a population of low-eccentricity planets that did not undergo scattering. We quantitatively explore this possibility in Section 3.4.

3.2 Spacing and Mutual Inclinations of Two-Planet Systems

As discussed above, scattering of three giant planets usually results in closely spaced systems, with $90\%$ of the fiducial two-planet systems satisfying $a_{\rm out}/a_{\rm in}<3$. Some observed systems are also quite closely spaced (see Fig. 3a), but there exists others that are much more hierarchical. Out of the 24 observed two-planet systems, 7 have $a_{\rm out}/a_{\rm in}<3$.

The mutual inclinations generated via scattering are usually low, with $90\%$ of systems having inclinations less than $12^{\circ}-18^{\circ}$, and the median is $\sim 2^{\circ}-3^{\circ}$. With only a few exceptions, the observed giant exoplanet systems lack constraints on mutual inclinations at present.

3.3 Relative Numbers of One and Two-Planet Systems

As mentioned previously, our observational sample consists of 83 WJs without a detected giant companion, and 24 WJs with a characterized external companion. This yields relative numbers of one and two-planet systems $N_{1,\rm obs}/N_{2,\rm obs}\simeq 3.5$. However, 8 of the 83 solitary WJs have a linear trend in the RV curve, indicating the possibility of an external companion. If we assume the presence of a planetary companion in each of these 8 systems, then $N_{1,\rm obs}/N_{2,\rm obs}=75/32\simeq 2.3$.

The three-planet scattering experiments (i.e. fiducial, lognorm-mass, near-eq-mass, fiducial-K-3, fiducial-K-5) yield relative numbers of single WJs and WJs with external companions $N_{1,\rm sim}/N_{2,\rm sim}\simeq 0.5-1.5$, where the range arises from the choice of planet masses and initial spacing. As a result, in-situ scattering of three giant planets generally predicts an excess of WJs with a companion compared to observations. Four-planet scattering (4-planets) yields $N_{1,\rm sim}/N_{2,\rm sim}\simeq 3.3$, in good agreement with observations.

It is useful to consider the extreme hypothesis that the majority of WJs form in systems of three giant planets, which then undergo scatterings. In this scenario, there are three ways of reconciling the discrepancy between $N_{1,\rm obs}/N_{2,\rm obs}$ and $N_{1,\rm sim}/N_{2,\rm sim}$:

A fraction of the 2-planet systems generated through scattering are actually unstable. As discussed in Section 2.2, $\sim 10\%$ or more of the two-planet systems are expected to eventually undergo instabilities. 2. 2.

A fraction of the single WJs actually have an external companion, most likely having $a_{\rm out}/a_{\rm in}<3$. The possibility of undetected outer companions was already discussed in Section 3.1. 3. 3.

Finally, a third way to resolve the difference in multiplicity frequency is to modify the initial conditions. For practical purposes we have restricted the initial spacings to initial spacings $K=3-5$. Wider spacings may lead to different relative frequencies of collisions and ejections, thus affecting the relative numbers of one and two-planet systems. Such an exploration is beyond the scope of this work.

Alternatively, agreement between the predicted and observed relative frequency of one and two-planet systems is not expected if in-situ-scattering occurred in only a fraction of WJ systems. Indeed, the inability to separately match the eccentricity distributions of one and two-planet systems (see Section 3.1) probably implies in-situ scattering is not a universal outcome of planet formation. We discuss this possibility in the following section.

3.4 Two-Population Mixture Model for Solitary WJ Eccentricities

Figures 4 and 7 show that scattering can account for the moderate ($e\gtrsim 0.1$) to high ($e\sim 0.9$) eccentricities of solitary WJs, but fails to reproduce the peak at low eccentricities. This may indicate two populations of WJs. Here we present a mixture model for the solitary WJ eccentricities, augmenting our sample of one-planet systems obtained from fiducial with a population of low-eccentricity planets (that did not undergo scattering). For simplicity, we neglect the uncertainties of observed eccentricities.

Consider a mixture model such that the probability density function for eccentricities takes the form

[TABLE]

where $f_{\rm circ}$ is the injected probability density function of low-eccentricity planets, and $f_{\rm scat}$ is that obtained from fiducial. The free parameter $\mathcal{F}=[0,1]$ quantifies the low-eccentricity planet fraction.

We adopt a half-Gaussian for $f_{\rm circ}(e)$ peaked at $e=0$ and truncated at $e=1$, and explore a range of characteristic widths $\sigma$, in the range $[0-0.2]$. This accounts for the fact that modest eccentricities may be generated by factors other than scattering, e.g. by planet-disk interactions. To construct a smooth function for $f_{\rm scat}(e)$, we apply a Gaussian Kernel-density estimate to our simulated eccentricities. Armed with $f_{\rm circ}(e)$ and $f_{\rm scat}(e)$, we explore a grid of values for the parameters $(\mathcal{F},\sigma)$, and calculate the likelihood function for obtaining the observed eccentricities $e_{i}$,

[TABLE]

with $n_{\rm obs}=83$. Maximizing $\mathcal{L}$ yields the parameter estimates $\mathcal{F}\approx 0.35$ and $\sigma\approx 0.07$. This estimate, alongside $68\%$ and $95\%$ contours is shown in Fig. 11. Taking the $95\%$ contour as the uncertainty, we find that $\mathcal{F}$ lies in the range $0.18-0.54$.

In other words, the observed eccentricity distribution of solitary WJs is consistent with scattering having occurred in roughly half or more of systems. This finding also holds for the 4-planets scattering simulations: Performing an identical analysis for the 4-planets scattering results, we obtain similar parameter estimates, with $\mathcal{F}\approx 0.39$ and $\sigma\approx 0.07$.

4 Conclusion

In this paper we have undertaken a comprehensive study of giant planet scattering at sub-AU distances ($0.1-1$ AU) from the host star. Our initial system of multiple giant planets in nearly circular orbits is consistent with WJs that formed either in-situ or by disk migration without substantial eccentricity excitation. Our goal is to catalogue the scattering pathways and outcomes, and to identify the extent to which scattering may have contributed to the observed population of eccentric WJs. We focus primarily on systems of three initially unstable planets, but also briefly consider systems of four planets.

4.1 Key Results

• •

For our explored parameter space (see Table LABEL:table1), three-planet scatterings result in roughly equal proportions of one-planet and two-planet systems. These are produced through a combination of collisions and ejections (see Fig. 2). The first (or only) planet “loss” is a planet-planet collision $60\%-80\%$ of the time. For systems that later undergo a second planet loss, ejections occur over $90\%$ of the time. Two-planet systems are produced almost entirely from planet-planet collisions, whereas one-planet systems arise from a combination of collisions and ejections. Branching ratios separated into two-planet and one-planet systems are depicted in Tables LABEL:table_2p and LABEL:table_1p.

• •

The outcomes of three-planet scattering depend on the semi-major axes of the initial systems (Fig. 5). The fraction of single planet systems produced ($f_{\rm 1p}$) increases from $\sim 38\%$ to $\sim 60\%$ as the initial $a_{1}$ increases from $0.1$ AU to $1$ AU, with a corresponding increase in the average planetary eccentricities.

• •

Scattering tends to produce closely spaced two-planet systems. For systems of initially three unstable planets, $90\%$ of the surviving two-planet systems have a semi-major axis ratio $a_{\rm out}/a_{\rm in}\lesssim 3$. For systems of initially four unstable planets, $90\%$ of two-planet systems have $a_{\rm out}/a_{\rm in}\lesssim 5$.

• •

The combined eccentricity distribution of all WJ systems (consisting of single WJs and the inner planet of the two-planet systems) produced by scattering agrees well with the observed distribution (see the left panel of Fig. 7). However, splitting the results into one and two-planet systems (see the middle and right panels of Fig. 7), we find some discrepancies between the simulations and observations: (i) Examining the two-planet systems, we find that scattering produces too many low-eccentricity planets, due to the fact that all two-planet systems arose from planet-planet collisions. This is inconsistent with the observed, much flatter eccentricity distribution. This discrepancy is alleviated, but not fully resolved by addition of a fourth planet to the initial scattering setup (see Fig. 8). (ii) Examining the one-planet systems produced by scattering and comparing with observed solitary WJs (those lacking detection of a giant planet companion), we find that scattering well reproduces the substantial tail of modest and high eccentricities, but does not reproduce the observed peak at low eccentricities.

• •

The observed eccentricity distribution of solitary WJs can be reproduced by two populations of planets, with roughly $60\%$ having undergone violent scattering, and the remaining having experienced a quiescent dynamical history (see Section 3.4 and Fig. 11).

• •

Planet-planet collisions have a range of impact parameters, with grazing collisions dominating over head-on (see Fig. 12).

4.2 Discussion

Our study shows that in-situ scattering of multiple giant planets is effective in producing eccentric WJs. However, while scattering provides a good match to the entire observed WJ eccentricity distribution (consisting of both one and two-planet systems), it yields an excess of circular two-planet systems, and a dearth of circular one-planet systems (see Fig. 7). A possible resolution to this discrepancy is if many observed solitary WJs have an undetected outer companion, thereby shifting the excess of circular planets from the two-planet systems to the one-planet systems. However, given the large RV amplitudes of the outer planets obtained from our scattering experiments ($10-100$ m/s, see Fig. 9), such companions should be readily detectible. Thus, observed solitary WJs are unlikely to have enough undetected outer companions to resolve the discrepancy in eccentricity distributions (see Fig. 10).

Besides the disagreement in eccentricity distributions, the spacings of two-planet systems generated through in-situ scattering are also inconsistent with many of the observed two-planet systems. Over $90\%$ of the two-planet systems produced from the fiducial simulations have semi-major axis ratio $a_{\rm out}/a_{\rm in}\lesssim 3$, whereas many observed systems are more hierarchical. The tendency for scattering to result in compact two-planet systems was pointed out by Raymond et al. (2009). Many of the observed eccentric two-planet systems remain challenging to explain: They are too compact to have originated from a traditional high-eccentricity migration scenario (Antonini et al., 2016), yet are too hierarchical to have formed via in-situ scattering. Highly eccentric, hierarchical systems (with $a_{\rm out}/a_{\rm in}\gtrsim 10$) could in principle have originated from a scattering event which initially occurred at several AU distances or beyond, as explored by Mustill et al. (2017). A thorough study of the parameter space for initial systems consisting of a single WJ and two or more unstable planets much farther out is an important, but computationally expensive problem.

Our study (Section 3.4) indicates that a two-population model well reproduces the observed eccentricity distribution of solitary WJs, consisting of a $\sim 60\%$ scattered component and $\sim 40\%$ quiescent component. The quiescent component is consistent with the finding by Huang et al. (2016) that a substantial fraction ($\sim 50\%$) of WJs have low-mass neighbors (and therefore have experienced no “dynamical violence”). Two populations of WJs have been previously proposed by Dawson & Murray-Clay (2013), consisting of low-eccentricity systems around low-metallicity stars, and higher eccentricity systems around high-metallicity stars. Scattering is consistent with such a metallicity trend, under the expectation that multiple giant planets on unstable orbits form more easily around higher metallicity stars. Our success in reproducing many of the WJ eccentricities with planet-planet scattering (see also Marzari & Nagasawa, 2019; Frelikh et al., 2019) indicates that formation of multiple giant planets at sub-AU distances, or arrival at sub-AU orbits due to disk migration in a compact configuration, may be common. However, extrapolating our results to the HJ population would imply an inconsistency with observations. Planet scattering interior to $\sim 0.1$ AU would result in more planet-planet collisions (see Fig. 5 and also Petrovich et al. 2014), leading to an even larger fraction of compact two-planet systems. Since HJs lack close, massive companions, scattering of three or more giant planets within $\sim 0.1$ AU is in contradiction with observations (although scattering of only two planets would remain viable). This may imply that HJs and WJs often have distinct formation and dynamical histories.

Our three-planet scattering simulations show that similar numbers of one and two-planet systems are produced, while our limited four-planet scattering experiments suggest the resulting one-planet systems could outnumber the two-planet systems by a factor of three. Only a subset of the observed two-planet systems are compact enough to be consistent with the scattering predictions from this paper (7 out of 24 systems have $a_{\rm out}/a_{\rm in}<3$, see also Section 3). Given the observed sample of 83 solitary WJs, the “scattered + quiescent” model would predict a somewhat larger number of compact two-planet systems than is currently observed.

Better consistency between the observed relative planet multiplicities and those obtained from the “scattering + quiescent” model may be obtained if some solitary WJs have a close giant planet companion, or if some of the simulated two-planet systems later become unstable. As discussed in Section 3.1, the RV amplitudes of the companions obtained from scattering are large, so the second explanation may be more favorable. The relative number of one and two-planet systems will increase with time, as some of the two-planet systems undergo instabilities. Follow-up long-term integrations of a small subset of the two-planet systems (spanning $10^{9}$ orbits of the inner planet) indicated that at least $\sim 10\%$ of the two-planet systems will eventually go unstable, but the actual number may be higher over Gyr timescales. Furthermore, physical effects not included in our scattering calculations (e.g. perturbations from a distant companion) may act to destabilize some of the compact two-planet systems.

Acknowledgments

KRA is grateful to Tom Loredo for assistance with the calculation presented in Section 3.4, and thanks Ruth Murray-Clay, Eric Ford, and Renata Frelikh for useful discussions. We also thank the anonymous referee for insightful feedback and comments. This work has been supported in part by NASA grants NNX14AG94G and NNX14AP31G, and NSF grant AST-1715246. KRA and BP were supported by NASA Earth & Space Science Fellowships in Planetary Science. This research has made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org. Simulations in this paper made use of the rebound code which can be downloaded freely at http://github.com/hannorein/rebound.

Appendix A Treatment of Planet-Planet Collisions

As in previous N-body works, this paper has assumed completely inelastic collisions, so that once the distance between two planets becomes less than the sum of their physical radii, the planets merge conserving mass and momentum. Such a treatment is clearly a simplification. In reality, giant planets will tidally interact upon close encounters, possibly leading to mass loss. In addition, head-on versus grazing collisions may result in different outcomes. Although expanding on this treatment is beyond the scope of this paper, here we discuss some of the collisional dynamics and results from the fiducial simulations. These results will be useful as initial conditions for future hydrodynamical simulations of planet-planet collisions.

We define the collision impact parameter, $b_{\rm coll}=|{\bf r}\times\hat{{\bf v}}|$ where ${\bf r}$ is the relative distance between the centers of mass of the two planets, and $\hat{{\bf v}}$ is the relative velocity unit vector. The top left panel of Fig. 12 shows that many planet-planet collisions are grazing, with $b_{\rm coll}/(R_{1}+R_{2})\sim 1$, where $R_{1}$ and $R_{2}$ are the planetary radii.

Next we discuss some aspects of the collision of two planets $m_{1}$ and $m_{2}$. Once the planet-planet gravitational interaction becomes sufficiently strong so that the tidal gravity from the host star can be ignored, the scattering process up to just before merging can be modeled approximately as an isolated two-body encounter. Denote $v_{\rm\infty}$ as the relative velocity of two planets “at infinity”, and $v_{\rm coll}$ as the relative velocity just before collision. Energy conservation yields

[TABLE]

where $v_{0}^{2}=2GM/(R_{1}+R_{2})$, with $M$ the sum of the planet masses. Angular momentum conservations yields

[TABLE]

where $b_{\infty}$ is the impact parameter at “infinity”.

The upper right panel of Fig. 12 shows $v_{\rm coll}/v_{0}$. As expected, the majority of collisions have $v_{\rm coll}/v_{0}\sim 1$. The lower left panel of Fig. 12 shows $v_{\infty}^{2}/v_{0}^{2}$. Typically $v_{\infty}^{2}/v_{0}^{2}\ll 1$. Negative values of $v_{\infty}^{2}$ indicate that the two planets approach each other on an elliptical orbit. Finally, the lower right panel of Fig. 12 shows $b_{\infty}$ in units of the Hill radius, $R_{\rm H}=a(M/3M_{\star})^{1/3}$.