Microlensing Results Challenge the Core Accretion Runaway Growth Scenario for Gas Giants

TL;DR

Microlensing observations challenge core accretion models by showing more intermediate-mass gas giants than predicted, suggesting alternative formation processes or environmental variations.

Contribution

This study compares microlensing data with core accretion models, revealing a significant discrepancy in predicted and observed planet mass distributions.

Findings

Microlensing detects more intermediate-mass gas giants than models predict.

Core accretion models underestimate planets with mass ratios of 10^{-4} to 4×10^{-4}.

Discrepancies suggest alternative formation processes or environmental factors.

Abstract

We compare the planet-to-star mass-ratio distribution measured by gravitational microlensing to core accretion theory predictions from population synthesis models. The core accretion theory's runaway gas accretion process predicts a dearth of intermediate-mass giant planets that is not seen in the microlensing results. In particular, the models predict $\sim10\,\times$ fewer planets at mass ratios of $10^{-4} \leq q \leq 4 \times 10^{-4}$ than inferred from microlensing observations. This tension implies that gas giant formation may involve processes that have hitherto been overlooked by existing core accretion models or that the planet-forming environment varies considerably as a function of host-star mass. Variation from the usual assumptions for the protoplanetary disk viscosity and thickness could reduce this discrepancy, but such changes might conflict with microlensing results at larger or smaller mass ratios, or with other observations. The resolution of this discrepancy may have important implications for planetary habitability because it has been suggested that the runaway gas accretion process may have triggered the delivery of water to our inner solar system. So, an understanding of giant planet formation may help us to determine the occurrence rate of habitable planets.

Tables

Table 1. Table 1: Poisson Probability for N ≥ N obs 𝑁 subscript 𝑁 obs N\geq N_{\rm obs} for Mass-ratio Bins with − 4.45 < log ​ q < − 2.95 4.45 log 𝑞 2.95 -4.45<{\rm log}q<-2.95

	I & L model		Bern model
$\log q$	Migration	No	Migration	No
$-4.2$	$0.14$	$2.0\times {10}^{-2}$	$1.5\times {10}^{-3}$	$4.4\times {10}^{-2}$
$-3.7$	$1.1\times {10}^{-5}$	$6.0\times {10}^{-6}$	$1.7\times {10}^{-8}$	$2.3\times {10}^{-4}$
$-3.2$	$3.6\times {10}^{-6}$	$2.6\times {10}^{-2}$	$1.8\times {10}^{-9}$	$1.3\times {10}^{-3}$

Equations

$$R_{\mathrm{E}}=D_{\mathrm{L}}\sqrt{\frac{4G{M}_{\mathrm{L}}}{c^2}\left(\frac{1}{D_{\mathrm{L}}}-\frac{1}{D_{\mathrm{S}}}\right)}=4.04\mathrm{AU}\sqrt{\frac{M_{\mathrm{L}}}{M_{\odot }}\frac{D_S}{8\mathrm{kpc}}4x(1-x)},$$

Full text

gravitational lensing: micro — planetary systems — planets and satellites: formation — planet-disk interactions

1 Introduction

The core accretion model (Pollack et al., 1996) of planet formation predicts a deficit of planets between the masses of Neptune and Saturn (Ida & Lin, 2004a; Mordasini et al., 2009). This desert is a consequence of the runaway accretion of hydrogen and helium gas onto protoplanetary cores that have attained a critical mass of $\sim 10\,M_{\oplus}$. These cores are preferentially formed outside the snow line. The rapid growth of gaseous envelopes around the cores is quenched by either the severe gas depletion throughout the planetary disk before the runaway accretion can begin or the growth of \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}100\,M_{\oplus} gas giants that clear the gaps in the disk. This process is expected to produce numerous $\sim 10\,M_{\oplus}$ “failed gas giant cores”, as well as gas giants of \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}100M_{\oplus}, but few intermediate-mass giant planets of $\sim 20$-$80\,{M_{\oplus}}$. In this Letter, we test this scenario with gravitational microlensing observations.

