Are the observed gaps in protoplanetary discs caused by growing planets?

TL;DR

This study investigates whether the gaps observed in protoplanetary discs are caused by growing planets, using population synthesis models to compare theoretical predictions with observations, and finds discrepancies suggesting alternative explanations or missing physics.

Contribution

The paper introduces a comprehensive planet population synthesis model including pebble and gas accretion, migration, and compares the results with observed disc gaps and exoplanet populations.

Findings

High pebble mass needed for gap formation, only possible in metal-rich stars.

Formed planets predominantly become gas giants, inconsistent with observed super-Earths.

Discrepancies suggest gaps may not always be caused by planets or missing physics in models.

Abstract

Recent detailed observations of protoplanetary discs revealed a lot of sub-structures which are mostly ring-like. One interpretation is that these rings are caused by growing planets. These potential planets are not yet opening very deep gaps in their discs. These planets instead form small gaps in the discs to generate small pressure bumps exterior to their orbits that stop the inflow of the largest dust particles. In the pebble accretion paradigm, this planetary mass corresponds to the pebble isolation mass, where pebble accretion stops and efficient gas accretion starts. We perform planet population synthesis via pebble and gas accretion including type-I and type-II migration. In the first stage of our simulations, we investigate the conditions necessary for planets to reach the pebble isolation mass and compare their position to the observed gaps. We find that in order to match the gap structures 2000 M E in pebbles is needed, which would be only available for the most metal rich stars. We then follow the evolution of these planets for a few My to compare the resulting population with the observed exoplanet populations. Planet formation in discs with these large amounts of pebbles result in mostly forming gas giants and only very little super-Earths, contradicting observations. This leads to the conclusions that either (i) the observed discs are exceptions, (ii) not all gaps in observed discs are caused by planets or (iii) that we miss some important ingredients in planet formation related to gas accretion and/or planet migration.

Equations

$${\tau }_{\mathrm{f}}=\frac{{\rho }_{\bullet }R}{{\rho }_{\mathrm{g}}{H}_{\mathrm{g}}},$$

$$M_{\mathrm{t}}=\sqrt{\frac{1}{3}}\frac{\left(\eta {\upsilon }_{\mathrm{k}}\right)}{G{\Omega }_{\mathrm{k}}}.$$

$${\dot{M}}_{\mathrm{planet}}=2{\left({\tau }_{\mathrm{f}}/0.1\right)}^{2/3}{r}_{\mathrm{H}}{v}_{\mathrm{H}}{\Sigma }_{\mathrm{peb}({\mathrm{r}}_{\mathrm{p}})}.$$

$${\Sigma }_{\mathrm{peb}}=\sqrt{\frac{2{\dot{M}}_{\mathrm{peb}}{\Sigma }_{\mathrm{g}({\mathrm{r}}_{\mathrm{p}})}}{\sqrt{3}\pi {ϵ}_{\mathrm{p}}{r}_{\mathrm{p}}{\upsilon }_{\mathrm{k}}}},$$

$${\dot{M}}_{\mathrm{peb}}=2S_{\mathrm{peb}}\pi r_{\mathrm{g}}\frac{d{r}_{\mathrm{g}}}{dt}\left(Z_{\mathrm{peb}}{\Sigma }_{\mathrm{g}({\mathrm{r}}_{\mathrm{g}})}\right).$$

$$r_{\mathrm{g}}={\left(\frac{3}{16}\right)}^{\frac{1}{3}}{\left(G{\mathrm{M}}_{\star }\right)}^{\frac{1}{3}}{\left({ϵ}_{\mathrm{D}}{Z}_{\mathrm{peb}}\right)}^{\frac{2}{3}}{t}^{\frac{2}{3}}$$

$$\frac{d{r}_{\mathrm{g}}}{dt}=\frac{2}{3}{\left(\frac{3}{16}\right)}^{\frac{1}{3}}{\left(G{\mathrm{M}}_{\star }\right)}^{\frac{1}{3}}{\left({ϵ}_{\mathrm{D}}{Z}_{\mathrm{peb}}\right)}^{\frac{2}{3}}{t}^{-\frac{1}{3}},$$

$$M_{\mathrm{iso}}$$

$${\dot{M}}_{\mathrm{gas}}$$

$${\dot{M}}_{\mathrm{gas},\mathrm{low}}=0.83{\Omega }_{\mathrm{k}}{\Sigma }_{\mathrm{g}}{H}^2{\left(\frac{r_{\mathrm{H}}}{H}\right)}^{\frac{9}{2}}$$

