# Primordial Radius Gap and Potentially Broad Core Mass Distributions

of Super-Earths and Sub-Neptunes

###### Abstract

The observed radii distribution of Kepler exoplanets reveal two distinct populations: those that are more likely to be bare rocks () and those that are more likely to be gas-enveloped (). There exists a clear gap in the distribution of radii that separates these two kinds of planets. Mass loss processes like photoevaporation by high energy photons from the host star have been proposed as natural mechanisms to carve out this radius valley. These models favor underlying core mass function of sub-Neptunes that is sharply peaked at 6–8 but the radial-velocity follow-up of these small planets hint at a more bottom-heavy mass function. By taking into account the initial gas accretion in gas-poor (but not gas-empty) nebula, we demonstrate that the observed radius valley can be reconciled with core mass functions that are broad extending well into sub-Earth regime. The maximally cooled isothermal limit prohibits cores lighter than 1–2 from accreting enough mass to appear gas-enveloped. The rocky-to-enveloped transition established at formation produces a gap in the radius distribution that shifts to smaller radii farther from the star, similar to that observed. For the best agreement with the data, our late-time gas accretion model followed by photoevaporative mass loss favors dust-free accretion in hotter disks with a core mass function that is as broad as .

0000-0002-1228-9820]Eve J. Lee

## 1 Introduction

In galactic and stellar astronomy, the initial mass function of stars is one of the most fundamental quantity that influences the structural and chemical evolution of the interstellar medium and the galaxy on average. Obtaining an analogous mass function for exoplanets is challenging. Sub-Neptunes and super-Earths dominate the population with many of them at orbital periods beyond 10 days (e.g., Fressin et al., 2013; Petigura et al., 2013; Burke et al., 2015), where we lose sensitivity to measure their masses with e.g., radial velocity surveys (e.g., Weiss & Marcy, 2014). Mass measurements using transit timing variations are available for only a handful of planets in multi-planetary systems, being favorable to those near mean motion resonances (e.g., Wu & Lithwick, 2013; Hadden & Lithwick, 2014).

Theoretically, Malhotra (2015) derived a log-normal distribution of total mass (i.e., core + envelope mass) function peaked at 4–10 using the observed period ratio distribution and applying the condition for dynamical stability given by Hill spacing. Wu (2019) searched for a log-normal distribution of core masses that best-fits photoevaporation model to the observed distribution of planetary radii. They argued that a mass distribution sharply peaked at 8 was necessary to reproduce the shape of the “radius valley”, a gap in the radius distribution at 1.3–1.6 predicted by mass loss theory (Owen & Wu, 2013) and later confirmed by the California-Kepler Survey (Fulton et al., 2017; Fulton & Petigura, 2018) and asteroseismology (Van Eylen et al., 2018). Rogers & Owen (2020) performed a more sophisticated hierarchical inference analysis fitting photoevaporation model to the observed radius-period distribution and concluded a similarly peaked mass distribution (with mean at 6) is required.

Such high masses are at odds with the radial velocity follow-up of Kepler planets which reports peak masses as low as 1 (Weiss & Marcy, 2014). Furthermore, the true radius/mass distribution may be more bottom-heavy than previously thought (Hsu et al., 2019).

In this paper, we assess whether a power-law core mass distribution that extends to the sub-Earth masses is consistent with the observed radius distribution as well as the shape of the gap in the radius-period space. Instead of assuming a distribution of initial envelope mass fraction that is independent of core mass, we calculate the expected envelope mass from nebular accretion in the late stages of disk evolution, a gas-poor environment deemed favorable for preventing runaway gas accretion to ensure the formation of super-Earths and sub-Neptunes (Lee et al., 2014; Lee & Chiang, 2016).

## 2 Methods

### 2.1 Underlying core mass distribution

We begin with the ansatz that the underlying sub-Neptune/super-Earth core mass distribution follows a power-law distribution:

(1) |

where is the mass of the core and we choose ; is the best-fit power-law slope to the distribution of peak posterior masses of sub-Neptunes from the radial-velocity follow-up by Marcy et al. (2014). We note that in logarithm of , is top-heavy, is neutral, and is bottom-heavy. We also experimented with exponential distribution in linear and logarithm of and found them to provide poor match to the data. The minimum and the maximum core masses are set to 0.01 and 30.

### 2.2 Initial envelope mass fraction