A unique strength of the microlensing method is its sensitivity to low-mass planets (Bennett & Rhie, 1996) in Jupiter-like orbits (Gould & Loeb, 1992), beyond the snow line. Microlensing is most sensitive to planets located at a projected separation similar to the Einstein radius, which is given by

[TABLE]

where $M_{\rm L}$ is the lens mass, $G$ is the gravitational constant, $c$ is the speed of light, $D_{\rm L}$ and $D_{\rm S}$ are the lens and source distances, and $x=D_{\rm L}/D_{\rm S}$. For a typical value of $4x(1-x)=0.75$ and $D_{\rm S}=8\,$kpc, this gives $R_{\rm E}=3.5\,{\rm AU}\sqrt{M_{\rm L}/M_{\odot}}$, which is larger than the snow line at $\sim 2.7\,{\rm AU}(M_{\rm L}/M_{\odot})$ for stars of a solar mass or less.

In this Letter, we compare the cold-planet mass-ratio function derived by the MOA collaboration (Suzuki et al., 2016) (hereafter S16), using 30 microlens planets, to predictions based on the core accretion theory (Pollack et al., 1996), which was originally constructed to study the formation of our own solar system. This theory involves many competing physical processes including the formation of protoplanetary embryos as progenitor of terrestrial planets, gas accretion onto cores of proto-gas giants, the clearing of gaps in the protoplanetary disks, and planetary migration induced by planets’ interaction with their natal disks. None of these processes can be calculated with certainty. Nevertheless, it is possible to introduce a set of approximations of these processes to simulate the statistical distribution for an ensemble of emerging systems from protostellar disks with an assumed evolving distribution of H/He gas and planetesimals. These constructs generate population syntheses that can be compared with observations. We compare the microlensing results to the population synthesis models of two different groups: Ida & Lin (2004a) and the Bern group (Mordasini et al., 2009).

2 Population Synthesis Calculations

Population synthesis models generate a set of simulated planetary systems that have a distribution of planets in the semimajor-axis, planet-mass plane. The planet distribution is thought to depend on the host-star mass, and so we have run simulations with different host-star masses ranging from $0.08$ to $1.0\,M_{\odot}$. We select a range of host-star masses in logarithmically uniform bins to span this mass range, with fewer bins for the Bern models, because the Bern code uses more CPU time. For each synthesis model, we generate planetary systems with and without the planetary migration effect because, as discussed below, the models with migration predicted many fewer planets than found by microlensing over a large range of mass ratios. The details of these population synthesis models are explained below.

2.1 IL Model

The detailed description of the Ida & Lin (IL) model is written in a series of papers (Ida & Lin, 2004a, b, 2005, 2008a, 2008b, 2010; Ida et al., 2013). The synthesis model includes classical models of planetary growth and migration from small planetesimals, combining planetesimal accretion, gas accretion onto the planet, type I and II migration, as well as planet-planet scattering between all planets. Planetary embryos are set with an initial mass of $10^{20}\,{\rm g}$ with orbital separations of $\sim 10\,\times$ the Hill radius of the classical isolation mass (Kokubo & Ida, 1998) at separations of $0.05-20\,{\rm au}$. The self-similar disk model is used for disk evolution. The gas surface density at 10 au is distributed in a range of $[0.1,10]$ times of the minimum-mass solar nebula model (Hayashi et al., 1985) with a log-normal distribution. The initial metallicity of the disk is distributed in a range of $[-0.2,0.2]$ dex with a normal distribution. The IL simulations were done with nine different host masses given by $M_{h}=10^{\gamma}M_{\odot}$, where $\gamma=-1.3,-1.15,-1.0,-0.85,-0.7,...,-0.1$, and the average weights used for these simulations were 0.067, 0.047, 0.063, 0.077, 0.093, 0.111, 0.135, 0.163, and 0.243, respectively. Exactly 10,000 simulations were run for each mass ratio.

