New Insights on Planet Formation in WASP-47 from a Simultaneous Analysis of Radial Velocities and Transit Timing Variations

We used the publicly available N-body integrator TTVFast (Deck et al. 2014) with the python wrapper ttvfast-python ¹³ to forward-model the transit times of the inner three planets. Unlike in the B15 analysis, we allowed all of the initial osculating elements, particularly the orbital periods and initial times of transits, to vary. Thus, the variables for each planet k are: the mass of the planet M_k, the orbital period P_k, the first time of transit tt_k, and the eccentricity parametrization $\sqrt{{e}_{k}}$ cos${\omega }_{k}$, $\sqrt{{e}_{k}}$ sin${\omega }_{k}$. We limited $e\lt 0.06$ for the three inner planets, in accordance with the 10 Myr stability analysis in B15. For each adjacent pair of planets, we satisfied the Hill criteria for stability for low-eccentricity orbits, as described in Equations (24) and (28) of Gladman (1993):

$$\begin{eqnarray}&&p-q\gt 2.4{({\mu }_{1}^{2}+{\mu }_{2}^{2})}^{1/3},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&p-q\gt \sqrt{8/3({e}_{1}^{2}+{e}_{2}^{2})+9\times \max {({\mu }_{1},{\mu }_{2})}^{2/3}},\end{eqnarray} \tag{ 2 }$$

where the subscript 1 refers to the inner planet and 2 refers to the outer planet, q is the apoapse distance of the inner planet $q=1+{e}_{1}$, p is the periapse distance of the outer planet $p=(1-{e}_{2})\tfrac{{a}_{2}}{{a}_{1}}$, e is the eccentricity, and μ is the planet-to-star mass ratio. We also required the planet masses and orbital periods to have positive values. Because the orbits of the planets are very nearly coplanar (Becker et al. 2015; Almenara et al. 2016), we fixed the orbital inclination and longitude of the ascending node in a manner consistent with coplanar, edge-on orbits. See Table 1 for a summary of the priors and constraints.

Table 1.  Priors on Dynamical Parameters

Parameter	Priors
M	M > 0, Hill criterion
P	P > 0, Hill criterion
TT	None
$\sqrt{e}$ cosω	$e\lt 0.06$, Hill criterion
$\sqrt{e}$ sinω	$e\lt 0.06$, Hill criterion
$\sqrt{{e}_{c}}$cos${\omega }_{c}$	$e\lt 1$, Hill criterion
$\sqrt{{e}_{c}}$sin${\omega }_{c}$	$e\lt 1$, Hill criterion

Download table as:  ASCIITypeset image

Table 2.  Dynamical Parameters from the Best N-body Fit to TTVs Only (TTVFast)

Parameter	Median ± Std. Dev.	Units
Me	176 ± 118	${M}_{\oplus }$
Mb	549 ± 252	${M}_{\oplus }$
Md	16.1 ± 3.8	${M}_{\oplus }$
Pe	0.78964 ± 0.00002	days
Pb	4.150 ± 0.006	days
Pd	9.12 ± 0.05	days
TTe	2146.7639 ± 0.0008	BJD–2454833
TTb	2149.969 ± 0.006	BJD–2454833
TTd	2155.40 ± 0.05	BJD–2454833
$\sqrt{{e}_{e}}$cos${\omega }_{e}$	−0.006 ± 0.14	⋯
$\sqrt{{e}_{e}}$sin${\omega }_{e}$	−0.008 ± 0.14	⋯
$\sqrt{{e}_{b}}$cos${\omega }_{b}$	0.001 ± 0.08	⋯
$\sqrt{{e}_{b}}$sin${\omega }_{b}$	0.006 ± 0.07	⋯
$\sqrt{{e}_{d}}$cos${\omega }_{d}$	−0.09 ± 0.11	⋯
$\sqrt{{e}_{d}}$sin${\omega }_{d}$	0.07 ± 0.09	⋯