For each core, its initial envelope mass fraction is calculated using the analytic scaling relationship derived by Lee & Chiang (2015) appropriate for gas accretion by cooling (equivalent to Phase II of the core accretion theory, Pollack et al. 1996; see also Ginzburg et al. 2016). We modify the expressions for the weak dependence on the nebular density (Lee & Chiang, 2016) and for the expected decrease in the bound radius due to three-dimensional hydrodynamic effects (Lambrechts & Lega, 2017; Fung et al., 2019). Shrinking the outer bound radius decreases the rate of accretion in a linear fashion (Lee et al. 2014; see also Ali-Dib et al. 2020 for understanding this effect in terms of entropy delivery). We verify that the expressions we provide here match the numerical calculations.

First, cores need to be sufficiently massive to accrete gas. We calculate the envelope mass only for cores that satisfy

(2) |

where (Valencia et al., 2006), is the outer radius of the bound envelope, is a numerical factor that takes into account the effect of three-dimensional advective flows, is the Hill radius, is the Bondi radius, is the disk temperature, is the orbital distance, and is a numerical coefficient and a free parameter. We note that for these small cores, inside 1 AU.

For dusty accretion, the envelope mass fraction

(3) |

where is the mass of the gaseous envelope, is the accretion time,

(4) |

We express equation 4 with the disk temperature . More precisely, the relevant temperature is that at the envelope radiative-convective boundary. The outer layers of dust-free envelopes are nearly isothermal so adopting obtains the same answer. Although equations 3 and 4 are derived assuming , adjusting for makes no significant difference.

Throughout this paper, (solar metallicity), , and is drawn from a logarithmically uniform distribution that range 0.01 and 1 Myr, consistent with the late-time formation scenario (Lee & Chiang, 2016). Motivated by Figure 11 of Fung et al. (2019), we explore and 0.2. We choose throughout, prompted by the required level of gas depletion to reproduce the observed peaks in period ratios just outside of first order mean-motion resonances (Choksi & Chiang, 2020).

For a given core mass, the maximum possible envelope mass that can be accreted is given by a fully isothermal profile (e.g., Lee & Chiang, 2015). No cores are allowed to accrete more than this maximally cooled isothermal mass:

(5) |

where is the local nebular volumetric density, is the Keplerian orbital frequency, is the local disk sound speed, is the Boltzmann constant, and is the mass of the hydrogen atom. The nebular mean molecular weight is assumed to be the same as that of the envelope.

### 2.3 Estimating radii

While the masses of sub-Neptunes are dominated by the cores, their radii are largely determined by their envelope mass fraction (Lopez & Fortney, 2014). We follow closely the procedure devised by Owen & Wu (2017) in converting envelope mass fractions to radii. Only the essential elements are shown here.

First, we assume that after the disk gas is completely dissipated and planets are laid bare to stellar insolation, their outer layers become isothermal and volumetrically thin (6 scale height above the radiative-convective boundary; Lopez & Fortney 2014). From the density profile given by the inner adiabat

(6) |

the total envelope mass

(7) |

where is the density at the radiative-convective boundary (rcb), is the adiabatic gradient, is the adiabatic index of the interior, is the gravitational constant, is the sound speed evaluated at the location of the planet, is the equilibrium temperature of the planet, is the radius at the radiative-convective boundary, and is the structure integral that follows the form

(8) |

To eliminate , we use temperature gradient at the rcb so that

(9) |

where is the Stefan-Boltzmann constant, is the opacity at the rcb, and is the cooling luminosity, which can be written as

(10) |

where is the Kelvin-Helmholtz cooling time of the envelope, and again follows the structure integral given by equation 8. Substituting equation 10 into equation 9,

(11) |

By re-arranging equation 7, we find another equation for :

(12) |

We numerically solve for that obtains satisfying both equations 2.3 and 2.3, using the root_scalar function from SciPy optimize package.
Throughout the paper, we adopt ,^{1}^{1}1We note that at formation, the inner adiabat follows more closely as the energy is spent on dissociating hydrogen molecules. It is expected that approaches 7/5 as the envelope cools below the dissociation temperature 2500 K but this is yet to be verified with detailed, self-consistent calculation that tracks planets from their formation through post-disk evolution., , , and (Rogers & Seager, 2010).^{2}^{2}2These values for opacity are obtained by fitting to the tabulated opacity by Freedman et al. (2008), which is designed for dust-free atmospheres. In the absence of post-disk pollution by nearby small grains or giant impact, it is reasonable to consider the upper envelope to be drained out of grains (the gravitational settling timescale of a micron-sized grain is about 1 Myr).
To save computation time, we set for any that gives , motivated by the 5% error in Kepler transit depth measurement (e.g., Fulton & Petigura, 2018). This limit can be found easily by taking the limit of and confirming numerically:

