WASP-107b's Density Is Even Lower: A Case Study for the Physics of Planetary Gas Envelope Accretion and Orbital Migration

We collected a total of 60 radial-velocity measurements of WASP-107 over 56 nights between 2017 and 2020 using HIRES (Vogt et al. 1994) on the Keck I telescope. The observations were conducted by the California Planet Search (CPS; PI: Andrew Howard). The bulk of the observations were obtained as part of a K2-NASA key project between 2017 and early 2018. Following this program, we began remonitoring the system in 2019 to better constrain the orbit of a possible second planet. The data were acquired using the "C2" decker, integrating until the exposure meter reached 60,000 counts (signal-to-noise ratio (S/N) ∼ 100/pixel, ∼15 minutes). Wavelength calibrations were done using the iodine cell (Butler et al. 1996). The HIRES data reduction followed the standard procedures of the CPS described in Howard et al. (2010). Typical HIRES RV uncertainties are ∼1.5 m s⁻¹. For two of the nights, in which we obtained three measurements, we bin the data in units of 0.1 days, yielding a total of 57 HIRES RVs (Table 1).

We supplement our radial-velocity data set with 31 CORALIE measurements obtained between 2011 and 2014, as reported in Anderson et al. (2017). The CORALIE measurements have a typical uncertainty of ∼10 m s⁻¹ (see Table 1), 6 times larger than the typical uncertainties from HIRES, but improve our baseline. With the 31 CORALIE RVs, kindly provided to us by D. R. Anderson, our combined data set consists of 91 RVs spanning ∼3400 days, with a 1072 day gap between the two (Figure 1). These RVs are included in the machine-readable version of Table 1.

Zoom In Zoom Out Reset image size

We use isoclassify (Huber et al. 2017) to provide updated estimates for the stellar radius, mass, age, and luminosity. We adopt as inputs to isoclassify the T_eff and [Fe/H] values from the Keck/HIRES spectra, constrained by the SpecMatch-emp tool (Yee et al. 2017). SpecMatch-emp provides constraints on the the stellar parameters by comparing the Keck/HIRES spectra with an empirical spectral library. We use JHK 2MASS magnitudes, while the parallax, right-ascension, and decl. are taken from Gaia DR2 (see Table 2).

Table 2. Stellar Parameters of WASP-107

Parameter	Value	Source
Identifying Information
EPIC ID	228724232
αJ2000 R.A. (hh:mm:ss)	12:33:32.84	Gaia DR2 (Gaia Collaboration et al. 2018)
δJ2000 Decl. (dd:mm:ss)	−10:08:46.22	Gaia DR2 (Gaia Collaboration et al. 2018)
Photometric Properties
B (mag)	14.62 ± 1.20	Høg et al. (2000)
V (mag)	11.47 ± 0.20	Høg et al. (2000)
G (mag)	11.1740 ± 0.0009	Gaia DR2 (Gaia Collaboration et al. 2018)
J (mag)	9.378 ± 0.021	2MASS (Cutri et al. 2003)
H (mag)	8.777 ± 0.026	2MASS (Cutri et al. 2003)
K (mag)	8.637 ± 0.023	2MASS (Cutri et al. 2003)
Spectroscopic and Derived Properties
μα (mas yr−1)	−96.647 ± 0.110	Gaia DR2 (Gaia Collaboration et al. 2018)
μδ (mas yr−1)	−9.483 ± 0.058	Gaia DR2 (Gaia Collaboration et al. 2018)
Parallax (mas)	15.4175 ± 0.0617	Gaia DR2 (Gaia Collaboration et al. 2018)
Distance (pc)	64.7 ± 0.3	Bailer-Jones et al. (2018)
Spectral Type	K6	Anderson et al. (2017)
Teff (K)	4425 ± 70	This Paper (SpecMatch-emp)
[Fe/H]	+0.02 ± 0.09	This Paper (SpecMatch-emp)
$\mathrm{log}\,{g}_{\star }$	4.633±0.012	This Paper (isoclassify)
M⋆ (M⊙)	0.683${}_{-0.016}^{+0.017}$	This Paper (isoclassify)
R⋆ (R⊙)	0.67±0.02	This Paper (isoclassify)
L⋆ (L⊙)	0.132 ± 0.003	This Paper (isoclassify)
Prot (days)	17 ± 1	Anderson et al. (2017)
Age (Gyr)	8.3 ± 4.3	Močnik et al. (2017); Isochronal (bagemass)
	${6.9}_{-3.4}^{+3.7}$	This Paper; Isochronal (isoclassify)
	3.4 ± 0.7	This Paper; Gyrochronological

