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

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 $\gtrsim 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 $\gtrsim 100\,{M}_{\oplus }$, but few intermediate-mass giant planets of ∼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

$$\begin{eqnarray}\begin{array}{rcl}{R}_{{\rm{E}}} & = & {D}_{{\rm{L}}}\sqrt{\displaystyle \frac{4{{GM}}_{{\rm{L}}}}{{c}^{2}}\left(\displaystyle \frac{1}{{D}_{{\rm{L}}}}-\displaystyle \frac{1}{{D}_{{\rm{S}}}}\right)}\\ & = & 4.04\,\mathrm{au}\sqrt{\displaystyle \frac{{M}_{{\rm{L}}}}{{M}_{\odot }}\displaystyle \frac{{D}_{S}}{8\,\mathrm{kpc}}\,4x(1-x)},\end{array}\end{eqnarray} \tag{ 1 }$$

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\,\mathrm{au}\sqrt{{M}_{{\rm{L}}}/{M}_{\odot }}$, which is larger than the snow line at $\sim 2.7\,\mathrm{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).

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

2.2. Bern Model

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 Figures 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\gt 0.01$, which gives rise to a slight decrease in sensitivity for q values just above 0.01.

Zoom In Zoom Out Reset image size

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.6{M}_{\odot }$, as explained in S16.

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 timescale (${\tau }_{\mathrm{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).

Zoom In Zoom Out Reset image size

Comparison between the population synthesis results and the S16 data shows some significant differences in the planetary occurrence rate for all mass ratios $q\geqslant {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}\leqslant q\leqslant {10}^{-3}$, and it decreases to a factor ∼5 for higher mass-ratios, ${10}^{-3}\lt q\leqslant 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}\leqslant q\leqslant {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 ∼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\mbox{--}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.

Table 1.  Poisson Probability for $N\geqslant {N}_{\mathrm{obs}}$ for Mass-ratio Bins with $-4.45\lt \mathrm{log}q\lt -2.95$

	I & L Model		Bern Model
$\mathrm{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}$

Download table as:  ASCIITypeset image

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.6\,{M}_{\odot }$, this mass-ratio desert in the population synthesis models corresponds to planets with masses $20\mbox{--}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 timescale 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 timescale, 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 (α) 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.

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) α 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 α (below 10⁻⁴). 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 α as the default models) around $0.5\,{M}_{\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 2(b)) yield compatibility with the microlensing data in the ${10}^{-4}\lt q\lt {10}^{-3}$ mass-ratio range. Lowering of protoplanets' tidal truncation mass also leads to the formation of too many planets with $q\lt 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\gt {10}^{-3}$. These issues with the low-α 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 a 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 α, h, and ${\tau }_{\mathrm{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.