(13) |

The photospheric radius—the observable—is a few scale height above . Correction for the photosphere is made using

(14) |

where is the density at the photosphere and is the surface gravity.

### 2.4 Envelope mass loss

Once the disk gas dissipates and the planets are laid bare to stellar insolation, those that are closest to the star are expected to lose their gaseous envelopes, either by photoevaporation (e.g., Owen & Wu, 2013) or by Parker wind (e.g., Ikoma & Hori, 2012; Owen & Wu, 2016; Ginzburg et al., 2018). The key difference between the two mechanisms is the source of insolation: whereas the former depends on the high-energy flux, the latter depends on the bolometric flux. As lower mass stars stay active for longer, photoevaporation model expects the radius-period gap to extend to longer orbital period, a hint of which is observed by Fulton & Petigura (2018, see their Figure 11). There is a discernible shift in the position of the gap towards larger radius around more massive host stars (Fulton & Petigura, 2018; Cloutier & Menou, 2020; Berger et al., 2020) . To reproduce this feature, photoevaporative model requires stellar-mass dependent core mass distribution (Wu, 2019) whereas this is a natural prediction of Parker wind, core-powered envelope mass loss model (Gupta & Schlichting, 2020). For solar-type stars, the two mechanisms predict similar location and shape of the gap in the radius-period distribution. Since the goal of this paper is to assess the likelihood of bottom-heavy core mass function for a fixed mass of the host star, we limit our analysis to photoevaporative mass loss for simplicity. We discuss potential effect of varying stellar mass in Section 4.

Following Owen & Wu (2017), we evolve the envelope mass over 5 Gyrs according to the energy-limited mass loss (e.g., Lopez & Fortney, 2013)

(15) |

where is the mass loss efficiency factor, and is the high-energy luminosity of the star (e.g., Ribas et al., 2005; Jackson et al., 2012)

(16) |

Orbital periods are drawn from the empirical distribution following Petigura et al. (2018)

(17) |

and then converted to orbital distance assuming solar mass host star.

## 3 Results

### 3.1 Primordial Radius Valley from Late-time Gas Accretion

We first show that late-time gas accretion alone produces a gap in the radius distribution (see Figure 1). The amount of gaseous envelope a core can accrete drops sharply below 1 as their gas masses are limited by the maximally cooled isothermal state. The exponential dependence of this isothermal envelope mass to the core mass (equation 5) implies a bimodal distribution of envelope mass fractions and therefore a bimodal distribution of radii, for a smooth, underlying core mass function (see Figure 2).

Figure 1 demonstrates that the location of the primordial “radius valley” shifts to smaller radii farther from the star. As the disk gets colder, planet’s Bondi radius increases and so the isothermal limit rises. Figure 2 illustrates this behavior where the rocky-to-enveloped transition shifts to smaller core masses at longer orbital periods. This negative slope of the valley in the radius-period space is reminiscent of that observed (Fulton et al., 2017; Van Eylen et al., 2018). We see a larger separation between the rocky and the enveloped planetary population for dust-free gas accretion. As Figure 2 shows, this difference arises from both the generally more rapid accretion and weaker dependence on core mass for dust-free envelopes.

As we will show in the next section, gas accretion needs to be dust-free in order for the primordial radius gap (and the post-evaporation gap) to align with the observation. From a numerical fit, we find the rocky-to-enveloped transition mass from dust-free accretion to scale with the disk temperature as . Since and , we find the radius valley , consistent within an errorbar of Van Eylen et al. (2018) and Martinez et al. (2019).

### 3.2 Mass Loss and Underlying Core Mass Distribution

Although the observed gap in the radius distribution and its dependence on orbital periods can be reproduced by late-time gas accretion, envelope mass loss is a natural next step once the disk gas completely dissipates. Figure 3 demonstrates that the location of the radius valley carved out by photoevaporative mass loss is robustly situated at 1.8 regardless of the primordial population. As Owen & Wu (2017) cogently explain, gas-enveloped planets transform to bare rocky cores by photoevaporation when their envelope mass loss timescale 100 Myrs, the typical saturation timescale of high-energy luminosity of host stars. For our choice of parameters, this transition occurs for 4–10 and 0.0004–0.002, corresponding to 1.8.