Download table as:  ASCIITypeset image

Our reported value of $\mathrm{log}\,g$ in Table 2 was obtained using the grid-modeling method in isoclassify with a loose input constraint on $\mathrm{log}g$ (4.7 ± 0.2). The constraint on $\mathrm{log}g$ that we obtain from this first fit is used as an input to a second fit using the direct method.

The parameters that were derived by both the direct and grid-modeling methods (stellar radius and luminosity) are consistent within less than 1σ. We report the radius and luminosity determined using the direct method, because in this mode the radius is not affected by the uncertainty on stellar models. The values that we quote for the stellar mass and age are the results obtained using the grid-modeling method.

The values of all parameters yielded by isoclassify are consistent with Anderson et al. (2017) and Močnik et al. (2017), and our mass estimate is more precise, which we attribute to the higher precision of the T_eff and [Fe/H] obtained from our HIRES spectra.

The updated stellar parameters are listed in Table 2. We use the new stellar radius in order to update the radius and semimajor axis of WASP-107b; taking R_b /R_⋆ and a_b /R_⋆ from Dai & Winn (2017), we find R_b  = 0.96 ± 0.03R_J and a_b  = 0.0566 ± 0.0017 au.

The age of WASP-107 was previously estimated using several methods: an isochronal age estimate using bagemass (Maxted et al. 2015) yielded 8.3 ± 4.3 Gyr—in agreement with the constraint we obtain with isoclassify, and a conflicting gyrochronological age of 0.6 ± 0.2 Gyr (Močnik et al. 2017). This discrepancy can be explained by the failure of standard spin-down models to reproduce observations of reduced braking efficiency for late-K dwarfs relative to their F- and G-type counterparts, leading them to experience a period of stalled rotation periods past ∼700 Myr (Agüeros et al. 2018; Curtis et al. 2019, 2020, 2020).

To resolve the discrepancy, we use a new gyrochronology model that accounts for stalled magnetic braking (R. Angus et al. 2020, in prep) in order to estimate the age of WASP-107. This model was calibrated by fitting a Gaussian process (GP)—a semiparametric model that is flexible enough to capture the complex nature of stellar spin-down, to a number of asteroseismic stars and open clusters, including NGC 6811, which exhibits stalled magnetic braking (Curtis et al. 2019). We also use kinematic ages of Kepler field stars to calibrate the model for old K and early M dwarfs, where there is a dearth of suitable open cluster calibration stars. These kinematic ages also reflect the stalled magnetic braking behavior seen in open clusters (Angus et al. 2020).

Using our new gyrochronology model, we infer an age of 3.4 ± 0.3 Gyr for this star using the measured rotation period of 17 ± 1 days (Anderson et al. 2017). The 0.3 Gyr uncertainty is the formal uncertainty that results from the uncertainty on the star's rotation period, and does not account for uncertainty in the model. Quantifying the magnitude of the model uncertainty is beyond the scope of this paper, and we adopt a 20% value of 0.7 Gyr as a more reasonable estimate of the true age uncertainty (Table 2). Importantly, we resolve the conflict between the age estimates and conclude that WASP-107 is not as young as previously suggested by gyrochronology models.

We fit the RV measurements using RadVel, an open source Python package for fitting Keplerian orbits to radial-velocity data sets (Fulton et al. 2018). The posterior distributions of the parameters are sampled using the Markov Chain Monte Carlo (MCMC) affine-invariant sampler emcee (Foreman-Mackey et al. 2013). We fix the period P_b and time of conjunction T_conj,b of WASP-107b to the K2 transit ephemeris reported in Dai & Winn (2017). We use the conventional unbiased basis formed by the jump parameters $h=\sqrt{e}\cos \omega$ and $k=\sqrt{e}\sin \omega$ to fit for orbital eccentricities and arguments of periastron. The MCMC fitting basis is thus (K_b , $\sqrt{{e}_{b}}\cos {\omega }_{b}$, $\sqrt{{e}_{b}}\sin {\omega }_{b}$), along with (P_c , T_conj,c , K_c , $\sqrt{{e}_{c}}\cos {\omega }_{c}$, $\sqrt{{e}_{c}}\sin {\omega }_{c}$) for two-planet models. We introduce an offset term γ_i and a jitter term σ_i for each instrument (i = {H, C} for HIRES and CORALIE). These terms account respectively for different zero-points between HIRES and CORALIE and for additional noise (e.g., stellar jitter) not encapsulated by the single-measurement errors. The offset terms are computed using a linearized analytic solution. We invoke physical priors on the RV semiamplitude K > 0 and eccentricity e ∈ [0, 1). This brings the total of fitting parameters to 17 for a joint two-planet fit of the HIRES and CORALIE data.