Note. These are the initial astrocentric Keplerian orbital elements, reported at epoch BJD 2456979.4961. They are not the time-averaged orbital properties of the planets.

Download table as:  ASCIITypeset image

We used the Markov Chain Monte Carlo (MCMC) Python package emcee (Foreman-Mackey et al. 2013) to explore the posteriors of various combinations of the dynamical parameters. We ran 60 walkers 5 × 10⁵ steps each, throwing away the first 10⁵ steps as burn-in, and checked that the multivariate extension of the potential scale reduction factor (PSRF) statistic ($\hat{R}\lt 1.2$, Gelman & Rubin 1992; Brooks & Gelman 1998) converged. We also inspected the chains by eye to check for convergence. Our best fit¹⁴ to the transit times using TTVFast is shown in Figure 1. Table 2 summarizes our results from this N-body fit to the TTVs.

Zoom In Zoom Out Reset image size

We find that for planets e and b, the masses and eccentricities of the planets are highly covariant with the initial osculating orbital period and the initial transit time of neighboring planets (see Figure 2). This is because the initial osculating orbital period is translated to an instantaneous velocity and acceleration, and the planet's acceleration depends on the mass and position of its neighbor. Therefore, it is critical to allow the initial orbital periods and times of transit of all the planets to vary in order to explore the full range of possible planet masses. Furthermore, only one super-period of the TTVs is observed, so an average orbital period and a representative time of transit are not as well determined for planets b and d as they might be with the observation of multiple super-periods. Thus, we find that the TTVs do not constrain the masses of WASP-47 e or WASP-47 b as tightly as what is reported in B15. Whereas B15 find ${M}_{b}={341}_{-55}^{+73}\,{M}_{\oplus }$, we find ${M}_{b}=549\pm 252\,{M}_{\oplus }$, using the exact same TTV measurements. The mass for planet e reported in B15 stems from their choice of prior: they find ${M}_{e}\lt 22\,{M}_{\oplus }$, whereas we find ${M}_{e}=176\pm 118\,{M}_{\oplus }$. However, the TTVs of planet b place strong constraints on the mass of planet d. B15 find ${M}_{d}=15.2\pm 7\,{M}_{\oplus }$, and we find ${M}_{d}=16.1\pm 3.8\,{M}_{\oplus }$. While an analytic covariance between planet masses and the free eccentricity exists (Lithwick et al. 2012), our stability constraint that $e\lt 0.06$ minimizes the effects of this degeneracy.

Zoom In Zoom Out Reset image size

We used the publicly available analytic TTV package TTVFaster to model the transit times of the inner three planets observed in B15. The TTVFaster code analytically transforms the planet mass M_k and average orbital elements P_k, tt_k, e_kcos${\omega }_{k}$, e_ksin${\omega }_{k}$ ¹⁵ into a TTV pattern, to first order in eccentricity, to a user-specified precision in the disturbing function. We found that modeling to sixth order in the expansion of the Laplace coefficient ${b}_{1/2}^{j}$ ($j=\{0,1,2,3,4,5,6\}$), (Murray & Dermott 2000) was sufficient to reproduce the observed TTV signature with the same fidelity as that produced in the N-body analysis (see Figure 1). In general, TTVFaster is designed to work for planets that (1) are not extremely close to a mean motion resonance, (2) have low eccentricities, or (3) have low masses. Because WASP-47 b is a Jupiter-mass planet, we wanted to see if TTVFaster modeled the orbital dynamics correctly.

Incorporating the priors from Table 1, we used Python packages lmfit (Newville et al. 2014) and emcee to explore the posteriors of the dynamical parameters. We ran 100 walkers 20,000 steps, throwing away the first 4000 steps as burn-in. We note that the chains converged much faster (according to the PSRF statistic) when we used TTVFaster than when we used TTVFast, because the TTVs provide better constraints on the average orbital parameters than they do on the initial orbital parameters.