Where the initial conditions make a difference is in the depth and the width of the gap. As illustrated in Figure 3, the narrow valley and peak in the distribution of radii are more likely to appear in dust-free envelopes (blue lines) with smaller outer radius (smaller ) that are assembled in hotter disks (higher ), and built from less bottom-heavy core mass functions (smaller ).

The narrowness of the radius peak for dust-free envelopes as opposed to dusty envelopes can be understood from the weaker dependence of on (see equations 3 and 4 as well as Figure 2). For a given range of , the confines of possible envelope mass fractions and therefore radii are more limited.

Smaller reduces the maximum and so keeps the primordial radius peak closer to the valley. Since photoevaporative mass loss effectively carves out the large radii peak and add them to the lower radii, observations are better reproduced when the initial radius valley is narrower.

In hotter disks, the isothermal maximal shrinks so that the rocky-to-enveloped transition appears at higher core masses. The result is a positive shift in the location of the primordial radius valley. The gas accretion rate for dust-free envelopes also reduces (see equation 4) and so the primordial distribution of radii agrees well with the observation (see the faint blue line in the top middle panel of Figure 3). Since the locations of the valley are coincident with that expected from photoevaporative mass loss, we only observe slight reduction in the peak at 2.3 and a slight shallowing of the valley at 1.8.

We observe a loss of a peak in the radius distribution when the underlying core mass function is too bottom-heavy (). While we defer detailed formal fitting of models to the data for future analyses, it is already apparent that the allowed range of appears tightly constrained, under the ansatz that the core mass distribution follows a power-law. It may be possible to restore the radius peak even with with sufficiently high but we judge to be unlikely as it implies the disk is hot enough to melt iron at 0.1 AU.

The combination of parameters that provides the model radius distribution agreeing best with the observation are highlighted in Figures 4 and 5, corresponding to dust-free envelopes and (0.1, 3, 0.7) and (0.2, 3, 1), respectively. Between the primordial and evaporated population, we see a slight tilting of the slope in radius-period space but overall, the sign of this slope starts negative and ends negative, similar to that observed. To bring the primordial radius gap carved out in cooler disks to better alignment with the data, cores need to be slightly puffier. In Figure 6, we show a case with (0.1, 2, 0.7) with the core density set at 90% of the Earth, which we discuss in more detail in Section 4.

The observed radius valley closes at 10 days and widens towards 100 days (Fulton & Petigura, 2018). In photoevaporation models that assume all cores to have started with % by mass envelope, this traingular delta is hard to reproduce if the underlying core mass function is assumed flat (see, e.g., Owen & Wu, 2013). Figure 4 shows that the primordial population can recover the observed triangular shape of the radius gap. In hot disks, the rocky-to-enveloped transition mass rises while the envelope mass accreted by the core shrinks (see equation 4) so that the rocky and the enveloped populations “converge” at 10 days. This convergence erodes away for a logarithmically flat mass distribution (see Figure 5).

## 4 Discussion and Conclusion

We demonstrated that the underlying core mass distribution of sub-Neptunes can be broad with substantial population of sub-Earth mass objects while still reproducing the observed gap in the radius distribution and in radius-period space. A radius gap is already in place at birth as cores lighter than 1–2 can never accrete enough gas to be observed as gas-enveloped. The maximum envelope mass given by the maximally cooled isothermal state drops exponentially with core mass so that for a smooth distribution of core masses, a sharp radius dichotomy across 1–2 appears. Furthermore, this primordial radius gap shifts to smaller radii at longer orbital periods as the maximum isothermal mass rises and so the rocky-to-enveloped transition shifts to smaller cores.

Late-time formation of sub-Neptune is often attributed to producing a positive slope of the radius-period valley, based on the calculation of Lopez & Rice (2018). As these authors state, and we emphasize, their calculation is appropriate for formation in a gas-empty environment after a complete disk gas dispersal. The positive slope of the radius-period gap obtains from computing the expected core masses in a minimum mass extrasolar nebula (MMEN; Chiang & Laughlin, 2013) which produces rising masses (and therefore radii) at larger orbital distances (using the updated MMEN by Dai et al. (2020) will produce a similar result). The slope of the valley in the radius-period space may indeed turn positive around low mass stars (Cloutier & Menou 2020; but see Wu 2019). Our premise is distinct: we consider the formation of sub-Neptunes in gas-poor but not gas-empty nebula so that gas accretion, however limited, occurs. It is formation that is late-time in terms of the evolution of disk gas but not so late that there is no gas left (e.g., inner holes of transitional disks).