We perform a maximum a posteriori (MAP) fit and use the best-fit parameters to seed an MCMC that estimates the full posterior. We run the MCMC with 4 ensembles of 200 walkers, each of which is allowed to burn in until the Gelman–Rubin statistic is <1.01, checking for convergence after 2000 steps per walker. We then save every 5th step as a posterior sample until either 10,000 samples per walker are reached or the Gelman–Rubin statistic decreases below 1.005, at which point we consider the chains to be well mixed. A second MAP fit is then run, starting at the median values determined by the MCMC posteriors.

We also attempted to use GP regression in order to mitigate the effect of stellar activity on the RVs by training the GP on ancillary K2 time series, but found that adding the GP parameters to the RV fit was disfavored from a model comparison point of view (using Bayes factors; see Section 4.2) and did not improve our constraints on the planet parameters.

We perform a total of four fits: a one-planet and two-planet fit to the HIRES RVs alone, and the same two fits to the combined HIRES and CORALIE data set. Model comparison is performed within a robust Bayesian framework, using the Bayesian evidence and Bayes factors assuming equally probable hypotheses a priori. In that framework, a higher Bayesian evidence value for a model that includes more free parameters indicates that the additional model complexity is granted by the quality of the data. For each model of the RVs, we approximate the Bayesian evidence using the estimator presented in Perrakis et al. (2014) (as in, e.g., Cloutier et al. 2019). This method consists of using the product of marginalized posterior distributions of "blocks" of parameters (obtained after a rearrangement of the parameters in MCMC chains) as an importance sampler in the evidence estimation. Nelson et al. (2020) found that this estimator obtains results in good agreement with more computationally intensive methods (e.g., nested sampling) while requiring no additional sampling of the parameter space. On the other hand, they found that non-Bayesian methods, such as Bayesian information criteria (BIC), that are based on simplifying assumptions such as an infinitely narrow posterior only poorly approximate the Bayesian evidence and can lead to false negatives in planet detections using RVs. In this work, we use the implementation of this estimator in the bayev package, publicly available on GitHub. ¹⁷ Each block consists of only one of the fitted parameters. The importance sampler is obtained after a random rearrangement of the parameters, and we use kernel density estimations of each marginalized posterior distribution of samples. For each model, we perform 1000 such rearrangements, calculating the Bayesian evidence for each, and quote the median from these values.

The HIRES data exhibit a significant long-period trend on top of the signal from WASP-107b because of the presence of a second planet in a long-period orbit (Figure 1). Our Bayesian model comparison based on the HIRES data alone favors the two-planet model at 4.9 × 10⁷ to 1 in the Bayes factor (equivalently 6.3σ; Table 3). The evidence for the existence of this second planet, WASP-107c, is further strengthened by including the previously obtained CORALIE data set in the radial-velocity analysis (see Table 3). Overall, the CORALIE data set increases the baseline by a factor of ∼3, allowing us to rule out a 1-planet model in favor of a two-planet Keplerian solution, with a Bayes factor of 1.2 × 10⁵⁸ to 1 (or a 16.5σ significance). The best-fit two-planet model to the combined data set is shown in Figure 1. In what follows, we will report constraints on the orbital parameters based on the two-planet model and the combined data set.

Table 3. Keplerian Model Comparison

Model	$\mathrm{ln}p(D| {{ \mathcal M }}_{i})$	$\tfrac{p(D| {{ \mathcal M }}_{2})}{p(D| {{ \mathcal M }}_{1})}$	"Sigma" a
HIRES
1-Planet	−174.3	⋯	⋯
2-Planet	−156.6	4.9 × 107	6.3σ
HIRES + CORALIE
1-Planet	−290.1	⋯	⋯
2-Planet	−156.4	1.2 × 1058	16.5σ

Notes. The Bayes factors are reported relative to the single-planet model. In both cases, the two-planet model is favored.

^a The correspondence between the Bayes factors and the "sigma" significance in a frequentist framework is calculated using p-values (Trotta 2008).

Download table as:  ASCIITypeset image