$${\dot{M}}_{\mathrm{gas},\mathrm{high}}=0.14{\Omega }_{\mathrm{k}}{\Sigma }_{\mathrm{g}}{H}^2.$$

$${\Gamma }_{\mathrm{tot}}={\Gamma }_{\mathrm{L}}+{\Gamma }_{\mathrm{C}}.$$

$${\tau }_{\mathrm{II}}={\tau }_{\upsilon }\times max\left(1,\frac{M_{\mathrm{p}}}{4\pi {\Sigma }_{\mathrm{g}}{r}_{\mathrm{p}}^2}\right).$$

$$f=10^{-0.62\pm 0.22}{\left(\frac{M_{\mathrm{P}}}{M_{\mathrm{Sat}}}\right)}^{-0.73\pm 0.17}.$$

Full text

††pubyear: 2019††pagerange: Are the observed gaps in protoplanetary discs caused by growing planets?–Are the observed gaps in protoplanetary discs caused by growing planets?

1 Introduction

Recent observations of dust grains in protoplanetary discs have revealed a large level of detailed substructures (ALMA Partnership et al., 2015; Andrews et al., 2016). However, up to now it is unclear what mechanism causes the gaps in discs despite of the ever increasing observational evidence for substructures in the discs (Walsh et al., 2014; van der Marel et al., 2015; ALMA Partnership et al., 2015; Isella et al., 2016; Fedele et al., 2017; Huang et al., 2018; van der Marel et al., 2019). Proposed mechanisms to explain gaps in protoplanetary discs include dust growth in condensation zones or snow lines (Zhang et al., 2015), zonal flows (Flock et al., 2015), self-induced dust pile-ups (Gonzalez et al., 2017), assemblage sintering (Okuzumi et al., 2016), large scale instabilities (Lorén-Aguilar & Bate, 2016), also notably the secular gravitational instabilities (Takahashi & Inutsuka, 2016) and, gap opening via clearing by embedded planets (Pinilla et al., 2012b, a; Pinilla et al., 2015).

In this work we will follow the assumption that all gaps in discs are caused by growing planets. Low mass planets only perturb the disc slightly and do not open deep gaps in the disc. However, even low mass planets can perturb the disc large enough to generate a small pressure bump exterior to their orbits. This pressure bump results in reversing the gas-drag such that small dust grains and pebbles are halted in their inward drift (Whipple, 1972; Paardekooper & Mellema, 2006; Pinilla et al., 2012a; Morbidelli & Nesvorny, 2012; Lambrechts et al., 2014; Bitsch et al., 2018b). This results in an accumulation of particles exterior to the planet which can explain the formation of rings and gaps in the disc observations (Dullemond et al., 2018). Additionally, as the pebbles are trapped exterior to the planet they can not be accreted onto the planet any more. The planet has thus reached pebble isolation mass and can start to accrete gas (e.g. Morbidelli & Nesvorny 2012; Lambrechts et al. 2014; Bitsch et al. 2018b). As the pebbles pile up exterior to the planet in rings, a planet that causes such a change in the disc structure must thus be located in a gap of the dust distribution.

Planet population synthesis generates a synthetic planet population by randomizing the initial conditions (e.g. metallicity, planetary embryo starting position, disc lifetime) relevant for planet formation. Previous planet population syntheses considered solid growth via planetesimal accretion (Ida & Lin, 2004; Alibert et al., 2009; Mordasini et al., 2009). In these syntheses models, the solid component of the planets is assumed to grow by accreting planetesimals, although numerical simulations show the inefficiency of this process (Tanaka & Ida, 1999; Levison et al., 2010), especially for large planetesimals of several 10km in size. The accretion of small mm-cm sized pebbles, on the other hand, can significantly enhance the growth of planets (Ormel & Klahr, 2010; Lambrechts & Johansen, 2012; Johansen & Lacerda, 2010; Lambrechts & Johansen, 2014; Bitsch et al., 2015b; Levison et al., 2015). However, most of these simulations have only looked at specific problems like the formation of super-Earths (Lambrechts et al., 2019; Izidoro et al., 2019; Bitsch et al., 2019a) or just gas giants (Levison et al., 2015; Chambers, 2016; Bitsch et al., 2019b).