The mass and eccentricity distributions we determined from the analytic solution to the observed TTVs are consistent with the N-body model (Figure 3). Since fitting an analytic model to the TTVs is orders of magnitude faster than a full N-body analysis (especially when the long time baseline for RVs is required), we use the analytic modeling technique for the rest of this paper.

Zoom In Zoom Out Reset image size

Table 4 summarizes the relative information in various dynamical analyses of the planet masses and eccentricities. B15 modeled the K2 TTVs in a three-planet N-body analysis in which the orbital periods and initial times of transit were fixed, resulting in narrow posteriors for the mass of planet b. The mass constraint for planet e comes from the choice of prior, rather than the TTVs. B15 also forward-modeled the system for 10 Myr using Mercury (Chambers 1999) to ensure stability, which resulted in the tight eccentricity constraints for the inner planets: ${e}_{k}\lt 0.06$.

Table 4.  Relative Information in Literature Dynamical Analyses

Parameter	B15	A16.1	A16.2	S16	TTVs–Nbody	RVs+TTVs
${M}_{\star }[{M}_{\odot }]$	1.04 ± 0.08	${1.11}_{-0.49}^{+0.89}$	1.029 ± 0.031	0.99 ± 0.05	0.99 ± 0.05	1.00 ± 0.05
Me [${M}_{\oplus }$]	$\lt {22}^{P}$	${9.1}_{-2.9}^{+5.5}$	${9.1}_{-2.9}^{+1.8}$	9.11 ± 1.17	176 ± 118	9.1 ± 1.0
Mb [${M}_{\oplus }$]	${341}_{-55}^{+73}$	${383}_{-120}^{+190}$	363.8 ± 8.6	356 ± 12	549 ± 252	358 ± 12
Md [${M}_{\oplus }$]	15.2 ± 7	${16.8}_{-7}^{+12}$	15.7 ± 1.1	12.75 ± 2.70	16.1 ± 3.8	13.6 ± 2.0
Mc [${M}_{\oplus }$]	⋯	${500}_{-190}^{+320}$	${470}_{-100}^{+200}$	411 ± 18	⋯	416 ± 16
ee	$\lt 0.06$	$\lt 0.11$	⋯	$={0}^{P}$	$\lt {0.06}^{P}$	$\lt {0.06}^{P}$
eb	$\lt 0.06$	$\lt 0.01$	⋯	$\lt 0.013$	$\lt 0.05$	$\lt 0.01$
ed	$\lt 0.06$	$\lt 0.024$	⋯	$={0}^{P}$	$\lt 0.044$	$\lt 0.025$
ec	⋯	0.36 ± 0.12	⋯	0.27 ± 0.04	⋯	0.28 ± 0.02

Note. B15—Becker et al. (2015), K2 TTVs, fixed P, TT for each planet, 10 Myr stability-enforced. A16.1—Almenara et al. (2016), photodynamical analysis of K2 TTVs and 71 RVs. A16.2—including stellar models. S16—Sinukoff et al. 2016, 118 RVs. TTVs–Nbody—TTVFast N-body analysis (presented herein), K2 TTVs. RVs+TTVs—simultaneous analysis of K2 TTVs and 118 RVs. All upper limits are at 95% confidence. P—results come from a prior.

Download table as:  ASCIITypeset image

Almenara et al. (2016, A16 hereafter) performed a photodynamical analysis of the K2 TTVs and 71 literature RVs from the PFS and CORALIE spectrographs. The high signal-to-noise of the transits of planet b allowed them to determine stellar-limb-darkening parameters and the planet impact parameter very precisely, which, in combination with the small eccentricity of planet b, led to a very precise determination of the stellar density through asterodensity profiling (Kipping 2014). This led to a model-independent estimate of the planetary masses, presented in column A16.1 of Table 4. By including constraints from the Dartmouth stellar isochrone models (Dotter et al. 2008), the authors were able to constrain the star and planet masses more precisely, but at the expense of accuracy. The model-dependent star and planet masses are shown in column A16.2.