2.2 Bern Model

The Bern model is a global planet formation (Alibert et al., 2005) and evolution (Mordasini et al., 2012) model that is based on core accretion. It predicts the properties of emerging planetary systems like the masses and orbits of the planets based directly on the properties of the parent protoplanetary disk (such as the disk mass, dust-to-gas ratio, and lifetime). As described in recent reviews (Benz et al., 2014; Mordasini et al., 2015), the model couples self-consistently several standard astrophysical submodels for planetary formation to compute the evolution of the gas and planetesimals disks, the accretion of gas and solids by the protoplanets, as well as interactions between the protoplanets (gravitational N-body interaction) and between the gas disk and the protoplanets (orbital migration). More specifically, the model consists of the following elements:

a protoplanetary gas disk modeled by the numerical solution of the classical 1D viscous evolution equation for the gas surface density in an axisymmetric constant $\alpha$-viscosity disk with stellar irradiation (Hueso & Guillot, 2005) and photoevaporation (Hollenbach et al., 1994); 2. 2.

the protoplanets’ accretion rates of solids modeled by a Safronov-type rate equation from planetesimals of a single size in the oligarchic growth regime (Inaba et al., 2001); 3. 3.

the disk of planetesimals modeled as a surface density with a dynamical state (Fortier et al., 2013); 4. 4.

the protoplanets’ gas accretion rate and planetary interior structure obtained from solving the standard 1D radially symmetric hydrostatic planet interior structure equations; 5. 5.

planetary orbital migration modeled as gas disk-driven non-isothermal Type I and Type II migration (Dittkrist et al., 2014); and 6. 6.

dynamical interactions between the protoplanets modeled with the explicit N-body integrator Mercury (Chambers, 1999).

In population syntheses models, the initial conditions of the model, which are the properties of the protoplanetary disks, are varied according to observed distribution of protoplanetary disk properties, and the global model is run typically several hundred times in order to synthesize the predicted population of model planets. The Bern group population synthesis calculations are much more computationally intensive than the IL calculations, so these were done at host masses of $M_{\rm h}=10^{\gamma}\,M_{\odot}$, where $\gamma=-0.903,-0.602,-0.301$, and 0, and the average weights for these masses were 0.268, 0.211, 0.308, and 0.213, respectively. The number of simulations run for each mass ratio was 1387, 4805, 4093, and 1392, respectively, for the standard simulations and 1150, 1717, 1181, and 1918, respectively, for the migration-free simulations.

3 Cold Exoplanet Mass-ratio Function from Microlensing

Our analysis is based on the S16 statistical analysis (Suzuki et al., 2016) of the combined sample of the 1474 well-characterized microlensing events from the MOA survey sample combined with earlier, smaller microlensing samples (Gould et al., 2010; Cassan et al., 2012). This sample includes 30 planets with mass ratios, $q$, and separations, $s$, that are displayed in Figure 1 and 2. Figure 1 shows the survey sensitivity for the combined survey, which is the sum of the detection efficiencies for all the events in the survey. The planets in this S16 sample are indicated by red dots, and planets with uncertainty in the separation, $s$, due to the close-wide degeneracy of high-magnification events are indicated as open red circles connected by red lines. The contours in this figure indicate the number of planets that would be detected if each event had a planet with the specified $s$ and $q$ values. The sensitivity to planets generally increases at larger $q$ values, but the older microlensing samples (Gould et al., 2010; Cassan et al., 2012) did not consider planets with mass ratios $q>0.01$, which gives rise to a slight decrease in sensitivity for $q$ values just above 0.01.