Our model of late-time gas accretion followed by photoevaporative mass loss best reproduces the observed location, width, and depth of the radius gap when the sub-Neptune cores follow mass functions shallower than or equal to and accrete dust-free gas in hot disks.^{3}^{3}3We find a potentially good agreement using a bottom-heavy core mass function for puffy cores but only in one-dimensional radius histogram. The triangular shape of the radius-period valley is challenging to reproduce with bottom-heavy core mass functions. The rate of accretion in dusty environment is too sensitive to core mass so that the final distribution of envelope mass fractions and therefore radii is too broad compared to that observed. The coagulation and the rain-out of dust grains may be an efficient process in sub-Neptune envelopes (Ormel, 2014).

### 4.1 Dependence on Disk Temperature and Stellar Mass

The required disk temperatures may be uncomfortably high. For accretion disks, the mid-plane temperature at 0.1 AU can be as high as 2000K (see D’Alessio et al., 1998, their Figure 3), consistent with our . Copious amount of dust in the upper layers of the disk could potentially increase the mid-plane temperature even further. Assuming the disk is optically thick, a factor of 5 enhancement in opacity (by e.g., high local dust-to-gas ratio) could be consistent with .

Even in colder disks, the location of the primordial radius valley can match the observation if the cores are slightly less dense, e.g., 90% of the Earth composition, consistent with what is reported for short-period super-Earths by Dorn et al. (2019) and more generally by Rogers & Owen (2020). Figure 6 demonstrates this agreement for , , and . We also note that the maximum core mass is set to 20 here. Shrinking the maximum core mass sharpens the radius peak at 2 but does not affect the location of the gap. We note that some of the model parameters that produce a broad peak may be narrowed by taking into account core-envelope interaction, in particular, the dissolution of gas into the magma core (Kite et al., 2019). Assessing the effect of core-envelope mixing at formation is a subject of our ongoing work.

We note that the primordial radius valley is expected to shift towards larger sizes around higher mass stars, assuming their disks are hotter. For disks heated by stellar irradiation, where is the effective temperature of the star and is the radius of the star, all evaluated for fully convective, pre-main sequence stars. For these passive disks, . Since the rocky-to-enveloped transition mass , and , the radius valley . For disks heated by accretion, . Taking (Calvet et al., 2004), we find , which corresponds to . Both estimates are within the 1- error bar estimate from Gaia-Kepler catalog by Berger et al. (2020). More accurate comparison will require better understanding of the thermal structure of the protoplanetary disks and their dependence on the host stellar mass.

### 4.2 Primordial vs. Mass Loss

Late-time gas accretion alone can reproduce the observed shape of the super-Earth/sub-Neptune radius-period distribution. Furthermore, the fact that cores smaller than 1–2 cannot accrete enough nebular gas to appear as enveloped open up the possibility that the underlying core mass distribution can be broader than previously reported, extending well into the sub-Earth regime.

Nevertheless, mass loss processes are the natural outcome after complete dispersal of disk gas, whether by photoevaporation or by core-powered envelope mass loss via Parker wind. Precise characterization of planet-hosting stars with Gaia find a growth in the super-Earth population in old (1 Gyr) vs. young (1 Gyr) stars (Berger et al., 2020), suggesting long-term mass loss processes continue to shape the overall exoplanet radius distribution.

There remain uncertainties in the exact magnitude of the mass loss for both photoevaporation and core-powered mass loss models. In the picture of photoevaporation, the unknown strength of planetary magnetic fields can shield against high-energy stellar photons (Owen & Adams, 2019). Furthermore, there is an order of magnitude variation in the magnitude and time evolution of stellar EUV and X-ray luminosity (Tu et al., 2015). In the picture of core-powered envelope mass loss, the amount of gas mass that can be lost via wind depends on the structure of the outer envelope subject to uncertain opacity sources. Even if the cores hold enough thermal energy to unbind the entire envelope, the timescale of heat transfer depends on the unknown viscosity and Prandtl number of the magma/rocky core (e.g., Stamenković et al., 2012). Further advances in both theory and observations should iron out these uncertainties.