S17 combined 47 new HIRES RVs with the 71 literature RVs. While the CORALIE RVs had provided a long enough baseline to detect the non-transiting planet c (Neveu-VanMalle et al. 2016) and the PFS RVs had provided sufficiently high precision to detect a marginal RV signal from planet e (Dai et al. 2015), the HIRES RVs provided both a long baseline and high precision in a single data set. The HIRES data combined with the other RV data sets resulted in smaller uncertainties for all the planet masses than what had been reported in previous RV studies.

Our TTV-only N-body analysis (TTVs–Nbody) and simultaneous RV and TTV analysis (RVs+TTVs) are shown in Table 4. The TTVs–Nbody column illustrates how much information is contained in the TTVs. The RVs+TTVs column illustrates how much information is gained from a joint analysis of the TTVs and RVs. Our RVs+TTVs analysis confirms the precise values obtained by A16 when they include constraints from stellar models (A16.2).

How much information about planet masses comes from the TTVs? As discussed in the TTV-only analysis (see Section 3), the TTVs do not provide much information about the transiting planet masses, with the exception of the mass of planet d, which is constrained through the TTVs of planet b.

Simultaneously modeling the TTVs and RVs of planet d yields a more precise determination of the mass of planet d than can be obtained from either analysis alone: the uncertainties shrink from $3.8\,{M}_{\oplus }$ (TTVs) and $2.7\,{M}_{\oplus }$ (RVs) to $2.0\,{M}_{\oplus }$ (TTVs+RVs). The TTVs provide no information about the mass of planet c, which has a very long period compared to the inner planetary system and thus has no effect on the TTVs. Thus, the RVs provide the majority of the information about planet masses, although the TTVs contribute substantially to the mass measurement of planet d.

Information about planet eccentricities comes from stability constraints, the TTVs, and the RVs. The eccentricity of planet e is not constrained by either the TTVs or the RVs, so its eccentricity varies from 0 to 0.06 (the upper limit from stability requirements). The RVs constrain ${e}_{b}\lt 0.013$ (95% confidence, S17). The TTVs constrain ${e}_{b}\lt 0.02$ (95% confidence), and the combined RVs+TTVs further constrain ${e}_{b}\lt 0.011$ (95% confidence). The RVs alone do not provide a strong constraint for the eccentricity of planet d (S17 fixed e_d = 0). The TTVs alone constrain ${e}_{d}\lt 0.05$ (95% confidence), and the combined TTVs+RVs constrain ${e}_{d}\lt 0.025$ (95% confidence). Thus, the TTVs provide additional information about the small eccentricities of planets b and d. The eccentricity of planet c is determined entirely from RVs because the planet is dynamically decoupled from the inner planetary system.

Thus, combining the TTVs and RVs provides more information about masses and eccentricities than either data set does alone. We discuss the physical importance of having precise measurements for the masses and eccentricities below.

WASP-47 is unusual in that the masses of its planets span almost two orders of magnitude. The low-mass planets e and d are shown in density–radius and mass–radius plots for small planets (see Figure 8). The lowest-mass planet, WASP-47 e, is $9.1\pm 1.0\,{M}_{\oplus }$. At $1.8\,{R}_{\oplus }$, this planet has a density of $9.22\pm 1.06\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. Planets larger than approximately $1.5\,{R}_{\oplus }$ are unlikely to be rocky (Weiss & Marcy 2014; Rogers 2015). Yet, planet e is small and dense enough that a rocky composition is likely based on an extrapolation of the empirical relationship for rocky planets smaller than 1.5 Earth radii (Weiss & Marcy 2014) and theoretical predictions of planet mass and radius for an Earth-like composition (Seager et al. 2007). However, the density of planet e is also consistent with slightly lower densities that might correspond to a rocky interior overlaid with a thin, low-mass envelope of high mean molecular weight materials. Such a composition has also been hypothesized for 55 Cancri e, which has a very similar mass, radius, and bulk density to WASP-47 e (Lopez 2016, S17). Like 55 Cnc e, WASP-47 e is also an ultra-short period planet cohabiting an orbital system with giant planets, which underscores the question of how planetary system architecture and planet compositions are related.