Recently, pebble-based planet population synthesis approaches appeared in the literature (Bitsch & Johansen, 2017; Ndugu et al., 2018; Chambers, 2018; Brügger et al., 2018; Johansen et al., 2019). These recent models allowed planetary seeds to accrete pebbles until pebble isolation mass, when gas accretion can start. These models also include type-I and type-II migration.

The detailed statistics of observed gaps in discs from the Disk Substructures at High Angular Resolution Project (hereafter, DSHARP) survey (Huang et al., 2018) provide us with an excellent statistics of the number and location of rings and gaps in protoplanetary discs. We adopt the pebble-based planet population synthesis approach (Bitsch & Johansen, 2017; Ndugu et al., 2018) to constrain core accretion models to the observed gap statistics (Huang et al., 2018), radial velocity observations of exoplanets (Johnson et al., 2010) and to the microlensing observations of exoplanets (Cassan et al., 2012).

This paper is organized as follows: In subsection 2.1 we present the used disc model, explain the planet formation model, and discuss the initial conditions of our simulations in subsection 2.2 and 2.3, respectively. In section 3, we match the synthetic gap opening position data to the observed gap statistics of the DSHARP disc survey. In section 4 we show the match of the synthesized planet population at the end of the disc lifetime with the radial velocity and microlensing data. Section 5 summarizes the main finding in the paper and put them into future investigative perspectives.

2 Theoretical Models

2.1 Disc Model

Our disc model is derived from the detailed radiative 2D hydrodynamical simulations of Bitsch et al. (2015a) which feature viscous and stellar heating. The Bitsch et al. (2015a) disc model is based on radially constant accretion rates, where the stellar accretion rate decreases in time for a fixed $\alpha$-viscosity parameter of $0.0054$. Reduction of accretion rates translate into reduced gas surface densities, which in turn reduce viscous heating and thus the disc’s temperature. The same disc model was used in the pebble-based planet formation simulations (Bitsch et al., 2015b; Bitsch & Johansen, 2016; Ndugu et al., 2018) and in N-body simulations (Izidoro et al., 2017, 2019; Bitsch et al., 2019a). The disc model used in our work has a lifetime which can span all the way up to 10 My as expected by observation (Hartmann et al., 1998; Mamajek, 2009).

However, there is now growing evidence that the bulk of discs are inviscid because the magneto-rotational instability is quenched at all scale-heights in discs (Turner et al., 2014). The accretion is thus thought to be driven by disc winds launched from the surface (Bai, 2013, 2014; Gressel et al., 2015; Bai et al., 2016; Suzuki et al., 2016). However, there is no clear census yet on how to implement these new features into planet population synthesis models. Recent simulations, for example the planet population synthesis simulations by Ida et al. (2018) and Johansen et al. (2019) adopted similar modifications to their disc models. Following this approach we use $\alpha=0.0054$ to calculate the disc structure, but use a different viscosity parameter to mimic planet migration and the pebble scale height, $\alpha_{\rm mig}$ which can be as low as $0.0001$. Throughout this paper, we note that $\alpha$ is the viscosity that sets the disc’s structure (temperature, aspect ratio and gas surface density) and the gas accretion rate onto the planet, while $\alpha_{\rm mig}$ is the bulk viscosity that is reflected in the growth parameters (migration and pebble accretion) of the planetary cores.

2.2 Planet formation and migration models

Our planet formation model requires accretion of pebbles onto a core, which eventually triggers runaway gas accretion once the pebble isolation mass is reached and the gaseous envelope contracted. Pebbles with the Stokes number

[TABLE]

are accreted onto cores that have reached the pebble transition mass (Lambrechts & Johansen, 2012, 2014)

[TABLE]

Here, $\rho_{\rm{\bullet}}$ is the pebble density, $R$ is the pebble radius, $\rho_{\rm g}$ is the gas density, $H_{\rm g}$ is the gas scale height, ${\eta}$ is the pressure support, $G$ is the gravitational constant, $\upsilon_{\rm k}$ is the Keplerian speed and $\Omega_{\rm k}$ the Keplerian frequency. Initially, the core accretes pebbles in a 3D fashion because the planetary Hill radius is smaller than the pebble scale height (Lambrechts & Johansen, 2012; Morbidelli et al., 2015). As the core agglomerates more pebbles, it becomes more massive and accretes in a 2D manner because $r_{\rm H}>H_{\rm peb}$ (Morbidelli et al., 2015) with a growth rate of