MCMC parameter estimates and uncertainties are displayed in Table 4 along with derived planet parameters. The corresponding joint and marginalized posterior distributions are shown in Figure 2. We find a revised mass for WASP-107b of M_b  = 30.5 ± 1.7 M_⊕ and an eccentricity of e_b  = 0.06 ± 0.04. This lower mass associated with a Jupiter radius makes WASP-107b an outlier in the mass–radius diagram (Figure 3). Combined with the transit measurement for the radius, we find the bulk density of WASP-107b to be ${\rho }_{b}={0.134}_{-0.013}^{+0.015}$ g cm⁻³, even lower than the previous estimate of 0.19 ± 0.03 g cm⁻³ (Anderson et al. 2017).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. MCMC Posteriors for the WASP-107 System

Parameter	Unit	Credible Interval
Orbital Parameters
Pb	days	≡5.7214742
Tconj,b	BJDTDB	≡2457584.329897
eb		0.06 ± 0.04
ωb	deg	${40}_{-60}^{+40}$
Kb	m s−1	14.1 ± 0.8
Pc	days	${1088}_{-16}^{+15}$
	yr	2.98 ± 0.04
Tconj,c	BJDTDB	${2458520}_{-70}^{+60}$
ec		0.28 ± 0.07
ωc	deg	$-{120}_{-20}^{+30}$
Kc	m s−1	${9.6}_{-1.0}^{+1.1}$
Derived Parameters
Mb	M⊕	30.5 ± 1.7
	MJ	0.096 ± 0.005
ρb	g cm−3	${0.134}_{-0.013}^{+0.015}$
${M}_{c}\,\sin i$	M⊕	115 ± 13
	MJ	0.36 ± 0.04
${S}_{c}^{{\prime} }$ a	mas	${26}_{-5}^{+8}$
Global Parameters
γH b	m s−1	≡1.52001
γC b	m s−1	≡3.55333
σH	m s−1	${3.9}_{-0.4}^{+0.5}$
σC	m s−1	${5.5}_{-2.9}^{+2.8}$

Notes.

^a Sky-projected angular separation. ^b Reference epoch for γ: 2457934.793533.

Download table as:  ASCIITypeset image

The HIRES and CORALIE data sets combined constrain the RV semiamplitude of the outer companion to ${K}_{c}={9.6}_{-1.0}^{+1.1}$ m s⁻¹, corresponding to a mass of ${M}_{c}\,\sin \,i=0.36\pm 0.04$ M_J. The orbital period of the outer planet (2.98 ± 0.04 yr) and its eccentricity (e_c  = 0.28 ± 0.07) are well constrained thanks to the detection of two steep rises in RV throughout the HIRES observing campaign. Continued monitoring of this system will enable further refinements of the orbit of this outer planetary companion.

Our measurement of WASP-107b's mass of 30.5 ± 1.7 M_⊕ is substantially lower than the 38 ± 3 M_⊕ reported by Anderson et al. (2017) from the analysis of the CORALIE data set. We argue that this discrepancy can be explained as arising from the improved precision of the HIRES RV measurements.

Our new mass estimate for WASP-107b is 20% lower than the previously published one, and lies  ≈2.2σ away from it. For a fitted semiamplitude K to a time series of N_obs radial velocities having a typical measurement uncertainty σ, K will be biased toward a higher value than the true one if K/σ or N_obs is small (Shen & Turner 2008). Thus, with N_obs = 31 and K/σ < 2, the CORALIE data set is expected to suffer from this bias more than the HIRES measurements, where N_obs = 60 and K/σ ≈11. Moreover, Anderson et al. (2017) fixed the eccentricity to zero (yielding smaller uncertainties on K) and used a larger stellar mass compared to our new determination, which further contributed to a larger estimated planet mass.

Our new mass measurement indicates that WASP-107b has an extraordinarily low density. We infer the gas-to-core mass ratio of WASP-107b by comparing our measurements to model predictions for various envelope mass fractions (Figure 3). We generate grids of thermal evolution models and predict 10 mbar radii for planets with varying incident fluxes S_inc, masses M_p , ages and fractions f_env[%] of their masses contained in the H/He envelopes surrounding their cores. Previously available model grids in the literature (e.g., Lopez & Fortney 2014) do not extend all the way to the mass and radius of WASP-107b and could not yet make use of the latest equation of state for hydrogen.