Zoom In Zoom Out Reset image size

By contrast, WASP-47 d, which is $3.6\,{R}_{\oplus }$, has a mass of $13.6\pm 2.0\,{M}_{\oplus }$. This is a slightly higher mass than what was reported in S17 because the TTVs add mass information. The additional information from the TTVs also narrows the mass posterior, shrinking the uncertainty from 2.7 to $2.0\,{M}_{\oplus }$. The density of WASP-47 d is $1.63\pm 0.23\,{\rm{g}}\,{\mathrm{cm}}^{-3}$, making it a high-density member of the population of sub-Neptune-sized planets with volatile envelopes. The mass–radius diagram in Figure 8 shows that the mass, radius, and density of WASP-47 d make it one of the most Neptune-like planets discovered to date. While its composition could be explained as a two-layer model of a H/He envelope atop a silicate-iron core, a Neptune-like composition that includes a thick layer of super-ionic water might also explain the bulk properties of WASP-47 d.

WASP-47 b is a Jupiter-mass planet ($358\pm 12\,{M}_{\oplus }$) that receives 440 ± 70 times as much incident stellar irradiation as the Earth does. At $13.11\pm 0.89\,{R}_{\oplus }$, the planet has a typical density ($1.02\pm 0.02\,{\rm{g}}\,{\mathrm{cm}}^{-3}$) for its mass and incident stellar flux (see Figure 9), consistent with various theories (Fortney & Nettelmann 2010; Batygin et al. 2011, and references therein) that stellar irradiation inflates the planet and/or prevents the planet from cooling.

Zoom In Zoom Out Reset image size

The orbits of the three inner planets are profoundly circular. B15 found that eccentricities of $\lt 0.06$ were required for stability. Here, we tighten the eccentricities to ${e}_{b}\lt 0.011$ and ${e}_{d}\lt 0.025$ (95% confidence). In the highest-eccentricity cases for planets b and d, they tend to be apsidally aligned. We compute ${e}_{d}\cos {\omega }_{d}-{e}_{b}\cos {\omega }_{b}=-0.001\pm 0.005$ and ${e}_{d}\sin {\omega }_{d}\,-{e}_{b}\sin {\omega }_{b}=0.0\pm 0.007$.

Although the tidal circularization timescale for planet e is only $\sim {10}^{4}-{10}^{5}$ years, depending on the tidal Q value for the planet, our N-body analysis revealed that the neighboring giant planet (b) perturbs the eccentricity of planet e on a timescale of 1.26 days. This timescale happens to be related to the orbital periods of both e and b by $1/{P}_{\mathrm{kick}}={(1/{P}_{e}+2/{P}_{b})}^{-1}$. Since there is no commensurability between the orbital periods of b and e, we expect that over long timescales, the kicks from planet b average out, so the argument of periastron of planet e should not have a preferred value. However, planet e could still have an average eccentricity that is higher than zero. Likewise, planet d should not be assumed to have zero eccentricity because it exhibits TTVs. Because planet d is near the 2:1 resonance with planet b, its argument of periastron circulates at the frequency of the TTV super-period.

The very low eccentricities of the WASP-47 planets are remarkably like those of the solar system planets (see Figure 10). The average eccentricity of the detected planets in WASP-47 is $\lt 0.09;$ in the solar system, the average eccentricity of the planets and Pluto is 0.08.

Zoom In Zoom Out Reset image size