[TABLE]

The pebble scale height is given as $H_{\rm peb}=\sqrt{\frac{\rm\alpha_{mig}}{\tau_{\rm f}}}H_{\rm g}.$ The pebble surface density at the planets position is given by

[TABLE]

where $r_{\rm p}$ denotes the planet’s semi-major axis and the pebble flux is given by

[TABLE]

Here $Z_{\rm peb}$ is the pebble metallicity, $S_{\rm peb}$ is the pebble flux scaling factor (see below) and $r_{\rm g}$ the pebble production line given by

[TABLE]

and

[TABLE]

where $\rm M_{\rm\star}$ is the stellar mass, $\epsilon_{\rm p}$ and $\epsilon_{\rm D}$ are the pebble and dust coagulation efficiency with values 0.5 and 0.05, respectively. This is the standard pebble accretion paradigm (Lambrechts & Johansen, 2014) , where $\epsilon_{\rm p}$ and $\epsilon_{\rm D}$ are values corresponding to detailed coagulation simulations (e.g. Birnstiel et al. 2011). Within the framework of Lambrechts & Johansen (2014), changes of the accretion rate could only happen due to changes in coagulation efficiency. The here used pebble growth model was described in detail in Bitsch et al. (2015b). However, Bitsch et al. (2015b) included an unphysical increase in the pebble formation process related to the evolution of the protoplanetary disc structure. This has been discussed in detail in Bitsch et al. (2018a). We thus here just scale the pebble flux by $S_{\rm peb}$ to allow faster planetary growth, similar to the approach used by Izidoro et al. (2019) and Bitsch et al. (2019b).

When the planet reaches pebble isolation mass, it accelerates the gas outside its orbit to super-Keplerian velocities, which generate pressure bumps halting pebble accretion, trapping inward moving pebbles and allows the planet to transition into a stage of gas accretion (Lambrechts et al., 2014). Through out our simulations, the new pebble isolation mass relation of Bitsch et al. (2018b) given by

[TABLE]

is used. $P_{\rm grad}$ corresponds to the radial pressure gradient in the unperturbed protoplanetary disc. This clearly shows the importance of the disc structure and viscosities for planetary growth and migration simulations.

At the pebble isolation mass, the gas envelope contracts over a long time while accreting some gas as long as $\rm M_{\rm env}<\rm M_{\rm core}$ at a rate given by (Piso & Youdin, 2014; Bitsch et al., 2015a)

[TABLE]

where $f$ is a scaling factor given by 0.2 (Piso & Youdin, 2014) and $\kappa_{\rm env}=0.05{\rm cm^{2}/g}$ is the planet’s envelope opacity very similar to the study by Movshovitz & Podolak (2008). $\rho_{\rm C}$ is the core density and we assumed $\rho_{\rm C}$ = 5.5 $\rm{g/cm^{3}}$. The planet continues to contract its envelope until $M_{\rm c}=M_{\rm env}$. When $M_{\rm c}<M_{\rm env}$, rapid gas accretion is triggered and the gas accretion rate is given by the minimum of (Machida et al., 2010)

[TABLE]

and

[TABLE]

We note that these gas accretion rates might be too generous and allow the formation of too many gas giants (e.g. Brügger et al. 2018; Bitsch et al. 2019b), however, until the planets reach pebble isolation mass gas accretion can be neglected. We will investigate different gas accretion rates in an upcoming publication and first tests imply that different gas accretion rates do not influence our results significantly.

Low mass planets radially drift inwards in the type-I fashion where their migration is modeled following the torque formula from Paardekooper et al. (2011), where the total torque, $\Gamma_{\rm tot}$ acting on the planet is given by

[TABLE]

$\Gamma_{\rm L}$ and $\Gamma_{\rm C}$ are the Lindblad and corotation torques, respectively. The Lindblad and corotation torques strongly dependent on the local radial gradients of gas surface density, $\Sigma_{\rm g}\propto r^{-\lambda}$, temperature $T\propto r^{-\beta}$, and entropy $S\propto r^{-\varepsilon}$, with $\varepsilon={\beta}+\left(\gamma-1.0\right)\lambda$ and $\gamma=1.4$ is the adiabatic index. When the planet grows massive enough, it carves gap in the disc (Crida et al., 2006) and eventually switches to the slower type - II migration (Baruteau et al., 2014):

[TABLE]