It is this survey sensitivity that is used to convert the microlensing results to the power-law mass-ratio functions investigated in S16, but our comparison to population synthesis models requires a slightly different approach. The population synthesis models provide a set of simulated planetary systems for each assumed host-star mass. The planet-star mass-ratio values are determined directly from these simulations, but the projected separations, $s$ in Einstein radii are not directly produced by the population synthesis calculations. We use a standard Galactic model (Han & Gould, 1995) to produce a probability distribution of primary lens masses for each of the 1474 microlensing events in the S16 sample. Then, for each of these events, we run 4000 random trials for the IL and Bern group simulations for each of the 1474. In each trial, we randomly select a lens distance and a host-mass bin based on the lens mass and distance probability distribution for the event under consideration, and we randomly select one of the simulated planetary systems from that mass bin. Next, we select a random orientation for that planetary system to determine the $s$ value for that event. Finally, for the trial for each of the S16 events, we apply the S16 detection efficiency as a function of $q$ and $s$ for that event to determine if the simulated planets are detected. This is equivalent to simulating 4000 trials of the S16 observations, and the total number of events simulated is $4000\times 1474=5.9\times 10^{6}$.

This procedure automatically selects planetary host-star masses from the distributions expected for our sample of microlensing events, using our assumed Galactic model. The microlensing rate imparts a weight that scales as $\sqrt{M_{\rm L}}$, and there is an additional weighting from the microlensing event and planet detection efficiencies. These effects push the median mass of the host stars produced by our survey up to about $0.6M_{\odot}$, as explained in S16.

4 Predicted Microlens Mass-ratio Function from Population Synthesis

Figure 2 compares the S16 results to the population synthesis models by IL (a) and the Bern group (b). Both sets of default models assumed the widely adopted values for the effective viscosities ($\alpha=10^{-3}$), aspect ratio ($h=H/r$), heavy element abundance, mass flux (${\dot{M}}$), and depletion time scale ($\tau_{\rm dep}$), for the natal disks of protoplanets (Hartmann, 1998). Despite some differences in their prescriptions for the growth of planetesimals and gas accretion onto gas giants, the default IL and Bern models similarly reproduce the mass-period distribution for gas giants around solar-type stars found by the early radial velocity (RV) surveys (Cumming et al., 2008). In order to account for a population of short-period gas giants found by these observations, these default models include the effect of planetary migration (Lin & Papaloizou, 1986).

Comparison between the population synthesis results and the S16 data shows some significant differences in the planetary occurrence rate for all mass ratios $q\geq 10^{-4}$, up to the brown dwarf boundary at $q=0.03$. This discrepancy is a factor of 10–25 at the lower mass-ratios, $10^{-4}\leq q\leq 10^{-3}$, and it decreases to a factor $\sim$5 for higher mass-ratios, $10^{-3}<q\leq 0.03$. One possible culprit for this discrepancy may be the loss of such planets due to both type I and II migration. This conjecture is verified by another series of simulations with identical disk-model parameters but without planetary migration. For mass ratios of $10^{-3}\leq q\leq 10^{-2}$, these migration-free simulations with the Bern models provide a good match to the S16 data, but the IL models produce an overretention of the same population by a factor of $\sim$6.

The migration-free approximation also provides a correction of an overestimation of the planetary migration rate (Lin & Papaloizou, 1986) and reduces an excess of hot Jupiters predicted by the population synthesis models (Ida & Lin, 2004a). Such a modification may be due to either diffusion across the gap (Fung et al., 2014), or relative low viscosity (Duffell & MacFadyen, 2013) (see below). However, these migration-free models continue to generate, respectively, 10 and 7 times fewer planets with $q=1-4\times 10^{-4}$ than the S16 data. Figure 2 contains four, seven, and eight planets in the mass-ratio bin centered at $q=6.3\times 10^{-5}$, $2.0\times 10^{-4}$, and $6.3\times 10^{-4}$, respectively. Based on their large deficits in Figure 2, we show in Table 1 the low Poisson probability for the predicted number of the planet detections from the IL and Bern models to match or exceed the actual number of planet detections in S16.