Compact systems like the WASP-47 inner planets might form in situ from protoplanetary disks that are more massive than the minimum mass solar nebula (Chiang & Laughlin 2013). Scenarios in which either gas-poor sub-Neptunes or gas-rich Jupiters form in situ have been proposed (Lee et al. 2014; Lee & Chiang 2015; Batygin et al. 2016; Lee & Chiang 2016). If all the WASP-47 planets formed in situ from the same nebular material, the challenge is to explain how WASP-47 b managed to achieve runaway growth, whereas its immediate neighbors remained gas-poor. Whether a Jupiter-mass or sub-Saturn mass planet forms depends on the core mass, atmospheric opacity, and the disk lifetime (Hori & Ikoma 2011; Venturini et al. 2015). Cores of mass $\sim 5-15{M}_{\oplus }$ and larger undergo runaway gas accretion (with the exact value of critical core mass depending on the atmospheric opacities, metallicities and the disk gas dissipation timescale; Batygin et al. 2016; Lee & Chiang 2016). However, gas-poor sub-Neptunes are formed instead of giant planets if the cores form and accrete their envelopes in a gas-poor (transitional) disk. The formation of sub-Neptune cores by giant impact requires at least four orders of magnitude of gas-depletion with respect to the minimum mass extrasolar nebula (Lee & Chiang 2016), leaving $\sim 0.3\,{M}_{\oplus }$ of gas in the disk. This gas budget is insufficient to form the $\sim 300\,{M}_{\oplus }$ of gas in WASP-47 b.

A modification of in situ formation is inside-out growth via pebble accretion, in which pebbles are transported inward through the disk until they are stopped by a pressure maximum at the boundary between the magneto-rotational instability (MRI) zone and the magnetic dead zone (Chatterjee & Tan 2014). At this boundary, the pebbles may become Toomre unstable or may coalesce via core accretion. As the gas disk clears and the boundary between the MRI and dead zone moves outward, the site of planet formation gradually moves out through the disk. Because the hot Jupiter is situated between two sub-Neptunes, the accretion rate of the disk would likely need to increase by roughly an order of magnitude, and then decrease again, to explain the high mass of the hot Jupiter planet compared to its neighbors. Furthermore, the efficacy of forming hot Jupiters via inside-out planet formation has not been well studied.