Here, $\tau_{\rm\upsilon}$ is accretion time scale of the protoplanetary disc. Combining the disc model, planet formation and migration models, we synthesize population of planets by sampling the initial conditions. For this, we slove equation 1-13 numerically. This is described in detail in the following subsection.

2.3 Initial conditions

Planet populations synthesis requires a randomization of important input parameters of discs and planet formation as these parameters are not a-priori known. These important input conditions are: (i) Disc metallicity and lifetime, which set the disc structure ($\Sigma_{\rm g}$, H/r) and thus influences the growth and migration speeds. The disc structure and its evolution are described in detail in Bitsch et al. (2015a). The gas surface density follows in the outer parts of the disc a power law with radius of -14/15 and the aspect ratio follows a 2/7 power law in radius. (ii) Implantation time and position which sets the growth and migration trajectories of cores. Our input parameters are sampled in similar manners as in Ndugu et al. (2018) but with the exception of subjecting the disc populations to exponentially decaying lifetimes as predicted by observations (e.g. Hartmann et al. 1998; Mamajek 2009). The disc lifetime is exponentially decaying with a half lifetime of 2.5 My, where we fixed the minimal and maximal lifetime of the disc to 0.1 and 10 My, respectively. Throughout the simulation, our exponentially-decaying surface density prescription is not derived from any underlying physical model but purely empirically derived from observations. The initial starting locations follow a linear distribution in semi-major axis. We show results for a logarithmic starting configuration in the appendix, but this does not affect our overall conclusions. Our seeds are placed up to 160 au. The starting time of the planetary seeds corresponds to the time when the pebble production line (Eq.6) reaches the starting position plus 100 orbital times, as this corresponds roughly to the timescale of the streaming instability to form planetesimals (e.g. Johansen et al. 2015a). We then start our planetary embryos at the pebble transition mass. This approach might not be realisitic, as planetesimals might take a long time to coagulate into planetary embryos at large distances, it thus gives us the maximal possible growth time by pebble accretion for embryos in the outer disc.

3 Match to the DSHARP campaign

The DSHARP campaign limits the alpha-viscosity parameter within the range of $10^{-6}$ to $10^{-2}$. Our choices of alpha are thus within the limits of the observations. A higher viscosity results in higher pebble isolation masses, resulting in planets reaching the pebble isolation mass only in the inner regions of the disc. Lower viscosities, on the other hand, reduce the pebble isolation mass and allow planets to reach the pebble isolation mass also in the outer parts of the disc. We generated a population of cores that just reach pebble isolation mass following the initial conditions described in subsection 2.3. Our approach is the simplistic one-embryo-per-disc. We are aware that single embryo-per disc planet population synthesis models might present an over estimate in the growth efficiency of cores due to the lack of competition for materials and the neglection of eccentricity effects which might affect the accretional efficiency (Johansen et al., 2015b; Liu & Ormel, 2018; Ormel & Liu, 2018) or the migration rates (Bitsch & Kley, 2010; Cossou et al., 2013). Our result is therefore the upper limit representation of the accretional efficiency of single planets for materials in the discs.

We make the assumption that the gaps in the protoplanetary discs are opened at positions where the planetary embryos have reached pebble isolation mass. To compare our simulations to the observations, we stop planet growth and migration once the planet reaches pebble isolation mass and track at what time the planet reached pebble isolation mass. We follow the evolution of 3000 planetary embryos for each set of simulations, where the positions of the embryos that have reached pebble isolation mass is compared to the position of the gaps seen in observations. The observed gap positions are extracted from Fig. 2 of Huang et al. (2018) which is normalised and compared to the normalised occurrence of the simulated gap positions as displayed in Fig.1.

From Fig.1 we see that matching the DSHARP observed gaps requires a very high (and probably unrealistic) pebble flux scaling factor, low alpha for migration and particle scale height. Otherwise the observed gaps at the outskirts of the disc can not be reproduced. For low pebble fluxes and larger alpha parameters, the outer gaps can not be matched, because growth is too slow in the outer disc. Low alpha means slower migration and very thin particle scale height, therefore cores radially migrate slower to the inner discs and transit earlier into the faster 2D pebble accretion fashion and can therefore reach the high pebble isolation mass in the outskirts of the protoplanetary disc earlier. At high alpha specifications, even higher pebble flux scaling factors are required to match the observations of the protoplanetary discs.