We compute the thermal evolution models as described in Thorngren et al. (2016). The atmosphere models are interpolated from the Fortney et al. (2007) solar-metallicity grids and serve as a boundary condition to interior structure models. We use updated equations of state (EOS) for the solar-ratio H/He gas (Chabrier et al. 2019) and choose the Thompson (1990) EOS for metals. Because the composition of the core of WASP-107b is unknown, we marginalize over the uncertainty on the fraction of water within its core by producing two grids for planets with either an Earth-like core (67.5% rock, 32.5% Fe by mass) or a pure water composition. Therefore, the true envelope mass fraction of the planet should lie in the range cornered by the two derived envelope mass fractions: the same planet radius and mass will yield a larger inferred H/He mass fraction if the core composition is rocky than if it is water-rich. The two-layer models include no mixing from the core into the envelope and in the pure-water core case, we include convection in the water component. Our grids span the range 10 to 10⁵ S_⊕ for incident flux, 15 to 40 M_⊕ for planet mass, 0.1 to 14 Gyr in terms of system age, and H/He envelope mass fractions of 0.01 to 100%. The curves of constant f_env in the mass–radius diagram computed using the median parameters for WASP-107b already suggest an extremely high envelope mass (Figure 3).

More quantitatively, we solve the inverse problem of inferring an envelope mass fraction from the planet's incident flux, mass, age, and radius using the smint (Structure Model INTerpolator) interpolation and envelope mass fraction fitting package, which we made publicly available on GitHub. ¹⁸ From any set of parameters (S_inc, M_p , age, f_env), smint returns a planet radius by performing linear interpolation over a grid of f_env, ${\mathrm{log}}_{10}{M}_{p}/{M}_{\oplus }$, system age, and ${\mathrm{log}}_{10}{S}_{{inc}}/{S}_{\oplus }$ setup using one of our model grids. We run an MCMC that fits for the combination of (S_inc, M_p , age, f_env), which best matches the observed planet radius. We adopt Gaussian priors on S_inc, M_p , and the system's age informed by the parameters derived from the RV fits and the gyrochronology model. We use a uniform prior on the planet's envelope mass fraction over the entire range spanned by the grid. Each of the 100 chains is run for 10,000 steps, 60% of which are discarded as burn-in. The autocorrelation times of the fitted parameters are all ≤52 steps (less than 15% of the length of the chains past the burn-in phase), securing that our chains are converged and sample well the posterior distributions. We perform two fits using the grids for the Earth-like and pure H₂O core compositions.

We infer upper limits on the core mass of 3.7 and 4.6 M_⊕ at 3σ (lower masses are excluded with 99.7% confidence) for the Earth-like and pure H₂O core respectively, corresponding to envelope mass fractions of 88% and 85% for a solar-metallicity envelope (Figure 4). As expected, we find a correlation between the system's age and the inferred envelope mass fraction due to the contraction and cooling of the planet as it ages. The inferred core mass is much smaller than traditionally assumed to be required for the triggering of runaway gas accretion. Interestingly, the small core still must have been able to accrete almost 30 M_⊕ in gas from the surrounding nebula once it had coagulated.

Zoom In Zoom Out Reset image size

Our estimate also serves as a conservative lower limit on the initial accreted envelope mass, given that planets undergo atmospheric loss throughout their existence (e.g., Owen & Wu 2013, 2016). In fact, the puffiness and close-in orbital distance of WASP-107b places it at the cusp of the sub-Jovian desert (e.g., Mazeh et al. 2016; Owen & Lai 2018). In particular, we find that assuming a zero Bond albedo, WASP-107b's core would need to be more massive than 4.1 M_⊕ (1σ upper limit accounting for the uncertainties on the stellar and orbital parameters) for the atmosphere to survive photoevaporation at its present orbital location (Ginzburg et al. 2016). This limit is higher than our abovementioned estimates of the core mass, suggesting that WASP-107b migrated to its present orbital location only recently (see Section 6.5).

The large H/He mass fraction of WASP-107b strikingly contrasts with those of the ice giants of our solar system (5%–15% H/He for both Uranus and Neptune, Podolak & Marley 1991; Guillot 2005). Instead, WASP-107b's structure most closely resembles those of Jupiter and Saturn, composed of more than 90% H/He by mass. WASP-107b is an extreme case even among the vastly diverse exoplanet interiors illustrated by the wide span in radius of the super-Neptunes, which highlights the great diversity of formation pathways within this population. WASP-107b, with its Jupiter-like composition, resembles K2-161b and Kelt-11b, which have large envelope mass fractions (>75%) associated with their low densities (Pepper et al. 2017; Brahm et al. 2019). Meanwhile, thermal evolution models predict for instance a Neptune-like structure for the 23 M_⊕ planet HAT-P-11b, with about 15% H/He by mass (Petigura et al. 2017).