Alternatively, the planets could have formed elsewhere in the disk and then migrated via interactions with the disk to their present locations. In both Type I and Type II migration, the gas disk damps planet eccentricity, allowing planets to maintain circular orbits in a manner consistent with the nearly circular orbits of the three inner planets. However, slow migration can trap the planets in mean motion resonances. While WASP-47 b and d are near the 2:1 mean motion resonance, they are not trapped. Thus, their migration history must either include a mechanism to prevent planets b and d from entering the 2:1 resonance, or remove them from the resonance (e.g., Adams et al. 2008; Goldreich & Schlichting 2014). In particular, Deck & Batygin (2015) find that if the inner planet of a pair near a first-order mean motion resonance is more massive (as is the case for WASP-47 b and d), escape from resonance is unlikely. Furthermore, migration in the disk does not explain the eccentricity of planet c ($0.28\pm 0.02$. Also, planets e and d would need to migrate through the disk without accreting gas.

Planet–planet scattering and Kozai–Lidov oscillations are big-body (as opposed to gas-and-dust) migration mechanisms. After the gas disk has dissipated, these mechanisms can introduce moderate to high eccentricities in the orbits, potentially accounting for the moderate eccentricity of planet c. Kozai–Lidov oscillations are initiated only when two planets have mutual inclinations of at least 40°. The planets swap angular momentum, causing dramatic variations in the inclinations and eccentricities of the planets over time (Kozai 1962; Lidov 1962). The strongest lines of evidence for past Kozai–Lidov interactions would be (1) an observed large mutual inclination between planet c and the inner solar system, and/or (2) a non-zero obliquity between planet b and the stellar spin axis. Sanchis-Ojeda et al. (2015) ruled out a significant spin–orbit misalignment for planet b, but the degeneracy between stellar $v\sin i$ and the planet-star obliquity allows moderate misalignments ($| \lambda | \lt 48^\circ$ with 1σ confidence). If planet c is coplanar and the spin–orbit alignments of the planets are small, a coplanar, high-eccentricity migration (HEM) scenario such as the one proposed in Petrovich (2015) might best explain the present locations of the planets. However, the current nearly circular, coplanar orbits of the inner planetary system are not consistent with a past modest eccentricity and/or inclination for planet b. Recall from the stability analysis in B15 that the nearly all of the scenarios in which the eccentricity of planet b exceeds 0.06 become unstable within 10 Myr. Therefore, although high-eccentricity migration mechanisms can explain the present positions of planets b and c, such mechanisms would have destroyed the small, close-in planets. Indeed, simulations of giant planets migrating inward find that the giant planets collide with low-mass inner planets in both high-eccentricity (Mustill et al. 2015) and low-eccentricity (Batygin & Laughlin 2015) scenarios.

One way to explain the present orbits and compositions of the WASP-47 planets is that the planets did not all form at the same time. As a point of reference, consider the solar system. In our own solar system, the giant planets must have formed early (within 1 to 10 Myr), when the protoplanetary disk was still gas-rich (de Pater & Lissauer 2001). In contrast, the formation of terrestrial-mass planets by planetesimal accretion can take as long as ∼10⁸ years in N-body simulations and is more efficient in gas-poor disks where the planetesimal eccentricities can grow, leading to more collisions (Lissauer 1987; Pollack et al. 1996; Ida & Lin 2004; Lee & Chiang 2016). Thus, it is possible that the solar system giant planets formed and migrated to their present locations before the terrestrial planets formed. The Nice Model demonstrates that early formation and migration of Jupiter and Saturn can reproduce the current orbital architectures of the gas giants and ice giants while also accounting for the period of Late Heavy Bombardment on the terrestrial planets, and Trojan asteroids (Gomes et al. 2005; Morbidelli et al. 2005; Tsiganis et al. 2005). The Grand Tack Model shows that a reversal in the direction of Jupiter's migration due to torques from Saturn can explain the small size of Mars and detailed compositional features of the asteroids (Walsh et al. 2011; Morbidelli et al. 2012). Both of these models feature two stages of planet formation:

• 1.  
  Early giant planet growth, combined with disk and/or planet-induced migration, produces the current compositions of giant planets and sculpts the radial distribution of planetesimals.
• 2.  
  When or after the gas disk clears, planetesimal accretion proceeds in the sculpted planetesimal disk, resulting in the formation of low-mass, predominantly rocky planets.

Likewise, the formation and evolution of the planets in WASP-47 might be best explained by a two-stage process. The giant planets in the WASP-47 system might have formed early in the gas-rich disk and exchanged energy with each other, additional bodies, or the disk to migrate to their present locations. If the gas disk cleared while WASP-47 b and c were finishing their migration, the final years of the giant planets' migration could have influenced the solid material in the disk, exciting protoplanetary solids to higher eccentricities and inducing a second stage of core accretion. Such an epoch could have formed WASP-47 e and d, and perhaps additional yet-undetected low-mass planets dominated by rocky interiors.

In WASP-47, the orbital architecture and compositions of the low-mass planets provide insight into the formation history of the giant planets. Because the low-mass planets need a small amount of gas to form (about $0.3\,{M}_{\oplus }$), their formation must occur before the gas disk completely dissipates. This puts a deadline on when planet b can finish its migration; it must park at its current location before planets e and d form (otherwise it destroys them). In other words, planet b must arrive at its present location before the gas disk dissipates.

Because a gas disk damps planet eccentricities on the timescale of 10³ years (Papaloizou & Larwood 2000), high-eccentricity migration cannot proceed in a gas disk. Furthermore, the timescale for type I migration in a gas disk is much faster than the timescale for HEM,  so gas disk migration dominates while the gas disk is present. If we consider a mostly gas-depleted disk, we can try to migrate WASP-47 b inward via HEM in a short window of time before the disk completely dissipates. However, some force must circularize planet b's orbit quickly enough to allow the neighboring small planets to form. In the absence of gas, the only circularizing force is tides raised on the planet by the star. The timescale for eccentricity damping due to tides (Murray & Dermott 2000, Equation (4.198)) is

$$\begin{eqnarray}&&{\tau }_{e}=-\displaystyle \frac{e}{\dot{e}}=\displaystyle \frac{4}{63}\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{\star }}{\left(\displaystyle \frac{a}{{R}_{{\rm{p}}}}\right)}^{5}\displaystyle \frac{{\mu }_{\mathrm{pl}}{Q}_{\mathrm{pl}}}{2\pi /P},\end{eqnarray} \tag{ 7 }$$