For a planet to reach the pebble isolation mass in the $S_{\rm peb}=10$ case in the outer disc (red curve in Fig.1) for a disc with [Z/H]=0.5 (see Fig.2), a total of about 2000 Earth masses of pebbles have to pass the planet. This amount of pebbles is about a factor of 10-20 higher than the pebbles trapped in the observed rings from DSHARP survey (Dullemond et al., 2018). In our model, [Z/H]=0.0 corresponds to a dust-to-gas ratio of 0.015, indicating that a disc with [Z/H]=0.5 would feature a dust-to-gas ratio of 0.047. This implies that the gas mass of a disc with [Z/H]=0.5 containing 2000 Earth masses in pebbles corresponds to $\approx$ 0.14 Solar masses, which is on the higher mass ranges of protoplanetary discs. We also note that in order to match the gaps seen by the DSHARP survey our simulations require a pebble production line that extends all the way to around 300 au. This is clearly beyond the the outer boundaries of the disc according to observations.

Therefore given the limitations of our model, reproducing the observed gaps in the DSHARP survey requires low alpha values, a very extended pebble production line and a pebble scaling factor of around $S_{\rm peb}=10$ resulting in very large amounts of pebbles needed to form planets at the position of the gaps.

Setting the pebble scale height, $H_{\rm peb}$ lower by a factor of $50$ but keeping $\alpha_{\rm mig}=0.0054$, we investigate how the match of our simulations to the gaps seen by the DSHARP survey changes. Lowering the pebble scale height helps the growth, but the gaps at the very outskirts of the discs are not matched. This is because although the transition into effective 2D peble accretion is faster, the seeds drift inwards before generating pressure bumps.

We show the time at which planets reach pebble isolation mass as a function of the host star metallicity in Fig.2. This figure shows a much wider spread in time of gap opening (0.01:10 My) than the specifications in the DSHARP campaign (1:2 My).The planets in the inner regions of the disc grow very fast, because there pebbles are produced immediately and the accretion rates are very high, enhancing the growth. Fig.2 tells us that making gaps in the outskirts of the discs require large amount of materials and additionally a large metallicity to achieve efficient planet growth in the outer regions of the protoplanetary discs. Lower metallicity implies that planet formation takes longer than the observed disc ages (1-2 My). This would imply that there is a correlation between host star metallicity and distance of giant planets. The more metal rich a host star is, the farther away giant planets could be located.

Following this time specifications for planetary embryos reaching isolation mass, full simulation runs were executed until the disc reaches its final lifetime of up to 10 My. These simulations were obtained for the model with $\alpha_{\rm mig}=0.0001$ and $S_{\rm peb}=1$ and for the model that gives good match to the gap occurrence rates of the DSHARP survey (model with $\alpha_{\rm mig}=0.0001$ and $S_{\rm peb}=10$). These models are indicated by the black and the red curve in Fig.1. We show the results of these planet population synthesis models in Fig.3 as function of the disc lifetime and in Fig.4 as a function of host star metallicity. These plots clearly show the importance of disc lifetimes and metallicity in the synthesised planet population simulations. These simulations clearly predict a large population of giant planets that could be seen by direct imaging surveys, for example, the HR8799 system, but rarely by the RV method. The observed gas giants in reality are mainly on eccentric orbits, therefore they most likely underwent gravitational scattering (Jurić & Tremaine, 2008; Sotiriadis et al., 2017). This means the large population of synthesised giant planets in the direct imaging field of view could be scattered by gravitational interaction which might refill the dearth of the planets in the RV field of view should N-body interaction between the protoplanets been incorporated. However, the planets in our simulation seem to be too massive compared to the observations. The comparison of subsets of the synthesised population of planets to the observational limits of the RV and the microlensing technique is shown in section 4.

4 Match to radial velocity and microlensing observation

Fig.5 shows the comparison of the synthesised gas giant sub-population to the RV measurements of Johnson et al. (2010) in a manner similar to the approach of Ndugu et al. (2018). The plot shows that we do not produce enough gas giants compared to the RV observations for both models, even though the low alpha and high $S_{\rm peb}$ model shows a very good match to the gap statistics. One clear reason for this mismatch is the lower migration speeds which deprives the gas giants from reaching the RV field of view before the disc dissipates. We show simulations with a logarithmic starting configuration of the embryos in the appendix, which essentially recover our here discussed results. To broaden the scope of our comparison, we match the synthesis data to the microlensing survey of Cassan et al. (2012). The microlensing survey we used considered planets inside the purple rectangle in Fig. 3 and Fig.4 (see Cassan et al. (2012) for the detailed limits). From the survey in Cassan et al. (2012), the planet frequency (for masses and distances as indicated in Fig.3 and Fig.4) scales as