## References

- Ali-Dib et al. (2020) Ali-Dib, M., Cumming, A., & Lin, D. N. C. 2020, MNRAS, 494, 2440, doi: 10.1093/mnras/staa914
- Berger et al. (2020) Berger, T. A., Huber, D., Gaidos, E., van Saders, J. L., & Weiss, L. M. 2020, arXiv e-prints, arXiv:2005.14671. https://arxiv.org/abs/2005.14671
- Burke et al. (2015) Burke, C. J., Christiansen, J. L., Mullally, F., et al. 2015, ApJ, 809, 8, doi: 10.1088/0004-637X/809/1/8
- Calvet et al. (2004) Calvet, N., Muzerolle, J., Briceño, C., et al. 2004, AJ, 128, 1294, doi: 10.1086/422733
- Chiang & Laughlin (2013) Chiang, E., & Laughlin, G. 2013, MNRAS, 431, 3444, doi: 10.1093/mnras/stt424
- Choksi & Chiang (2020) Choksi, N., & Chiang, E. 2020, MNRAS, 495, 4192, doi: 10.1093/mnras/staa1421
- Cloutier & Menou (2020) Cloutier, R., & Menou, K. 2020, AJ, 159, 211, doi: 10.3847/1538-3881/ab8237
- Dai et al. (2020) Dai, F., Winn, J. N., Schlaufman, K., et al. 2020, AJ, 159, 247, doi: 10.3847/1538-3881/ab88b8
- D’Alessio et al. (1998) D’Alessio, P., Cantö, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411, doi: 10.1086/305702
- Dorn et al. (2019) Dorn, C., Harrison, J. H. D., Bonsor, A., & Hands, T. O. 2019, MNRAS, 484, 712, doi: 10.1093/mnras/sty3435
- Freedman et al. (2008) Freedman, R. S., Marley, M. S., & Lodders, K. 2008, ApJS, 174, 504, doi: 10.1086/521793
- Fressin et al. (2013) Fressin, F., Torres, G., Charbonneau, D., et al. 2013, ApJ, 766, 81, doi: 10.1088/0004-637X/766/2/81
- Fulton & Petigura (2018) Fulton, B. J., & Petigura, E. A. 2018, AJ, 156, 264, doi: 10.3847/1538-3881/aae828
- Fulton et al. (2017) Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109, doi: 10.3847/1538-3881/aa80eb
- Fung et al. (2019) Fung, J., Zhu, Z., & Chiang, E. 2019, ApJ, 887, 152, doi: 10.3847/1538-4357/ab53da
- Ginzburg et al. (2016) Ginzburg, S., Schlichting, H. E., & Sari, R. 2016, ApJ, 825, 29, doi: 10.3847/0004-637X/825/1/29
- Ginzburg et al. (2018) —. 2018, MNRAS, 476, 759, doi: 10.1093/mnras/sty290
- Gupta & Schlichting (2020) Gupta, A., & Schlichting, H. E. 2020, MNRAS, 493, 792, doi: 10.1093/mnras/staa315
- Hadden & Lithwick (2014) Hadden, S., & Lithwick, Y. 2014, ApJ, 787, 80, doi: 10.1088/0004-637X/787/1/80
- Hsu et al. (2019) Hsu, D. C., Ford, E. B., Ragozzine, D., & Ashby, K. 2019, AJ, 158, 109, doi: 10.3847/1538-3881/ab31ab
- Ikoma & Hori (2012) Ikoma, M., & Hori, Y. 2012, ApJ, 753, 66, doi: 10.1088/0004-637X/753/1/66
- Jackson et al. (2012) Jackson, A. P., Davis, T. A., & Wheatley, P. J. 2012, MNRAS, 422, 2024, doi: 10.1111/j.1365-2966.2012.20657.x
- Kite et al. (2019) Kite, E. S., Fegley, Bruce, J., Schaefer, L., & Ford, E. B. 2019, ApJ, 887, L33, doi: 10.3847/2041-8213/ab59d9
- Lambrechts & Lega (2017) Lambrechts, M., & Lega, E. 2017, A&A, 606, A146, doi: 10.1051/0004-6361/201731014
- Lee & Chiang (2015) Lee, E. J., & Chiang, E. 2015, ApJ, 811, 41, doi: 10.1088/0004-637X/811/1/41
- Lee & Chiang (2016) —. 2016, ApJ, 817, 90, doi: 10.3847/0004-637X/817/2/90
- Lee et al. (2014) Lee, E. J., Chiang, E., & Ormel, C. W. 2014, ApJ, 797, 95, doi: 10.1088/0004-637X/797/2/95
- Lopez & Fortney (2013) Lopez, E. D., & Fortney, J. J. 2013, ApJ, 776, 2, doi: 10.1088/0004-637X/776/1/2
- Lopez & Fortney (2014) —. 2014, ApJ, 792, 1, doi: 10.1088/0004-637X/792/1/1
- Lopez & Rice (2018) Lopez, E. D., & Rice, K. 2018, MNRAS, 479, 5303, doi: 10.1093/mnras/sty1707
- Malhotra (2015) Malhotra, R. 2015, ApJ, 808, 71, doi: 10.1088/0004-637X/808/1/71
- Marcy et al. (2014) Marcy, G. W., Isaacson, H., Howard, A. W., et al. 2014, ApJS, 210, 20, doi: 10.1088/0067-0049/210/2/20
- Martinez et al. (2019) Martinez, C. F., Cunha, K., Ghezzi, L., & Smith, V. V. 2019, ApJ, 875, 29, doi: 10.3847/1538-4357/ab0d93
- Ormel (2014) Ormel, C. W. 2014, ApJ, 789, L18, doi: 10.1088/2041-8205/789/1/L18
- Owen & Adams (2019) Owen, J. E., & Adams, F. C. 2019, MNRAS, 490, 15, doi: 10.1093/mnras/stz2601
- Owen & Wu (2013) Owen, J. E., & Wu, Y. 2013, ApJ, 775, 105, doi: 10.1088/0004-637X/775/2/105
- Owen & Wu (2016) —. 2016, ApJ, 817, 107, doi: 10.3847/0004-637X/817/2/107
- Owen & Wu (2017) —. 2017, ApJ, 847, 29, doi: 10.3847/1538-4357/aa890a
- Petigura et al. (2013) Petigura, E. A., Howard, A. W., & Marcy, G. W. 2013, Proceedings of the National Academy of Science, 110, 19273, doi: 10.1073/pnas.1319909110
- Petigura et al. (2018) Petigura, E. A., Marcy, G. W., Winn, J. N., et al. 2018, AJ, 155, 89, doi: 10.3847/1538-3881/aaa54c
- Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62, doi: 10.1006/icar.1996.0190
- Ribas et al. (2005) Ribas, I., Guinan, E. F., Güdel, M., & Audard, M. 2005, ApJ, 622, 680, doi: 10.1086/427977
- Rogers & Owen (2020) Rogers, J. G., & Owen, J. E. 2020, arXiv e-prints, arXiv:2007.11006. https://arxiv.org/abs/2007.11006
- Rogers & Seager (2010) Rogers, L. A., & Seager, S. 2010, ApJ, 712, 974, doi: 10.1088/0004-637X/712/2/974
- Stamenković et al. (2012) Stamenković, V., Noack, L., Breuer, D., & Spohn, T. 2012, ApJ, 748, 41, doi: 10.1088/0004-637X/748/1/41
- Tu et al. (2015) Tu, L., Johnstone, C. P., Güdel, M., & Lammer, H. 2015, A&A, 577, L3, doi: 10.1051/0004-6361/201526146
- Valencia et al. (2006) Valencia, D., O’Connell, R. J., & Sasselov, D. 2006, Icarus, 181, 545, doi: 10.1016/j.icarus.2005.11.021
- Van Eylen et al. (2018) Van Eylen, V., Agentoft, C., Lundkvist, M. S., et al. 2018, MNRAS, 479, 4786, doi: 10.1093/mnras/sty1783
- Weiss & Marcy (2014) Weiss, L. M., & Marcy, G. W. 2014, ApJ, 783, L6, doi: 10.1088/2041-8205/783/1/L6
- Wu (2019) Wu, Y. 2019, ApJ, 874, 91, doi: 10.3847/1538-4357/ab06f8
- Wu & Lithwick (2013) Wu, Y., & Lithwick, Y. 2013, ApJ, 772, 74, doi: 10.1088/0004-637X/772/1/74