Both the simulated and observed S16 samples obey Poisson statistics, but the statistical uncertainty on the simulated planets are negligible compared to the Poisson uncertainties for the observed sample. So, we use the Poisson statistics for the S16 events only to compare the data to the simulations. These numbers obey Poisson statistics and are used for the computation of the Poisson probabilities reported in Table 1. As mentioned above, Table 1 shows the Poisson probabilities for the three $q$ bins centered on $q=2.0\times 10^{-4}$ for the standard and migration-free simulations by both groups. Population synthesis models fail this comparison for the bin centered at $q=2.0\times 10^{-4}$ for each of the IL and Bern population synthesis models. The two adjacent bins, centered at $q=6.3\times 10^{-5}$ and $q=6.3\times 10^{-4}$ also fail this comparison in almost every case.

Since the median mass of host stars probed by the microlensing survey is $M_{\rm L}\approx 0.6M_{\odot}$, this mass-ratio desert in the population synthesis models corresponds to planets with masses $20-80{M_{\oplus}}$. According to the core accretion scenario, solid planetesimals coagulate into cores that begin to gradually accumulate gaseous envelopes after their mass exceeds $\sim 10{M_{\oplus}}$ (Pollack et al., 1996). Runaway gas accretion is initiated after the envelopes’ mass becomes comparable to that of the cores. The planets’ growth rate accelerates on a time scale inversely proportional to their mass until gas is severely reduced either throughout the disk or in the proximity of their orbit. Since the duration of runaway growth is generally shorter than the global disk depletion time scale, planets form less frequently in the $q$ desert, which is bounded by the critical-core mass for the onset of rapid gas accretion and the mass for the formation of tidally induced gap (Lin & Papaloizou, 1980). Gas giants’ asymptotic $q$ is determined by the magnitude of viscosity ($\alpha$) and aspect ratio ($h$) (Lin & Papaloizou, 1993). For the disk parameters we adopted in both the default and migration-free models, it exceeds $4\times 10^{-4}$ near the snow line.

5 Discussion and Conclusions

The incompatibility between the default population synthesis models and the exoplanet mass-ratio function measured by microlensing suggests that some of the assumptions we have adopted for the runaway gas accretion scenario might be incorrect or incomplete. We note that a similar contradiction between the core accretion prediction of a sub-Saturn-mass desert at a short-period orbit (Ida & Lin, 2004b) is also not seen in the Kepler data (Thompson et al., 2018), but the contradiction seems more significant at separations beyond the show line, where gas giant formation is thought to occur. Perhaps there are some physical processes that can suppress or quench the gas accretion rate. Planetesimals might be captured by the accreting gaseous envelope and heat it up, which would slow the accretion (Alibert et al., 2018). Gas giant planets’ asymptotic mass can also be lowered in disk regions with smaller (Dobbs-Dixon et al., 2007; Fung et al., 2014) $\alpha$ and $h$ than the values we have adopted in the default and migration-free models. During the advanced stages of their viscous evolution, protostellar disks may indeed have smaller $h$ values than that for the minimum mass nebula model (Garaud & Lin, 2007), especially around low-mass stars. Growing protoplanets may also undergo type III migration (Masset et al., 2006) or gravitational scattering to regions beyond the snow line where MHD turbulence may be suppressed by magnetic diffusivity (Bai, 2017) to reduce the magnitude of effective $\alpha$ (below $10^{-4}$). These regions may be manifested in the form of narrow and axisymmetric rings and gaps in protostellar disks commonly found by ALMA and they may provide nests for the low-mass protoplanets. These effects need to be explored further elsewhere.

The Bern group has performed some preliminary simulations for low-viscosity disks (with one-third the value of $\alpha$ as the default models) around $0.5M_{\odot}$ host stars. In addition, the tidal interaction of protoplanets with the protoplanetary gas disk can lead to gap opening, which in turn reduces the planet’s gas accretion rate by reducing the surface density of the gas in the vicinity of the planet (D’Angelo & Lubow, 2010). These results (represented by the gold histogram in Figure 2b) yield compatibility with the microlensing data in the $10^{-4}<q<10^{-3}$ mass-ratio range. Lowering of protoplanets’ tidal truncation mass also leads to the formation of too many planets with $q<3\times 10^{-5}$ such that center of the $q$ desert is shifted to lower $q$ values with less paucity. However, these model parameters also suppress the emergence of planets with $q>10^{-3}$. These issues with the low-$\alpha$ models are verified by analogous IL simulations that also fail to reproduce the population of Jupiter-mass planets found by the RV surveys around solar-type stars. However, this difficulty might be mitigated if some higher mass-ratio planets could be formed by mechanism other than core accretion, such as gravitational instability in the protoplanetary disk (Durisen et al., 2007; Forgan & Rice, 2013; Boss, 2017).