[TABLE]

The fit assume a Saturnian pivot-mass (${M_{\rm Sat}}$) of 95 Earth masses. From the microlensing data, low mass planets are expected to be more common than giant planets at distances of 0.5-10 au. The fits of their survey is represented by the gray band of Fig.6. From the comparison to the microlensing survey, both models overestimate the gas giant occurrence rate and underestimated the low mass planet occurrence rates even for the model which reproduced the observed gap occurrence rates in the DSHARP campaign very well. This presents yet another direct conflict with the match to the gap occurrence rates. Testing the behavior of this conflict for logarithmic starting position distribution (see, appendix), we found the same trends and mismatch between the models in matching the observed gap positions, RV field of view and the microlensing survey. This leaves us with the following summary of our simulations: If all the gaps and rings in the DSHARP survey are caused by growing planets formed in the core accretion scenario then extremely high (even unrealistically high) amounts of pebbles are needed. This implies that either not all gaps are caused by forming planets or that we are still missing some ingredients to the core accretion scenario. For example, naturally occurring pressure bumps could help in assisting core accretion in the outer parts of the disc (e.g. Coleman & Nelson 2016).

5 Discussions

Performing planet population synthesis from more physical photoevaporative disc models does not significantly differ from our approach. For example, in the photoevaporative disc model by Monsch et al. (2019), the hole of photoevaporation is around 1-2 au, which is interior to the region where the planets resides at the end of the gas-disc lifetime in our simulations with low viscosity. Ideally this mechanism would help us to keep planets in the RV bin, because planet migrates less in low alpha specifications. However, we basically have no hot Jupiters that could be trapped further out by changes in the disc structure due to photoevaporation at 1-2 au. This would only play a role if the viscosity and thus the type-II migration rate is higher. At higher migration rates, it is not possible from our simulation for cores to open gaps exterior to 100 au. Therefore this is in conflict with the DSHARP survey, where gaps are also observed in the outer discs.

Our outer disc is basically not different from other detailed radiative disc models as it is determined by stellar irradiation. For example, in Baillié et al. (2015), the outer disc structure looks very similar to our model, with H/r following a 2/7 power law with radius (Chiang & Goldreich, 1997). As the outer disc structure is similar in the disc models, because it is determined by stellar irradiation, planet formation should follow equivalent pathways.

The only difference appears in the location of snow regions, where several snow regions were reported within 0.1-1 au in Baillié et al. (2015). These differences arise from the different opacity model used in Baillié et al. (2015), who used the opacities from Semenov et al. (2003), which feature multiple condensation fronts for different chemical species, which are not included in our model.

However, the inner 1 au region is only reached in our model by gas giants, if the viscosity and thus the migration speed is high, which fails to explain the observed gap structure. In the case of low viscosity, outward migration around these ice lines would cease to exist, as the corotation torque saturates at low viscosity (e.g. Paardekooper et al. (2011)). Planets would thus also migrate inwards as in our model. This implies that even the difference in the inner structure of the disc does not have a big influence on our results.

On the other hand, if steeper and simpler power laws for the disc structure were used (e.g. an MMSN like power law disc (Weidenschilling, 2006) ), planet migration would be fast and inwards in the inner regions of the disc as noted in Bitsch et al. (2015b). However, the steeper slope of the surface density profile in the MMSN compared to our disc would greatly reduce planet formation in the outer disc, making it hard to form planets in the outer disc with this steep profile.

6 Summary

In this paper, we took the pebble-based planet population synthesis approach of Ndugu et al. (2018) to simulate statistics of gaps observed in protoplanetary discs (Huang et al., 2018). The planet population model is in the one embryo-per-disc fashion where single planetary embryos evolve to pebble isolation mass where we compare our synthesized planet population with the observations of gaps in discs. The comparison of the synthesized data to the observed gap statistics (Fig.1 and Fig.2) asserts the following:

• (i)

The rings and gaps in the observed disc distributions can be caused by planets, however, if these planets form by the core accretion scenario very high (maybe even unrealistically high) amounts of pebbles are needed, especially to explain the formation of the rings exterior to 100 au.