According to the core-nucleated accretion theory (Perri & Cameron 1974; Mizuno 1980; Pollack et al. 1996), which stands as the "standard model" for planet formation, a planetary core will first assemble from the solids present in its neighborhood within the protoplanetary disk. Past this phase of solid accretion, a second phase will start, dominated by gas accretion. The rate of envelope accretion is governed by cooling of the gas (Lee & Chiang 2015; Ginzburg et al. 2016) and by local and global hydrodynamic flows (Tanigawa & Tanaka 2016). The observed distribution of GCRs can effectively be reproduced if all theses mechanisms are taken into account (Lee 2019). For core masses below ≈20 M_⊕, the cooling timescale remains the longest until the runaway is triggered (Figure 2 of Lee 2019), beyond which the rapid dispersal of the background disk limits the delivery of gas onto planets. We note here that this rapid dispersal is valid only within ∼1 au, where the photoevaporative wind is launched, carving out a gap disconnecting the inner disk from the outer disk (Owen et al. 2011).

To reproduce the inferred GCR of WASP-107b, we need to identify the gas accretion scenario in which the runaway growth begins but is curbed prematurely within the typical disk lifetime ∼1–10 Myr (e.g., Haisch et al. 2001; Mamajek 2009; Alexander et al. 2014; Pfalzner et al. 2014). We explore two potential pathways, motivated by the results of Lee & Chiang (2015). First, we consider dust-free accretion at distances beyond the present-day orbit of WASP-107b. Second, we consider in situ (0.05 au) accretion of dusty, metal-enriched gas.

As dust grains coagulate and sediment, they no longer contribute to the atmospheric opacity (e.g., Mordasini 2014; Ormel 2014; Piso et al. 2015). The upper radiative layer becomes largely isothermal such that the envelope becomes colder at larger stellocentric distances. As the internal modes of gaseous molecules gradually freeze out, the opacity drops and the envelope cools faster and therefore grows faster. Dust-free accretion beyond ∼1 au was invoked as a way to explain Kepler-51b, a super-puff of a ∼2 M_⊕ core that accreted a ∼30% by mass envelope (Lee & Chiang 2016).

For WASP-107b, we adopt the procedure of Lee (2019)—where at each timestep, the rate of gas accretion is set by the minimum between cooling, local hydrodynamic, and global disk accretion—but replace the cooling calculation to that of dust-free accretion. One-dimensional structure equations are solved numerically following Lee et al. (2014), corrected for the decrease in the bound radius by three-dimensional advective flows (e.g., Lambrechts & Lega 2017). We confirm that the updated numerical solutions computed at a few orbital distances agree well with the analytic scaling relationship derived by Lee & Chiang (2015), and search for the range of potential formation locations by scaling with respect to the fiducial numerical calculation at 2 au (see top row of Figure 5). For a core that weighs 2 M_⊕, dust-free accretion should have occurred at ≳1 au for the runaway to be launched.

Zoom In Zoom Out Reset image size

Are there ways to form WASP-107b in situ? It appears impossible for solar-metallicity gas, but the rate of accretion can be boosted in highly metal-enriched environments. While this may be counterintuitive, because higher opacity would slow down the cooling of gas, this impediment by opacity endures only up to Z ∼ 0.2 (10 times the solar value). Beyond this point, the mean molecular weight of the envelope is so high that it becomes more prone to gravitational collapse, requiring more vigorous accretion to maintain hydrostatic equilibrium. This boost in accretion rate applies equally to dust-free and dusty accretion. Because dusty accretion is generally slower than dust-free accretion, we elect to show the former to delineate all the possible limiting cases.

We proceed the same way as in the calculation for dust-free accretion and follow the procedure of Lee (2019), accounting for the effect of metallicity in the rate of cooling, which can be written as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{M}_{\mathrm{atm}}}{{M}_{\mathrm{core}}} & \sim & 0.09\,{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{gas}}}{13{\rm{g}}{\mathrm{cm}}^{-3}}\right)}^{0.12}{\left(\displaystyle \frac{{M}_{\mathrm{core}}}{20{M}_{\oplus }}\right)}^{1.7}{\left(\displaystyle \frac{t}{0.1\mathrm{Myr}}\right)}^{0.4}\\ & & \times \,\left(\displaystyle \frac{0.02}{Z}\right){\left(\displaystyle \frac{\mu }{2.37}\right)}^{3.4}\mathrm{Exp}\left(\displaystyle \frac{t}{{t}_{\mathrm{run}}}\right)\end{array}\end{eqnarray} \tag{ 1 }$$

where Σ_gas is the disk gas surface density evaluated at the position of the core and at time t, ¹⁹ μ is the mean molecular weight, and t_run is the runaway timescale determined numerically as the time at which ${M}_{\mathrm{atm}}/{M}_{\mathrm{core}}$ reaches 0.5:

$$\begin{eqnarray}&&{t}_{\mathrm{run}}\sim 2.24\,\mathrm{Myr}\,{\left(\displaystyle \frac{20{M}_{\oplus }}{{M}_{\mathrm{core}}}\right)}^{4.2}\left(\displaystyle \frac{Z}{0.02}\right){\left(\displaystyle \frac{2.37}{\mu }\right)}^{8.5}.\end{eqnarray} \tag{ 2 }$$

Both Z and μ are evaluated at the radiative-convective boundary, assuming Z is uniform throughout the planetary envelope.

The required metallicity enhancement is severe, because of the need to trigger the runaway prior to disk dissipation and reach the observed GCR (see bottom row of Figure 5). We find that the required metallicity is Z ≳ 0.95, or approximately 48 times the solar value by mass.

We use the same type of interior structure models as in Section 5.4 in order to constrain the range of metallicities allowed by the mass and radius of WASP-107b. This time, the planet is modeled as a well-mixed, fully convective envelope with a radiative atmosphere, and we assume a 50–50 mixture of rock and ice for metals. We do not include a core to get a conservative (high) upper limit on the envelope metallicity (see Thorngren & Fortney 2019). Our models are sensitive to the assumed age of the system, because the planet cools and contracts over time. We set the prior on the system age to a normal distribution corresponding to the gyrochronological age (3.4 ± 0.7 Gyr; see Table 2).

We fit for the metal mass fraction using MCMC with 6 chains in parallel. In each chain, we burn in for 10,000 steps and then record every 10th step until we have amassed 100,000 samples. We infer a 3σ upper limit of Z ≲ 0.25 on the envelope's metal mass because higher metallicities would result in smaller radii. Interestingly, our constraint on the metallicity (<12.5 times the solar value by mass at 3σ) is in agreement with the metallicity constraint retrieved from the spectroscopic characterization of the 12.6 M_⊕ sub-Neptune GJ 3470b (Benneke et al. 2019) and suggests that the two planets may have had similar formation processes.

The 3σ upper limit of 0.25 on the envelope metallicity is 3.8 times smaller than required by the dusty gas accretion models and effectively rules out the high-metallicity accretion scenario: the planet could not be this puffy to this day if such large amounts of metals were present in its atmosphere. We highlight the fact that the rate of dusty accretion is only weakly dependent on orbital distance, under the assumption that the dust-to-gas mass ratio is spatially constant within the disk. Our constraints on the envelope metallicity therefore rule out the dusty accretion scenario at any orbital distance in disks with uniform dust-to-gas mass ratio. Furthermore, as mentioned above (Section 5.4), the observationally inferred GCR provides only a lower limit on the amount of gas accreted by the core, because mass loss powered by stellar irradiation and the planet's leftover heat from formation can pare down the envelope mass fraction and is believed to sculpt the radius distribution of small planets (Owen & Wu 2017; Ginzburg et al. 2018; Gupta & Schlichting 2020). Reaching a higher GCR at formation would necessitate an even higher metallicity.

This leaves dust-free accretion as the only plausible formation scenario for WASP-107b's massive envelope. With over 85% of its mass contained in its H/He envelope, WASP-107b effectively provides evidence for the possibility of a <4.6 M_⊕ core accreting extensive amounts of primordial gas.

Assuming that WASP-107b indeed formed as a dust-free world at ≳1 au, it must have undergone inward migration, whether through the underlying disk in the early, gas-rich environment or through perturbation by neighboring planets in late, gas-poor or gas-empty environment. We propose that the migration of WASP-107b happened early, such that it went close to experiencing the full extent of runaway gas accretion before it was deprived from its supply of gas.

The fact that WASP-107b accreted copious amount of envelope suggests that it assembled early, beyond ∼1 au. After its assembly, the planet is expected to undergo inward Type I migration (Lin et al. 1996; Papaloizou et al. 2007)—this has been proposed as the formation channel of super-puffs explaining why they are frequently observed as parts of resonant chains (Lee & Chiang 2016). Alternatively, or additionally, WASP-107b could have arrived at its present orbit through planet–planet scattering involving WASP-107c: this scenario would have the advantage of explaining the eccentricity of planet c. Still, we note that it is uncertain whether the runaway accretion can be stopped at the right time to produce WASP-107b beyond 1 au, which may suggest that WASP-107b began its inward transport just at the right time as it reached its envelope mass fraction of >85%.