Since the microlensing and RV surveys sample host stars with somewhat different, but overlapping, mass distributions, a stellar-mass-dependent $q$ distribution (due to variations in the values of $\alpha$, $h$, and $\tau_{\rm dep}$) remains a possibility. In principle, gaps in the $q$ distribution may be smoothed out when the survey samples include a range of stellar masses. We plan to test this possibility with high angular resolution follow-up observations of the S16 planetary microlensing events using adaptive optics (AO) observations on the Keck telescopes under a recently approved NASA Key Strategic Mission Support (KSMS) program and Hubble Space Telescope observations (Bhattacharya et al., 2018). When the host star is detected, its mass can usually be determined from a measurement of the host-star brightness combined with constraints from the microlensing light curve. In these follow-up observations, it is important to measure the separation of the lens and source stars (Batista et al., 2015; Bennett et al., 2015) to confirm the lens star identification because confusion with other stars, like a binary companion to the lens or source, is possible (Bhattacharya et al., 2017; Koshimoto et al., 2017).

One of the first targets observed under the Keck KSMS program is the S16 event OGLE-2012-BLG-0950, which has been found (Bhattacharya et al., 2018) to have a host mass of $M_{h}=0.58\pm 0.04\,M_{\odot}$ and a planet mass of $39\pm 8\,M_{\oplus}$. This planet is right in the middle of the mass gap expected from the runaway gas accretion process in the default and migration-free models. There are also several solar-type stars outside the S16 sample that have planets in this intermediate-mass range, including the microlens planet OGLE-2012-BLG-0026Lc (Beaulieu et al., 2016), with $M=46.0\pm 2.6\,M_{\oplus}$, and two planets from the HARPS survey (Mayor et al., 2011) just inside the snow line with $M\sin i\sim 50\,M_{\oplus}$. If the mass measurements of the other planets in the S16 sample reveal a number of other planets in this predicted mass gap, this would rule out the host-mass dependence as the reason for the discrepancy between the microlensing data and the expectations for the runaway gas accretion process. The resolution of this discrepancy may help us to understand habitability of inner planets because it has been suggested that the delivery of water to the inner planets of our solar system may be a consequence of the runaway gas accretion process (Raymond & Izidoro, 2017). Ultimately, the exoplanet survey (Bennett & Rhie, 2002) of the Wide Field Infrared Survey Telescope (Penny et al., 2018) will perform a much more comprehensive survey that will probe planets beyond the orbital separation of Venus with sensitivity down to the mass of Mars ($0.1\,M_{\oplus}$).

D.P.B., A.B., and C.R. were supported by NASA through grant NASA-80NSSC18K0274. The MOA project is supported by JSPS KAKENHI grant Nos. JSPS24253004, JSPS26247023, JSPS23340064, JSPS15H00781, JP16H06287 and JSPS17H02871. The work by C.R. was supported by an appointment to the NASA Postdoctoral Program at the Goddard Space Flight Center, administered by USRA through a contract with NASA. C.M. acknowledges the support from the Swiss National Science Foundation under grant BSSGI0_155816, “PlanetsInTime.” Parts of this work have been carried out within the frame of the National Center for Competence in Research PlanetS supported by the SNSF. Work by N.K. was supported by JSPS KAKENHI grant No. JP18J00897. N.J.R. is a Royal Society of New Zealand Rutherford Discovery Fellow.