• (ii)

The timescale to form planets that could correspond to the outermost gaps and rings in the observations is very close or even longer than the inferred lifetime of the observed discs. The outermost rings can only be match with very high pebble amounts and elevated metalicities, which could lead to the prediction that giant planets far way from the central star should be more frequent around high metallicity stars.

• (iii)

For our model to match the gap statistics observed by the DSHARP survey, a pebble production line at around 300 au would be needed. This is larger than the sizes of the observed discs. In order to match the pebble accretion scenario under the assumption that all gaps are caused by planets, the pebbles from the outer disc would have needed to drift inwards at the time of the disc observations.

In the second stage of our simulations, planet evolution is extended until the protoplanet migrates to the inner edge of the disc (0.004 au) or the end of the disc lifetime is reached. The synthetic planet population is then compared with a radial velocity observational survey (Johnson et al., 2010) and a microlensing survey (Cassan et al., 2012). We clearly see from these comparisons that planets formed in the environments that match the observed gaps grow basically all to giant planets and we have no small planets left (Fig.6). This is not in agreement with observations which show that close in super-Earths (Mayor et al., 2011; Mulders et al., 2018) and ice giants around the water ice line (Cassan et al., 2012; Suzuki et al., 2018) should be the most common type of planets. Our simulations thus predict, if the protoplanetary discs seen by the observations are typical discs in respect to their size and heavy element content, a large population of giant planets, which could be seen by direct imaging surveys. However, this population has not been found up to date. This leaves us with the following, not mutually excluding, conclusions:

• (i)

Not all rings can be caused by planets, because of the contradictions to the observations.

• (ii)

Planet formation simulations are missing some important ingredients regarding gas accretion (see also Nayakshin et al. (2019)) and planet migration.

• (iii)

The discs of the DSHARP survey are exceptions due to their observational bias (Andrews et al., 2018) which do not represent the norm of protoplanetary discs and planet formation.

The survey by Long et al. (2018) in the Taurus star-forming region and the survey by Cieza et al. (2019) in the Ophiuchus star-forming region, on the other hand, found a significant population of very compact discs (outer disc radii ¡ 30 au). These discs are not pictured within the DSHARP sample, due to the selection biases in the DSHARP campaign (Andrews et al., 2018). If these discs are indeed exceptions and only 1 % of all the disc produce these wide orbit gaps, our simulations would only apply to 1 % of the stars and could thus represent the observed number of giant planets seen by direct imaging campaigns.

We thus urge to investigate protoplanetary discs in more detail to get reliable statistics about their sizes and mass as function of their age to constrain the immediate steps of planet formation models.

Acknowledgments

N.N & E.J acknowledge financial support from the International Science Program (ISP). B.B, thanks the European Research Council (ERC Starting Grant 757448-PAMDORA) for their financial support.

Appendix A Influence of the starting position of the planetary cores

Unfortunately, no observational constraints exist about the initial distribution of planetary embryos in protoplanetary discs. To test how the initial starting position distribution of seeds affect our results, we present in this section the findings for simulations with a logarithmic distribution of the starting positions. Apart from the starting position distribution, we maintained the same parameter configurations as in our previous simulations. We use $\alpha_{\rm mig}=0.0001$ but with $S_{\rm peb}=1$ and $S_{\rm peb}=10$, respectively. We show the results of these simulations in Fig.A1, A2 and A3. We see that the starting configurations of the planetary embryos can have important consequences for the final configuration of the planetary systems. The logarithmic starting configuration shifts, by construction, the synthetic gap occurrence rates to the inner discs. This leads to a better match in the comparison to the RV detection of giant planets, as more planets are formed in the inner disc which can end up within the RV detection limits, compared to a linear starting configuration where most planets end up in the outer disc (Fig.3 and Fig.4). However, as stated in the main paper, scattering events if multiple planets grow in the disc could move planets from the outer population into the detection bin of RV surveys.

A comparison with the microlensing observations reveals a very similar trend as for the linear starting configurations of the planetary seeds. Our simulations over predict the amount of giant planets and under predict the number of super-Earths and ice giants. Both, the linear and logarithmic starting configuration, thus fail at reproducing the observational constraints.

The differences between the two different starting configurations are very small and do not change our comparison to observations. We thus conclude that the starting configuration does not play a major role in our assessment, but the mismatch to observations is rather explained by the conclusions made in Section 5.