Although our constraints on the eccentricity of WASP-107b are consistent with a circular orbit at 2σ, we find that it could have a moderate eccentricity. The timescale of tidal circularization (${\tau }_{\mathrm{circ}}=e/\dot{e}$) for WASP-107b is ∼66 Myr, much shorter than the estimated age of the system. ²⁰ Any remaining eccentricity of the inner planet resulting from dynamical interactions will therefore be quickly dampened by tidal interactions with its host star. Given this short circularization timescale, a present-day eccentricity of 0.06 (see Table 4) would imply recent dynamical interactions that excited the eccentricity to high values. In particular, we demonstrate using our constraints on the WASP-107b's orbit that tidal interactions alone could have brought the planet from ≳0.10 au to its current orbital location over less than 1 Gyr (Figure 6).

Using the occurrence of spot-crossing anomalies during planetary transits, WASP-107b was inferred to feature a high spin–orbit misalignment (40–140 degrees; Dai & Winn 2017; Močnik et al. 2017; Rubenzahl et al. 2021). Such high obliquity is often considered a sign of dynamic upheaval and could be attributed to nodal precession, as was proposed for HAT-P-11b (Yee et al. 2018). Alternatively, WASP-107b is sufficiently small in mass and sufficiently far away from the star that its natal disk could have been driven to substantial spin–orbit misalignment with respect to the central star. For a star-disk system with an outer planet that is misaligned with the disk, secular resonance between the disk's precession around this planet and the star's precession around the disk can drive a large star-disk misalignment (e.g., Batygin & Adams 2013; Lai 2014; Zanazzi & Lai 2018). This resonance is suppressed for rapidly rotating stars in the presence of a close-in massive inner planet (i.e., the star's precession around the inner planet is faster than the disk's precession around the outer planet). We confirm that for the WASP-107 system, star-disk misalignment by secular resonance is possible as long as WASP-107b is beyond ∼0.02 au, which it is (Equation (70) of Zanazzi & Lai 2018). Another scenario that could account not only for the spin–orbit misalignment but also for WASP-107b's moderate eccentricity is that of a resonance initiated during the dispersal of an otherwise coplanar protoplanetary disk (Petrovich et al. 2020). In this case, one would expect a large mutual inclination between planets b and c.

Alternatively, we examined whether Kozai–Lidov cycles induced by the presence of planet c could explain the high obliquity of WASP-107b and its moderate eccentricity. We calculate the Kozai timescale (Kiseleva et al. 1998) using our posterior chains and considering a random inclination for planet c. Comparing it to the period of general relativistic (GR) apsidal precession (see Figure 7), we find that the Kozai cycles are likely suppressed by GR precession, with τ_K  ≫ τ_GR.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The abovementioned scenario relies on the assumption that our two-layer models represent accurate end cases for the planetary structure. In this subsection, we consider the possibility that WASP-107b is instead a sub-Neptune with an inflated radius.

The amount of tides raised on a planet and their dissipation can increase dramatically for planets on eccentric, inclined orbits. In particular, it was suggested that WASP-107b could be a puffed-up sub-Neptune with an envelope mass fraction as low as ∼10% assuming an obliquity of ∼60°, a tidal quality factor of the planet ${Q}_{p}^{{\prime} }\sim {10}^{5}$, and an eccentricity of ∼0.13 (approximately 2σ higher than the limit we derive from our RV analysis; Millholland et al. 2020).

Furthermore, while we considered above the possibility that the eccentricity of WASP-107b's orbit is associated with recent dynamical interactions, it can also be explained by a large tidal quality factor for the planet. We estimate that ${Q}_{p}^{{\prime} }\gtrsim 5\times {10}^{6}$ is required to maintain non-zero eccentricity over the system's age against tidal circularization. Using this ${Q}_{p}^{{\prime} }$ and the updated eccentricity ∼0.06, we expect significantly less inflation by tides with a dissipation rate of ∼7 × 10²⁵ ergs s⁻¹ (approximately one order of magnitude less than the estimates of Millholland et al. 2020). Without running an evolutionary model, we cannot determine the expected envelope mass fraction in the presence of tidal dissipation, but we estimate that it would be greater than ∼30%–40% based on Figure 7 of Millholland et al. (2020), which is large enough to have initial radii of  ≳4R_⊕ (i.e., a sub-Saturn). In this case, WASP-107b would have a core that weighs 18–21M_⊕. We find that in this case, its formation can be reconciled with in situ dusty accretion with Z ∼  0.1–0.3.

While the order of magnitude of the tidal quality factor is unconstrained observationally for this planet and higher values would lead to a lesser impact of tides on the planetary radius, tidal inflation could provide an alternative explanation for the low measured density of WASP-107b.