where ${\mu }_{\mathrm{pl}}$ is the ratio of the elastic to gravitational forces in the planet $\approx {({10}^{4}\mathrm{km}/{R}_{{\rm{p}}})}^{2}$, and Q_pl is the tidal dissipation factor of the planet. For $a=0.05\,\mathrm{au}$ and ${Q}_{\mathrm{pl}}={10}^{6.5}$ (typical for hot Jupiters, Jackson et al. 2008), the circularization time is ${\tau }_{e}\approx {10}^{8}\,\mathrm{years}$, which is an order of magnitude longer than the disk lifetime. There is not enough time to move WASP-47 b to its present position via HEM and then circularize its orbit before the little planets form. Therefore, a history of HEM is highly unlikely for planet b.

On the other hand, planet b could also have formed in situ or nearly in situ in the gas-rich disk. Antonini et al. (2016) find through N-body simulations that, for the majority of observed Jupiter pairs (i.e., systems with one Jupiter-mass planet inside 1 au and one Jupiter-mass planet outside 1 au), the inner planet is unlikely to have formed beyond 1 au. Thus, it seems probable that WASP-47 b formed inside 1 au, and thus inside the snowline.

Since planet c has modest eccentricity ($0.28\pm 0.02$, it almost certainly had some sort of planet–planet interaction in the past. Although it is possible for instabilities in the gas disk to excite giant planet eccentricities, this effect only works for $e\lt 0.1$ and therefore cannot explain the eccentricity of planet c (Duffell & Chiang 2015). As described above, planet b is not a likely candidate for pumping the eccentricity of planet c. Thus, additional massive bodies that exchanged angular momentum with planet c in the past might still be present in the WASP-47 system. Alternatively, planet c might have ejected whatever companion pumped its eccentricity.

Based on the above arguments, evidence of two-stage planet formation in WASP-47 could include:

• 1.  
  a C/O ratio for WASP-47 b, e, and d consistent with formation inside the snowline;
• 2.  
  alignment of the stellar spin axis with the orbits of the transiting planets, consistent with a coplanar, mostly circular orbital history for planet b;
• 3.  
  an additional massive planet/sub-stellar/stellar companion with an eccentric orbit that exchanged angular momentum with planet c in the past.

A two-stage formation mechanism might also best explain the current architectures of other multi-planet systems with diverse planet compositions, such as 55 Cancri and Kepler-89 (KOI-94). An improved census of giant planets accompanying the low-mass planetetary systems discovered by Kepler will enable studies of the occurrence of planetary systems with diverse compositions like that of our solar system, illuminating the dominant physical processes in their formation.

Continued observations of this system and other systems with diverse planet compositions from the James Webb Space Telescope, the WFIRST mission, and other facilities will provide evidence for or against a two-stage planet formation scenario.