Can Uranus and Neptune form concurrently via pebble, gas, and planetesimal accretion?

ABSTRACT

The origin of Uranus and Neptune has long been challenging to explain, due to the large orbital distances from the Sun. After a planetary embryo has been formed, the main accretion processes are likely pebble, gas, and planetesimal accretion. Previous studies of Uranus and Neptune formation typically do not consider all three processes; and furthermore, do not investigate how the formation of the outer planet impacts the inner planet. In this paper, we study the concurrent formation of Uranus and Neptune via pebble, gas, and planetesimal accretion. We use a dust-evolution model to predict the size and mass flux of pebbles, and derive our own fit for gas accretion. We do not include migration, but consider a wide range of formation locations between 12 and |$40\, \textrm {au}$|⁠. If the planetary embryos form at the same time and with the same mass, our formation model with an evolving dust population is unable to produce Uranus and Neptune analogues. This is because the mass difference between the planets and the H–He mass fractions become too high. However, if the outer planetary embryo forms earlier and/or more massive than the inner embryo, the two planets do form in a few instances when the disc is metal-rich and dissipates after a few Myr. Furthermore, our study suggests that in situ formation is rather unlikely. Nevertheless, giant impacts and/or migration could potentially aid in the formation, and future studies including these processes could bring us one step closer to understanding how Uranus and Neptune formed.

planets and satellites: formation, planets and satellites: gaseous planets, planets and satellites: general

1 INTRODUCTION

Uranus has a mass of |$14.5\, \textrm {M}_{{\oplus }}$| and a semimajor axis of |$19.1\, \textrm {au}$|⁠, placing it around |$10\, \textrm {au}$| beyond the orbit of Saturn. The second ice giant Neptune has a mass of |$17.1\, \textrm {M}_{{\oplus }}$| and a semimajor axis of |$30.0\, \textrm {au}$|⁠, placing it at the inner edge of the Kuiper Belt and around |$10\, \textrm {au}$| beyond Uranus’ orbit. Unlike the Solar system gas giants, Uranus and Neptune are primarily composed of heavy elements¹ with a hydrogen–helium (H–He)-fraction |$\lesssim$|20 per cent (see Helled & Bodenheimer 2014 and references therein).

The formation of Uranus and Neptune has been challenging for the planet formation models since several decades. Safronov (1969) demonstrated that the time-scale for core formation via planetesimal accretion exceeds the lifetime of the disc at the orbits of Uranus and Neptune. This time-scale issue can be alleviated by instead considering the accretion of small mm–cm sized pebbles (e.g. Lambrechts, Johansen & Morbidelli 2014; Venturini & Helled 2017). However, when the accretion of gas is properly considered, the challenge reverses and the difficulty is to prevent the planets from undergoing runaway gas accretion and becoming gas giants (Helled & Bodenheimer 2014, Brouwers et al. 2021). This problem can be circumvented by removing the gas disc at the exact right time, which is known as the fine-tuning problem. Furthermore, both pebble accretion and planetesimal accretion have an efficiency which is generally decreasing with orbital distance, making it hard to explain why Neptune has a mass that is higher than Uranus.

Some of the above concerns can be alleviated if the planets have not formed in situ. In fact, in order to explain multiple properties of the Kuiper Belt and the Oort cloud, it has been suggested that Neptune must have migrated outwards to reach its current location (e.g. Malhotra 1995; Hahn & Malhotra 1999; Nesvorny 2015). In Nesvorny (2015), it was found that Neptune must have been located interior of |$25\, \textrm {au}$| and slowly migrated outwards over a time-scale longer than |$10\, \textrm {Myr}$|⁠. A formation location closer to the Sun implies a more efficient core accretion, making it possible to form the planets over a shorter time-scale. In the Nice model, Uranus and Neptune are formed at distances of |$\sim 12-20\, \textrm {au}$|⁠, before an instability occurs which causes them to move outwards and eventually end up at their current orbital locations (Tsiganis et al. 2005). During this instability, Uranus and Neptune also switched places in some of their simulations, a mechanism which can explain why Neptune has a mass higher than that of Uranus.

The aforementioned studies concerned late migration that occurred after the dispersal of the gas disc; however, planets can also undergo migration during the disc lifetime due to interactions with the surrounding gas disc (see Paardekooper et al. 2022 for a recent review on planet migration). The direction and magnitude of this migration depends on both planet and disc properties. Interactions with dust in the disc can further modify the migration properties (Guilera et al. 2023). There is no evidence that suggests that Neptune could not have had an early period of migration driven by interactions with the protoplanetary disc, prior to the late outward migration which placed it at its current location. The existence of the classical Kuiper belt beyond |$40\, \textrm {au}$| does put some constraints on the initial formation locations, since the planets could not have formed in or passed through that region without disturbing the planetesimal belt.

Taken together, the formation of Uranus and Neptune is still very unconstrained and suffers numerous challenges. Valletta & Helled (2022) used a pebble accretion model coupled with a realistic gas accretion model to study the possibility of forming Uranus and Neptune in situ, and found that they could indeed obtain planets with the right masses and H–He mass fractions. In this work, we develop this formation model further by considering the accretion of pebbles, gas, and planetesimals in an evolving disc. We use a state-of-the art dust evolution model and also account for the blocking of pebbles towards Uranus by accretion onto Neptune. Furthermore, we consider a wide range of formation locations and disc parameters. The aim of our study is to investigate whether both planets could form concurrently without invoking additional growth mechanisms such as giant impacts.

We present our models for disc evolution and planet growth and describe our simulation set-up in Section 2. The results for in situ formation are shown in Section 3, and in Section 4, we show the results from when the formation locations are varied. In Section 5, we discuss potential mechanisms that can change our results compared to the previous sections, and finally our main findings are summarized in Section 6.

2 OUR MODEL

We consider growth via pebble, gas, and planetesimal accretion in an evolving one-dimensional global disc model. We use two methods for determining the pebble flux: in the constant model the pebble flux is proportional to the gas flux through the disc, while in the evolving model pebble flux is calculated using a dust evolution code. We derive a fit for gas accretion based on mesa simulations (Paxton et al. 2011, 2013, 2015, 2018, 2019). Planetesimal accretion is modelled using a semi-analytic model and assuming that the initial planetesimal surface density is proportional to the gas surface density. The effect of planetary migration is not included; however, we do test a large range of formation locations.

2.1 Disc model

The evolution of the disc’s surface density is modelled using the analytic solution for an unperturbed thin accretion disc from Lynden-Bell & Pringle (1974),

$$\begin{eqnarray} \Sigma _\mathrm{\rm disc}(t) = \frac{\dot{M}_{\rm disc}}{3 \pi \nu _\mathrm{out} (r/r_\mathrm{out})^\gamma } \exp \left[{-\frac{(r/r_\mathrm{out})^{2-\gamma }}{T_\mathrm{out}}}\right], \end{eqnarray}$$

(1)

where T_out = t/t_s + 1 and |$t_\mathrm{ s} = 1/(3 [2-\gamma ]^2) \times r_\mathrm{out}^2/\nu _\mathrm{out}$|⁠. In the above expressions t is the time, |$\dot{M}_{\rm disc}$| is the disc accretion rate, ν_out is the kinematic viscosity at semimajor axis r = r_out, and γ is radial viscosity gradient. We adopt |$r_{\rm out}=50\, \textrm {au}$| and γ = 15/14. The kinematic viscosity is given by,

$$\begin{eqnarray} \nu = \alpha \Omega _\mathrm{ K} H_\mathrm{ g}^2, \end{eqnarray}$$

(2)

(Shakura & Sunyaev 1973), where α is the viscosity parameter, Ω_K is the Keplerian angular velocity, and H = c_s/Ω_K is the scale height of the gas disc. We use α = 0.005 throughout this study. The sound speed is calculated as c_s = (k_BT/[μm_H])^1/2 where K_B is the Boltzmann constant, μ is the mean molecular weight, m_H is the mass of the hydrogen atom, and |$T=150\, \textrm {K} (r/\textrm {au})^{-3/7}$| is the mid-plane temperature of the disc (Chiang & Goldreich 1997). The value of μ is chosen to be 2.34.

The evolution of the disc’s accretion rate is given by,

$$\begin{eqnarray} \dot{M}_{\rm disc} (t) = \dot{M}_{0,\rm disc} T_\mathrm{out}^{-(5/2-\gamma)/(2-\gamma))} \left[1-\left(\frac{t}{t_{\rm disc}}\right)^{3/2}\right], \end{eqnarray}$$

(3)

where the last term is implemented to mimic gradual disc dissipation until time t = t_disc (Ruden 2004). The radial velocity of the disc gas at semimajor axis r can be obtained from the continuity requirement using the expression,

$$\begin{eqnarray} v_{\rm R,disc} = -\frac{\dot{M}_{\rm disc}}{2\pi r \Sigma _{\rm disc}}. \end{eqnarray}$$

(4)

2.2 Pebble accretion model

Growth via pebble accretion is modelled using the exact monodisperse accretion rate from Lyra et al. (2023),

$$\begin{eqnarray} \dot{M}_{\rm pe} = \pi R_{\rm acc}^2 \rho _{\rm pe} \delta v \exp {(-\xi)}[I_0(\xi) + I_1(\xi)], \end{eqnarray}$$

(5)

where

$$\begin{eqnarray} \xi = \left(\frac{R_{\rm acc}}{2H_{\rm pe}}\right)^2. \end{eqnarray}$$

(6)

In the above equations, I_X are modified Bessel functions of the first kind of real order, R_acc is the accretion radius, ρ_pe is the density of pebbles in the mid-plane, δv is the approach speed, and H_pe is the pebble scale height.

We assume that all pebble accretion occur in the Hill regime, and thus the approach speed and accretion radius can be approximated as,

$$\begin{eqnarray} \delta v = \Omega _\mathrm{ K} R_{\rm acc} , \end{eqnarray}$$

(7)

and

$$\begin{eqnarray} R_{acc}=({\rm St}/0.1)^{1/3} R_\mathrm{ H}, \end{eqnarray}$$

(8)

where |$\rm St$| is the Stokes number of the pebbles and R_H is the Hill radius of the planet (Johansen & Lambrechts 2017). Our assumption of Hill accretion can sometimes overestimate the accretion rate when we consider small (⁠|$\lt0.1-1\, \textrm {M}_{{\oplus }}$|⁠) planetary masses and Stokes numbers, where accretion should occur in the Bondi regime (see e.g. fig. 5 in Lyra et al. 2023). Since the difference in accretion rate is small near the transition from Bondi to Hill accretion, and the Bondi regime is rather limited in our simulations, this should not have a large impact on our results. Furthermore, Lyra et al. (2023) demonstrated that the pebble accretion rate changes when considering a distribution of pebble sizes rather than a single size. The polydisperse accretion rate can then be significantly higher than the monodisperse rate in the Bondi regime, and slightly lower in the Hill regime.

The pebble density in the disc’s mid-plane is given by,

$$\begin{eqnarray} \rho _{\rm pe} = \frac{\Sigma _{\rm pe}}{\sqrt{2\pi }H_{\rm pe}}. \end{eqnarray}$$

(9)

We calculate the pebble scale height as |$H_{\rm pe} = H \sqrt{\alpha _{\rm T}/(\alpha _{\rm T} + {\rm St})}$| (Klahr & Henning 1997; Lyra & Lin 2013), where α_T is the turbulent parameter (note, that the turbulent parameter is not the same as the viscous parameter α that regulates the viscous evolution of the gas disc). The pebble surface density is obtained from the continuity requirement |$\Sigma _{\rm pe} = \dot{M}_{\rm pf}/(2\pi r v_{\rm R,pe})$|⁠.

We use two different models for determining the radial pebble flux past the planet location |$\dot{M}_{\rm pf}$| and the Stokes number. In the constant model, we adopt a constant Stokes number and |$\dot{M}_{\rm pf} = Z \times \dot{M}_{\rm disc}$|⁠, where Z is the disc metallicity. These are common assumptions used in the literature; however, as demonstrated by Drążkowska, Stammler & Birnstiel (2021) the growth of planets can change significantly when considering a more realistic dust evolution model. Motivated by this, we also construct an evolving model, where we use pebble fluxes and Stokes numbers from the code pebble predictor (Drążkowska et al. 2021). The pebble predictor is a semi-analytic model which predicts the pebble flux and flux-averaged Stokes number as a function of time and semimajor axis in an arbitrary unperturbed disc. This model is heavily dependent on the assumed fragmentation velocity v_frag and radial extent of the disc R_edge. We consider two different fragmentation velocities, |$v_{\textrm {frag}}=1$| and |$10\, \textrm {m}\mathrm{ s}^{-1}$|⁠, and three different disc extents, |$R_{\textrm {edge}}=50$|⁠, |$100$|⁠, and |$200\, \textrm {au}$|⁠. Further details on the pebble predictor are presented in Appendix  A. Furthermore, when calculating the radial pebble flux towards the inner planet, we remove the pebbles accreted onto the outer planet.

The radial drift velocity of pebbles with Stokes number St «1 can be approximated as,

$$\begin{eqnarray} v_{\rm R,pe} = -2 {\rm St} \eta v_{\rm K} + v_{\rm R,disc}, \end{eqnarray}$$

(10)

where

$$\begin{eqnarray} \eta = -\frac{1}{2}\left(\frac{H}{r}\right)^2 \frac{\partial \ln P}{\partial \ln r}. \end{eqnarray}$$

(11)

In the above expression, v_K is the Keplerian orbital velocity and P = Σ_discT/H is the disc’s pressure.

In our simulations, we allow the accretion of pebbles to continue until the total planetary mass exceeds the pebble isolation mass (M_iso), which can be calculated by (Bitsch et al. 2018),

$$\begin{eqnarray} M_{\rm iso} = 25\, \textrm {M}_{{\oplus }} \times f_{\textrm {fit}}, \end{eqnarray}$$

(12)

where

$$\begin{eqnarray} {f_{\textrm {fit}} = \left[\frac{H/r}{0.05} \right]^3 \left[0.34\left(\frac{\log (\alpha _3)}{\log (\alpha _{\rm T})} \right)^4 + 0.66 \right] \left[1-\frac{\partial \ln P/\partial \ln r + 2.5}{6} \right],} \\ \end{eqnarray}$$

(13)

and α₃ = 10⁻³.

2.3 Gas accretion model

The rate of gas accretion of a growing planet is highly uncertain, and there are various different prescriptions presented in the literature. Rather than picking one of these prescriptions, we chose to derive our own fit for gas accretion based on planet formation simulations performed with the mesa code that was updated to model planetary formation (Valletta & Helled 2020). These simulations are carried out using the same disc parameters and pebble-accretion prescription as our constant model. We simulate 30 random cases and self-consistently calculate the gas accretion rate. The resulting gas accretion can be well represented by an analytical fit,

$$\begin{eqnarray} \dot{M}_{\text{gas,1}} &=& 10^{a_{1}} \left(\frac{M_{\text{core }}}{1 \textrm {M}_{{\oplus }}} \right)^{b_{1}} \left(\frac{\dot{M}_{\text{solid }}}{10^{-7} \textrm {M}_{{\oplus }}/\mathrm{yr}} \right)^{c_1} \textrm {M}_{{\oplus }}/\mathrm{yr}, \end{eqnarray}$$

(14)

and the accretion rate can be obtained as,

$$\begin{eqnarray} \dot{M}_\mathrm{gas} =\dot{M}_\mathrm{gas,1} +\dot{M}_\mathrm{gas,2}. \end{eqnarray}$$

(16)

The parameters of the inferred fit are presented in Table 1. The fit to the gas accretion depends on the planetary core mass, envelope mass, and solid accretion rate. Although we only considered pebble accretion in the mesa simulations, in our model the solid accretion rate is taken as the sum of pebble accretion and planetesimal accretion. The fitting parameters further depend on the chosen grain opacity. We consider two different scaling factors for the contribution of grains in the opacity f_g: 0.1 and 1.0. Further details on the mesa simulations and the construction of the fit are given in Appendix  C.

Table 1.

Parameters of the fitting function for gas accretion.

	fg = 1.0	fg = 0.1
a1	−8.655389	−8.058656
b1	3.488167	3.262527
c1	−0.449784	−0.464667
a2	−10.725292	−11.188670
b2	3.989025	5.834267
c2	2.415257	2.880980
d2	−0.307779	−1.116815

Similar to Bitsch, Lambrechts & Johansen (2015), we limit the gas accretion rate onto the planet to 80 per cent of the disc accretion rate, as it has been shown that not even deep gaps can fully halt the radial flow of gas (Lubow & D’Angelo 2006). Since there is no significant gas accretion expected during the earliest phases of core formation, we turn on gas accretion after the core reaches a mass of |$1\, \textrm {M}_{{\oplus }}$|⁠. Furthermore, since neither Uranus nor Neptune have massive gas envelopes, we only consider gas accretion in the regime M_core > M_env. Whenever a planet reaches M_core < M_env in our simulations, we stop the simulation and record the total masses. Finally, we implement a floor value for the solid (i.e. heavy element) accretion rate of |$\dot{M}_{\rm solid}=10^{-10}\, \textrm {M}_{{\oplus }}\textrm {yr}^{-1}$|⁠, in order to prevent problems with the fit as the gas disc dissipates and |$\dot{M}_{\rm solid} \rightarrow 0$|⁠.

2.4 Planetesimal accretion model

The location and timing of planetesimal formation in the Solar nebula is highly disputed. Population studies often assume a smooth wide-stretched disc of planetesimals, whereas simulations of planetesimal formation tend to promote formation in specific regions of the disc. Examples of such locations are around ice-lines (Drążkowska & Alibert 2017; Schoonenberg & Ormel 2017), at local pressure bumps (e.g. Carrera et al. 2021), and at planetary gap-edges (Stammler et al. 2019; Eriksson, Johansen & Liu 2020). Formation in a wider region of the disc has been shown to be possible in the case of efficient gas removal via photoevaporation (Carrera et al. 2017). In this study, for simplicity, we assume that the initial planetesimal surface density Σ_pl is equal to the planetesimal metallicity (Z_pl) times the gas surface density at the start time of the simulation (t_start). We consider planetesimal metallicities of 0 (core growth only occurs via pebble accretion), 0.25Z, and 0.5Z. The pebble metallicity is then taken to be Z − Z_pl. The time t_start is varied in the parameter study, and since the gas surface density decreases with time, this implies that planets which begin to form late will be surrounded by a planetesimal disc of lower mass than planets which begin to form early.

The planetesimal accretion rate is calculated using the semi-analytic model of Chambers (2006):

$$\begin{eqnarray} \dot{M}_{\rm pl} = \frac{2\pi \Sigma _{\rm pl} R_\mathrm{ H}^2}{P_{\rm orb}} P_{\rm coll}, \end{eqnarray}$$

(17)

where P_orb is the planet’s orbital period and P_coll is the mean collision rate. We use the mean collision rates from Inaba et al. (2001). P_coll depends on the eccentricity and inclination of the planetesimals, and the capture radius of a protoplanet. The evolution of the eccentricity, inclination, and surface density of the planetesimal disc are modelled using the statistical approach (e.g. Inaba et al. 2001). In this study, we use the model developed by Fortier et al. (2013). This calculation takes into account gas drag as well as viscous stirring from the planet and the surrounding planetesimal disc. The collision radius is approximated using equation 7 from Valletta & Helled (2021) when the total mass is above |$2\, \textrm {M}_{{\oplus }}$| and the H–He mass fraction is above 1 per cent. For lower masses and H–He mass fractions, the collision radius is set to be the core radius, calculated using a density of |$1455\, \textrm {kg}\textrm {m}^{-3}$| (an average between the density of Uranus and Neptune). We consider planetesimals of |$100\, \textrm {km}$| in size, and with a bulk density of |$1000\, \textrm {kg}\textrm {m}^{-3}$|⁠. Further details of the planetesimal accretion model are presented in Appendix  D.

2.5 Simulation set-up

We adopt a linear time-grid with a time-step of |$\mathrm{ d}t=500\, \textrm {yr}$|⁠, and use a simple Euler method to update the planetary masses. We initiate our planetary embryos at the same time and with a mass of |$0.01\, \textrm {M}_{{\oplus }}$|⁠. As we discuss below, we also perform simulations where the outer planet is inserted at a higher mass, to mimic an earlier and/or more massive embryo formation. In our constant model, we vary the Stokes number, disc metallicity, fraction of solid mass in planetesimals versus pebbles, turbulent alpha, initial disc accretion rate, time of embryo formation, disc evaporation time, and grain opacity. In our evolving model, we calculate the Stokes number using the peeble predictor code, and instead vary the radial extent of the protoplanetary disc. In the evolving model, we further vary the fragmentation velocity. The values of all parameters that are varied in the parameter study are listed in Table 2. For each configuration we perform ∼500 000 simulations using the constant model, and ∼300 000 simulations when using the evolving model.

Table 2.

Values of used parameters in this study.

St (constant)	5e-3	6e-3	7e-3	8e-3	9e-3
	1e-2	2e-2	3e-2	4e-2	5e-2
Z	5e-3	6e-3	7e-3	8e-3	9e-3
	1e-2	2e-2	3e-2	4e-2	5e-2
Zpl	0	0.25Z	0.5Z	–	–
αT	1.0e-5	2.5e-5	5.0e-5	7.5e-4	1.0e-4
	2.5e-4	5.0e-4	7.5e-4	1.0e-3	–
|$\dot{\mathbf {M}}_\mathrm{\mathbf {0,disc}}$|	1e-8	2e-8	3e-8	4e-8	5e-8
[M⊙yr−1]	6e-8	7e-8	8e-8	9e-8	1e-7
tstart	1e5	1e6	2e6	–	–
[yr]
tdisc	3e6	5e6	10e6	–	–
[yr]
fg	0.1	1	–	–	–
Redge (evolving)	50	100	200	–	–
[au]
vfrag(evolving)	1	10	–	–	–
[ms−1]

In the first part of the paper, we consider in situ formation where Uranus and Neptune are located at |$19.1$| and |$30.0\, \textrm {au}$|⁠, respectively. In the second part of the paper, we vary the planetary formation locations, and also consider the option that the planets could have switched places after the dissipation of the gas disc. For simplicity, we assume that the planets do not migrate or shift location during the disc’s lifetime. Therefore, when varying the formation locations, we do not consider the potential mechanisms that eventually led to the current semimajor axes of Uranus and Neptune at 19.1 and |$30.0\, \textrm {au}$|⁠, respectively. This could have occurred due to a period of dynamical instability after the dispersal of the gas disc, as suggested by the Nice model.

3 IN SITU FORMATION

We first investigate the possibility of forming both planets at their current semimajor axes. The growth-tracks for one example simulation with the parameters listed in Table 3 are shown in Fig. 1, along with the corresponding time evolution of the pebble flux and Stokes number. For these given parameters, it takes |$\sim 10^5\, \textrm {yr}$| for the solid population to grow until maximum size in the evolving model. The growth takes approximately twice as long at Neptune’s location compared to Uranus’ location, resulting in a smaller maximum Stokes number and pebble flux. However, the radial drift is slower at Neptune’s location, which leads to a slower decline of the pebble flux with time. In this example, the pebble flux resulting from the evolving model declines at a slower rate than it does in the constant model. Furthermore, since the time evolution of the gas disc in neglected in the calculation of the Stokes number and pebble flux in the pebble predictor, the pebble flux from the evolving model does not equal zero at t = t_disc.

Figure 1.

The growth of Uranus and Neptune in situ. The top panel shows the growth-tracks for one example simulation, produced using the parameters given in Table 3. The evolving model results in a higher mass for Uranus than the constant model, and a larger mass difference between the two planets. The corresponding pebble flux and Stokes number evolution are shown in the two bottom panels.

Table 3.

Parameters used to produce the example simulation presented in Fig. 1.

‘Representative’ St (constant)	0.0129
Z	0.01
Zpl	0.5Z
αT	10−5
|$\dot{M}_{0,\textrm {disc}}$|	|$9\times 10^{-8}\, \textrm {M}_{\odot }\textrm {yr}^{-1}$|
tstart	|$10^5\, \textrm {yr}$|
tdisc	|$3 \times 10^6\, \textrm {yr}$|
fg	1.0
Redge(evolving)	|$200\, \textrm {au}$|
vfrag(evolving)	|$10\, \textrm {ms}^{-1}$|

In order to demonstrate how the growth-tracks change when switching between the constant model and the more realistic evolving model, we calculate a ‘representative’ Stokes number and use that in the constant model. We obtained this ‘representative’ Stokes number as follows. We calculate the time-average of the flux-weighted Stokes number at both planet locations; and take the average of the two numbers. The resulting growth-tracks are presented in the top panel of Fig. 1. The evolving model produces a rather good Uranus analogue; however, the corresponding mass of Neptune is just above |$1\, \textrm {M}_{{\oplus }}$|⁠. The constant model results in a much lower mass for Uranus, and in general a smaller mass difference between the two planets. The growth-tracks for a second example simulation are shown in Appendix  B, where we further demonstrate how the planet growth changes when considering: only pebble accretion; pebble and gas accretion; and finally pebble, gas, and planetesimal accretion.

In Fig. 2, we show the outcome of all the simulations when assuming in situ formation. The scatter points indicate the total masses of all the planets with a H–He mass fraction <20 per cent at the time of disc dissipation (we calculate the H–He mass fraction as the ratio of gas mass to total mass). Although we used a wide range of parameters, we do not produce any Uranus and Neptune analogues. Most models do manage to produce a Uranus analogue, although they are rare, but the corresponding mass of Neptune in these cases is always well below |$10\, \textrm {M}_{{\oplus }}$|⁠. This is not unexpected, since both pebble and planetesimal accretion tend to be more efficient at smaller semimajor axes.

Figure 2.

Plots showing the results of all our in situ simulations where the H–He mass fraction is <20 per cent at the time of disc dissipation. Simulations in which the H–He mass fraction becomes larger than 20 per cent are not shown on the plots. The masses of Uranus and Neptune are marked with an asterix. Although we consider a wide range of parameters, no Uranus and Neptune analogues are found.

When a fragmentation velocity of |$1\, \textrm {m}\textrm {s}^{-1}$| is being used (middle column), the dust-evolution model struggles to grow Neptune above |$1\, \textrm {M}_{{\oplus }}$|⁠. This problem becomes less severe when the fragmentation velocity is increased, and in this scenario there are a few cases when Neptune grows to become more massive than Uranus. This is either an effect of the dust-evolution model, and/or due to the blocking of pebbles by Neptune. When we use f_g = 0.1, the gas fractions become too high in the mass range relevant to Uranus and Neptune. As the fraction of mass in planetesimals versus pebbles increases, the total masses of planets with H–He mass fractions below 20 per cent in the constant simulation decreases. This is expected since pebble accretion tends to be more efficient than planetesimal accretion at large semimajor axes. This trend is harder to spot when looking at the results from the evolving model, since the simulation outcome in general is much more variable.

In Fig. 3, we show more detailed results from one of the simulation sets (corresponding to the scatter points labeled Zpl = 0.50Z in the bottom right panel of Fig. 2, in total ∼24 000 simulations). The colour of the scatter points represents H–He mass fractions of the planets at t = t_disc, and the grey crosses indicate the total masses of the planets at the time one of them reaches a H–He mass fraction above 50 per cent. None of these simulations resulted in a planet of Uranus mass or larger with H–He mass fractions below 20 per cent. Furthermore, no simulation leads to an outer planet with mass above |$5\, \textrm {M}_{{\oplus }}$|⁠, while keeping the H–He mass fractions of both planets below 50 per cent. In Fig. 3, we find planets with total masses as low as |$7\, \textrm {M}_{{\oplus }}$| that reach H–He mass fractions of 50 per cent. If allowed, most of these planets would grow to become gas-giants.

Figure 3.

Total masses and H–He mass fractions at the time of disc dissipation for all in situ simulations with parameters as indicated by the plot title. The crosses show simulations where one of the planets obtained a H–He mass fraction above 50 per cent before disc dissipation. If growth were allowed to continue, most of these planets would grow to become gas giants. The growth-tracks for the simulation marked Fig. 1 is shown in Fig. 1.

In summary, the results of this section suggest that it is unlikely that Uranus and Neptune formed in situ at their current locations, if the embryos formed at the same time and with the same mass. In Section 5.1, we investigate whether this result changes when we allow for different embryo masses.

4 VARYING THE FORMATION LOCATION

In the previous section we investigated the possibility of forming both planets at their current locations. However, the initial formation locations of the planets is unknown, and it is likely that some migration took place during and/or after the dissipation of the gas disc. Assuming that both planets formed beyond the current orbit of Saturn, the formation locations could have been somewhere between |$\sim 12-40\, \textrm {au}$|⁠. The inner boundary comes from requiring dynamical stability with Saturn, and the outer boundary comes from demanding that no planet disturbs the classical Kuiper Belt. We consider the following formation locations for the inner planet: 12, 15, 18, 21, 24, 27, and |$30\, \textrm {au}$|⁠. We then place the outer planet at 5, 10, 15, and |$20\, R_\mathrm{ H}$| beyond this location, where we used the current mass of Neptune to calculate the Hill radii. This leads to formation locations beyond |$40\, \textrm {au}$| in three cases, which we subsequently remove from the study. Due to the complicated nature of planetary migration, it is not included in our simulations. Instead, we assume that migration occurred after the dissipation of the gas disc. The possible impacts of planetary migration on our results are discussed in Section 5.2.

We consider successful Uranus and Neptune analogues to have masses within |$1.5\, \textrm {M}_{{\oplus }}$| of their current masses and H–He mass fractions <20 per cent. In Appendix  E, we show how the number of successful analogues changes when allowing the planetary masses to differ by |$3\, \textrm {M}_{{\oplus }}$| instead of |$1.5\, \textrm {M}_{{\oplus }}$|⁠. Since it is possible that Uranus and Neptune have switched places after the dissipation of the disc, we also consider this scenario when identifying successful analogues. In other words, we also consider simulations where the mass of the inner planet is within |$1.5\, \textrm {M}_{{\oplus }}$| of Neptune’s current mass, and the mass of the outer planet is within |$1.5\, \textrm {M}_{{\oplus }}$| of Uranus’ current mass, to be successful.

Fig. 4 shows a few examples from the parameter study, where we have used circles to highlight successful analogues. The top left panel shows our most compact configuration, with the inner planet at |$12\, \textrm {au}$| and the outer planet at |$13.5\, \textrm {au}$|⁠. In this scenario, the constant model can produce successful Uranus and Neptune analogues. There are also cases where the blocking of pebbles by the outer planet, leads to the outer planet being more massive than the inner planet. When switching to the more realistic evolving model, we find no successful analogues, since the mass of the outer planet becomes too low in the relevant mass range. In order for these planets to become successful analogues, the outer planet would need to acquire some additional mass via other mechanisms, such as giant impacts (see discussion in Section 5.3). When increasing the orbital distance between the planets, the mass difference increases, resulting in fewer successful analogues. In this scenario, multiple giant impacts would be required to produce Uranus and Neptune analogues.

Figure 4.

Plots showing the total masses of planets with H–He mass fractions <20 per cent, for simulations with four different semimajor axes configurations (the semimajor axes of the planets are indicated in the x and y-labels). The legends on the top left panels are the same in the other three panels. The top left panel shows results obtained with our most compact configuration, and the bottom right panel shows results obtained with a very non-compact configuration.

The number of successful Uranus and Neptune analogues obtained with each planet configuration is presented in Fig. 5. The number of simulations behind each grid-cell in the histograms is ∼80 000 for the constant model, and ∼24 000 for the evolving model. When using the constant model that has a constant Stokes number and a pebble-flux proportional to the disc accretion rate, we find that successful analogues can be produced when the orbital distance between the planets is equal to 5 Neptune Hill radii. The successful simulations have in common Z ≥ 0.03, |$\dot{M}_{0,\textrm {disc}} \ge 5\times 10^{-8}\, \textrm {M}_{\odot }\textrm {yr}^{-1}$|⁠, |$t_{\rm start}=1\, \textrm {Myr}$|⁠, and |$t_{\rm disc}=3\, \textrm {Myr}$|⁠. The number of successful analogues decreases with increasing semimajor axes and planetesimal-to-pebble mass fractions.

Figure 5.

Histograms showing the number of Uranus and Neptune analogues found for each planetary configuration. The top left panel is from simulations with no planetesimal accretion, the top right panel is from simulations with a planetesimal-to-solid mass fraction of 25 per cent, and the bottom panel is from simulations with a planetesimal-to-solid mass fraction of 50 per cent. The results show that the constant model can produce Uranus and Neptune analogues when the separation between the two planets is small. When considering the more realistic evolving model, no analogues are found.

Our formation model does not lead to Uranus and Neptune analogues when using the dust-evolution model, regardless of the formation locations that are being used. The main reasons for this are as follows: (1) the H–He mass fractions of planets in the Uranus and Neptune mass range are typically found to be higher than 20 per cent; and (2) the growth efficiency strongly depends on the semimajor axis, such that a small difference in semimajor axis still generates a relatively large difference in planetary mass. Therefore, our results suggest that it is unlikely for Uranus and Neptune to have formed solely via pebble, gas, and planetesimal accretion, if the embryos formed at the same time and with the same mass. However, we would like to stress that although we have introduced a rather advanced and comprehensive formation model, simplifications have still been made. For example, we do not consider how the heavy elements mix within the planetary atmosphere, which could affect the accretion rate of gas onto the planet (see Section 5.5). We also assume that the initial planetesimal disc is wide-stretched and that the planetesimals are single sized, whereas in reality the surface density of planetesimals could be in-homogeneous within the disc, and the planetesimals are expected to have a size distribution, which is also time dependent. Finally, there are additional processes we have not considered such as giant impacts and planetary migration which could affect our conclusions. We discuss these processes in the following section.

5 POTENTIAL SOLUTIONS

We consider the growth of Uranus and Neptune via pebble, gas, and planetesimal accretion. We can form both planets in our simulations when we assume a constant Stokes number, a pebble flux that is proportional to the disc accretion rate, and a small orbital separation between the planets. The first two assumptions are often used in pebble accretion simulations; but as we have demonstrated, the growth via pebble accretion can significantly change when considering a more realistic model where the dust population evolves with time and semimajor axis. When we switch to using the evolving model, we no longer find any successful Uranus and Neptune analogues, despite the large range of formation locations and parameters that were used. In this section, we discuss potential mechanisms that could assist in forming Uranus and Neptune simultaneously.

5.1 What if the outer planet had a head start?

In our main study, we assumed that the embryos of Uranus and Neptune formed at the same time and with the same mass. However, since the efficiency of pebble accretion strongly depends on the semimajor axis and the planetary mass, this leads to the outer planet being much less massive than the inner one. This mass difference could significantly decrease if the outer embryo formed earlier than the inner one, or similarly, if the outer embryo is more massive than the inner one. Since the process of planet formation is rather stochastic, it is certainly possible that embryos form at different times and with different masses.

We perform additional simulations where the mass of the outer embryo is increased by a factor of 10 and 100 compared to the inner embryo (with a mass of |$0.01\, \textrm {M}_{{\oplus }}$|⁠). This difference in mass could be due to an earlier embryo formation and/or a more massive embryo being formed. We limit this study to the case when all solid accretion occurs via pebble accretion (Z_pl = 0). Fig. 6 shows how the results for the case of in situ formation change when the mass of the outer embryo is increased. The number of successful analogues that are obtained at various formation locations are shown in Fig. 7.

Figure 6.

Total mass of planets that are forming in situ and has a H–He mass fraction <20 per cent at the time of disc dissipation. Blue scatter points show the results when the planetary embryos have the same mass; red scatter points show the results when the outer embryo is 10 times more massive than the inner one; and yellow scatter points show the results when the outer embryo is 100 times more massive than the inner one. These simulations are performed without planetesimal accretion.

Figure 7.

Histograms showing the number of Uranus and Neptune analogues that are found when the outer embryo is 10 times more massive than the inner embryo (left panel), and when the outer embryo is 100 times more massive than the inner one (right panel). The simulations are performed without planetesimal accretion. In contrast to when the embryos had the same mass (Fig. 5), we now produce a few analogues while using the evolving model.

When considering in situ formation, the constant model with the default grain opacity produces Uranus and Neptune analogues when the outer embryo is 10 times more massive than the inner one. When the difference in embryo mass is further increased, the outer planet becomes too massive compared to the inner one. When we use the more realistic evolving model, we do not obtain any successful analogues; however, we are much closer to doing so than we were when we used embryos of the same mass.

When we vary the formation locations, we find a few successful Uranus and Neptune analogues, both while using the constant model and the more realistic evolving model. The successful analogues have in common |$\textrm {R}_{\rm edge} \ge 100\, \textrm {au}$|⁠, Z ≥ 0.03, α_T < 10⁻⁴, |$\dot{M}_{0,\textrm {disc}} \ge 5\times 10^{-8}\, \textrm {M}_{\odot }\textrm {yr}^{-1}$|⁠, and |$t_{\rm disc}=3 \times 10^6\, \textrm {yr}$|⁠. Unlike when the embryos formed with the same mass, the required separation between the planets is now typically around |$15-20\, \textrm {R}_{\rm H}$|⁠, rather than |$5\, \textrm {R}_{\rm H}$|⁠. If we are more generous when searching for successful analogues, and allow the mass difference to be |$3\, \textrm {M}_{{\oplus }}$| instead of |$1.5\, \textrm {M}_{{\oplus }}$|⁠, the number of successful analogues is more than doubled (see Appendix  E). The parameters leading to successful analogues varies more in this scenario, but they still have in common |$t_{\rm disc}=3 \times 10^6\, \textrm {yr}$|⁠. This is the case for all successful analogues found throughout the study, and the reason is that when the disc lifetime is longer, continuous gas accretion results in H–He mass fractions above 20 per cent.

In summary, we can simulate the formation of Uranus and Neptune solely via pebble, gas, and planetesimal accretion if the outer planetary embryo forms earlier and/or more massive than the inner one. The required disc metallicity is above 1 per cent, and overall, very specific parameters are required in order to form Uranus and Neptune. Therefore, the formation of Uranus and Neptune remains a challenge to planet formation theories.

5.2 What if the planets migrated?

The effects of disc-driven migration on the formation of Uranus and Neptune are hard to predict, since the migration of the embedded planets can be directed both inwards and outwards depending on the surrounding disc conditions (Paardekooper et al. 2022). Let us assume that the planets migrate inwards during formation via classical Type-I migration (see e.g. equations 3–4 of Johansen, Ida & Brasser 2019). If the negative radial gradients of gas surface density and temperature are 15/14 and 3/7, respectively, and we assume an unperturbed gas surface density, the migration rate has a weak negative dependence on semimajor axis, and is directly proportional to the planetary mass. Since in our models the inner planet typically grows faster than the outer planet, the inner planet would thus migrate faster inwards than the outer planet. This would lead to a faster growth for the inner planet, since pebble accretion is more efficient at smaller semimajor axes, leading to an even larger inward migration, etc. In that case, the effect of migration would enhance the mass difference between the planets in our simulations, making it even more challenging to form Uranus and Neptune.

In the above line of argument we assumed an unperturbed disc. If the planets open up a gap in the gas surface density profile, it could affect the migration rate (see e.g. Kanagawa, Tanaka & Szuszkiewicz 2018). However, since a planet beyond |$10\, \textrm {au}$| with a mass |$\lt20\, \textrm {M}_{{\oplus }}$| is not expected to open up a deep gap, our assumption should be justified. Torques exerted by the dust disc and significantly different temperature and surface density structures could change the migration rate compared to the classical Type-I rate. The migration of planets can further be halted by trapping in resonances with other planets.

Planetary migration is further expected to affect the accretion efficiency of planetesimals. However, it is difficult to assess the magnitude of this effect. A migrating planet can increase the surface density of planetesimals Σ_pl in its vicinity, because planetesimals that are initially outside the feeding zone can enter the feeding zone (Tanaka & Ida 1999; Alibert et al. 2005; Shibata, Helled & Ikoma 2020, 2022; Turrini et al. 2021). However, mean motion resonances of a massive protoplanet like Jupiter were found to prevent planetesimals from entering the planet’s feeding zone (Shibata et al. 2020, 2022), which reduces Σ_pl. Mean motion resonances of a Neptune-sized planet are not as strong as those of a massive giant planet, but it could still change the planetesimal accretion rate. In addition, gravitational scattering of other protoplanets can affect the configuration of mean motion resonances (e.g. Tanaka & Ida 1997; Levison, Thommes & Duncan 2010). The mutual gravitational interaction of the forming Uranus and Neptune could affect the planetesimal accretion rate if they experienced radial migration.

To summarize, it is hard to predict how disc-driven migration would have affected the formation of Uranus and Neptune, but it is reasonable to assume for some sort of disc-driven migration to have taken place.

5.3 Giant impacts

Our results show that it is very hard to form Uranus and Neptune simultaneously, and in the few cases where we do succeed, the required disc metallicity is above the typically assumed 1 per cent. The main challenge concerns keeping the H–He mass fractions and mass difference between the two planets small. One possible solution is giant impacts. The H–He mass fractions of lower mass planets are typically below 20 per cent, and Uranus and/or Neptune could have formed by colliding two or more such planets together. For example, in the cases where we form an Uranus analogue, but the corresponding mass of Neptune is too low, a giant impact onto Neptune could deliver the missing heavy-element mass while not increasing the H–He mass fraction. The likelihood for such an event to have taken place that far out in the Solar system is unknown, although collisions are expected to be common (e.g. Izidoro et al. 2015; Chau et al. 2021).

5.4 The effect of disc substructures

In our model, we assume that the gaseous disc is smooth and evolves as a viscous accretion disc. However, observations of protoplanetary discs have revealed that many discs harbour substructures (e.g. Andrews et al. 2018). The most commonly observed substructure in the dust distribution is rings and gaps, where several observed rings have been shown to be consistent with dust trapping inside pressure maxima (e.g. Dullemond et al. 2018). Such pressure maxima could be the result of some fluid dynamics process, icelines or planet-disc interactions (see Bae et al. 2023 for a recent review on disc substructures).

Since the surface density and radial drift of pebbles in a structured disc differ considerably compared to a smooth disc (Eriksson et al. 2020), the growth of planets via pebble accretion would also be affected. Similarly, local variations in the temperature and gas surface density structure could affect the direction and speed of migration, and result in migration traps (Guilera & Sándor 2017). Because pressure maxima collects large amounts of pebbles, they can be efficient sites for planetesimal formation and further planetary growth (Guilera et al. 2020; Chambers 2021; Lau et al. 2022; Jiang & Ormel 2023).

Indeed substructures in the disc would affect the formation history of the growing planets. If Uranus and/or Neptune were to form at the location of a pressure maxima, their growth could occur at much shorter time-scales than inferred in this work. At the same time, the pebble flux interior to a pressure maxima could also be significantly smaller, potentially leading to slower growth for planets that reside closer to the star than the pressure maxima. The effect of disc substructure on the growth of Uranus and Neptune is complex and should be investigated in detail in future research.

5.5 Gas accretion rate uncertainties

Our calculations of the gas accretion rates are self-consistent and represent a significant improvement in comparison to the commonly used simplifying assumption that gas accretes as a constant fraction of the solid accretion rate. Nevertheless, there are several simplifications in our gas accretion model that could be improved in future work. For example, our models do not include the interaction of the solid material in the envelope. Pebbles and planetesimals, whether they are made of rock or ice are expected to vapourize and enrich the envelope with heavier elements (e.g. Pollack et al. 1986; Podolak, Pollack & Reynolds 1988). This pollution of the envelope with heavy elements can accelerate the planetary growth, making giant planet formation more efficient (e.g. Stevenson 1982; Hori & Ikoma 2011; Venturini et al. 2015; Valletta & Helled 2020). If this is the case, forming Uranus and Neptune would be more difficult in our current model.

However, this mechanism is not completely understood and has received attention primarily in the context of Jupiter’s formation. It remains possible that a more realistic gas accretion prescription actually leads to lower gas accretion rates. This is suggested from three-dimensional gas accretion models that are found to have lower gas accretion rates than one-dimensional models (Ormel, Shi & Kuiper 2015; Cimerman, Kuiper & Ormel 2017), although this is predicted for planets forming at much shorter orbital periods than Uranus and Neptune (e.g. Moldenhauer et al. 2021, 2022 calculated recycling rates at 0.1 au). Pebble enrichment of the nebular gas in the locations where the planets form could also contribute to recycling (Wang et al. 2023). It is therefore clear that inferring realistic gas accretion rates for planets forming at large orbital distances, such as Uranus and Neptune, are desirable. Currently, it is unknown how such enrichment could affect the formation of Uranus and Neptune.

5.6 Planetesimal accretion rate uncertainties

To calculate the planetesimal accretion rate, we adopt the statistical model. However, this model has been developed for modelling the formation of terrestrial planets and cores of giant planets. The effect of the gas accretion is not included in the current statistical model. For example, Shibata, Helled & Kobayashi (2023) showed that the current statistical approach cannot reproduce the results of N-body simulations once the protoplanet enter the runaway gas accretion. Uranus and Neptune have not entered the runaway gas accretion, but the steady gas accretion during the planetary growth could affect the planetesimal accretion rate, and we hope to address this in future research.

6 CONCLUSIONS

We have studied the formation of Uranus and Neptune via pebble, gas, and planetesimal accretion. We considered two different models for pebble accretion: a simple model with constant Stokes number and a pebble flux proportional to the disc accretion rate; and a more realistic model where the Stokes number and pebble flux are obtained from a dust-evolution model and vary with time and semimajor axis. We do not include migration, but test a wide range of formation location and a wide range of disc parameters. Our main conclusions can be summarized as follows:

• If the embryos form at the same time and with the same mass, our formation model with an evolving dust population is unable to produce Uranus and Neptune analogues, regardless of the assumed formation locations. When we use the simpler and less realistic constant model, Uranus and Neptune analogues can form when the orbital distance between the planets is small.

• If the outer embryo forms earlier and/or more massive than the inner embryo, we can form both planets simultaneously in a few instances where the disc is metal-rich and has a lifetime of a few Myr. Overall, very specific parameters are required to form Uranus and Neptune, and the formation of these planets remains a challenge to planet formation theory.

• Based on our results, it is unlikely for Uranus and Neptune to have formed in situ. When we use the evolving model, we do not produce any analogues regardless of the assumed formation locations and embryo masses. In the constant model, we find Uranus and Neptune analogues when the outer embryo is 10 times more massive than the inner embryo.

• The key challenge in forming Uranus and Neptune is keeping the H–He mass fractions below 20 per cent and keeping the planetary masses similar. When the grain opacity is low or the disc lifetime is longer than |$\sim 3\, \textrm {Myr}$|⁠, the H–He mass fractions become too large in the mass regime relevant for Uranus and Neptune. The mass difference between the planets increases with the orbital separation between the planets, since the pebble accretion rate strongly depends on the semimajor axis.

Our study demonstrates the complexity of modelling planet formation properly. In addition, it clearly shows how different model assumptions affect the results concerning the forming planets. This is not only important for improving our understanding of the origin of Uranus and Neptune, but also of the origin of the many intermediate-mass exoplanets detected in our galaxy.

Finally, we suggest that more accurate determinations of the H–He mass fractions in Uranus and Neptune would be valuable in constraining their formation path. We therefore look forward to a future mission to Uranus and Neptune as well as to the upcoming accurate measurements of mass, radius, and atmospheric compositions of intermediate-mass/size exoplanets.

ACKNOWLEDGEMENTS

The authors wish to thank the anonymous referee for providing helpful comments that lead to an improved manuscript. LE acknowledges funding by the Institute for Advanced Computational Science Post-doctoral Fellowship. LE would further like to thank Stony Brook Research Computing and Cyber Infrastructure, and the Institute for Advanced Computational Science at Stony Brook University for access to the SeaWulf computing system, made possible by grants from the National Science Foundation (number 1531492 and Major Research Instrumentation award number 2215987), with matching funds from Empire State Development’s Division of Science, Technology and Innovation (NYSTAR) programme (contract C210148).

MML, SS, and RH acknowledge that part of this work has been carried out within the framework of the National Centre of Competence in Research Planets supported by the Swiss National Science Foundation under grant numbers 51NF40_182901 and 51NF40_205606. The authors acknowledge the financial support of the SNSF. RH and SS also acknowledge support from SNSF under grant number 200020_188460.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

1

All elements heavier than helium.

APPENDIX A: DUST EVOLUTION WITH THE pebble predictor

The pebble predictor code was presented in Drążkowska et al. (2021), and is publicly available at https://github.com/astrojoanna/pebble-predictor. It is a semi-analytic model for predicting the flux-averaged Stokes number and total flux of pebbles at all locations and at all times in an arbitrary unperturbed disc. The pebble predictor takes as input the initial radial surface density profiles of gas (Σ_{disc, 0}) and dust (Σ_{dust, 0} = ZΣ_{disc, 0}), the radial temperature structure (T), turbulent strength (α_T), internal density of dust grains (ρ_•), and fragmentation velocity (v_frag). The initial size of the dust grains is assumed to be |$1\, \mu \textrm {m}$|⁠, and growth proceeds via turbulence driven collisions until either the fragmentation or radial drift barrier is reached. The time evolution of the gas disc is neglected in the pebble predictor model, and the evolution of the dust disc is approximated by keeping track of the dust movement. Despite these simplifications, the results of the pebble predictor can reproduce rather well the outcome of full coagulation simulations (see fig. 7 in Drążkowska et al. 2021).

We consider two different fragmentation velocities, |$v_{\textrm {frag}}=1$| and |$10\, \textrm {m}s^{-1}$|⁠, and assume |$\rho _{\bullet }=1000\ \, \textrm {kg} \ \textrm {m}^{-3}$|⁠. We use a logarithmic radial grid with inner edge at |$1\, \textrm {au}$|⁠, outer edge at 50, 100, or |$200\, \textrm {au}$|⁠, and 100, 200, or 400 grid points, respectively. The time grid has 10 000 logarithmically spaced grid points and stretches from |$1\, \textrm {yr}$| to |$10\, \textrm {Myr}$|⁠. In Fig. A1, we show the time evolution of the pebble flux and Stokes number at the current semimajor axes of Uranus and Neptune, for the set of disc parameters given in Table B1. We further show how the results vary with the assumed fragmentation velocity, and how they change when considering that half of the solid mass is locked up in planetesimals (bottom panels).

Figure A1.

The time evolution of the pebble flux (left) and Stokes number (right) that is obtained from the pebble predictor when using the parameters given in Table B1. The pebble flux obtained from taking the product of the disc metallicity and the disc accretion rate (the constant model) is shown with black lines in the left panels. The black lines in the right panels show a ‘representative’ Stokes number (see text for details).

When we use the higher of the two fragmentation velocities, the pebbles grow to significantly larger sizes. The corresponding pebble flux also peaks at higher values; however, because of this the pebble reservoir is also depleted more quickly. When comparing with the pebble flux that results in from taking the product of the disc metallicity and the disc accretion rate (the constant model, black lines), the pebble predictor in this case results in lower pebble fluxes at late times. In the constant model, we use a Stokes number that is constant with time and semimajor axes. The black line in the right panels show a ‘representative Stokes number’, which is obtained after: calculating the time-average of the flux-weighted Stokes number; and taking the average of the resulting values obtained at the two semimajor axes and with the two fragmentation velocities. We use this representative Stokes number to compare the growth tracks obtained with the a model and the evolving model in Appendix  B.

APPENDIX B: GROWTH-TRACKS FROM AN EXAMPLE SIMULATION

Fig. B1 presents the growth-tracks that are obtained when using the simulation parameters specified in Table B1. We show how these growth-tracks vary when considering growth via only pebble accretion (top row); growth via pebble and gas accretion (middle row); and growth via pebble, gas, and planetesimal accretion when assuming that 50 per cent of the solid mass is locked up in planetesimals (bottom row). Furthermore, we demonstrate how the results change when switching between the constant model (right panels) and the more realistic evolving model (left panels). A comparison between the pebble fluxes and Stokes numbers used in each of the two models is shown in Fig. A1.

Figure B1.

Examples of growth-tracks when the planets form in situ, produced using the parameters given in Table B1. The different panels demonstrate how the growth-tracks vary when different accretion mechanisms are considered, and when we switch between the constant and evolving models.

Table B1.

Parameters of the example case presented in Figs A1 and B1.

‘Representative’ St (Zpl = 0)	0.021
‘Representative’ St (Zpl = 0.5Z)	0.015
Z	0.02
αT	5 × 10−5
|$\dot{M}_{0,\textrm {disc}}$|	|$6\times 10^{-8}\, \textrm {M}_{\odot }\textrm {yr}^{-1}$|
tstart	|$10^5\, \textrm {yr}$|
tdisc	|$3 \times 10^6\, \textrm {yr}$|
Redge (evolving)	|$100\, \textrm {au}$|

Let us begin with analysing the case when all growth occurs via pebble accretion, and the pebble flux and Stokes number is obtained from the pebble predictor code (top left panel). Both planets grow several times more massive when we increase the fragmentation velocity from 1 to |$10\, \textrm {ms}^{-1}$|⁠. The mass of the outer planet is several times lower than the mass of the inner planet. When we switch to using the constant model with a constant Stokes number and a pebble flux proportional to the disc accretion rate (top right panel), the planets grow more massive. This is mostly because the pebble flux does not decrease as fast with time as it does in the more realistic evolving model. We also show how the growth of the inner planet becomes slower, when we remove the pebbles that are accreted onto the outer planet from the flux towards the inner planet.

When we take into account the accretion of gas onto the planetary cores (middle panels), the growth-tracks change drastically. The inner planet obtains a H–He mass fraction higher than 50 per cent and is removed from the simulation, regardless of the pebble model and grain opacity that is being used. The outer planet suffers the same fate when we use the lower grain opacity, despite the core mass being lower than |$5\, \textrm {M}_{{\oplus }}$| in the pebble predictor model. This is a common outcome of models where we use the lower grain opacity. This demonstrates the need for a proper treatment of gas accretion in planet formation simulations during the core accretion phase.

When we trap half of the solid mass in planetesimals (bottom panels), the planetary growth is less efficient. The total amount of accreted planetesimal mass onto the planets is much less than |$1\, \textrm {M}_{{\oplus }}$|⁠. The cause for the difference in planetary growth between the case with and without planetesimal accretion is due to the difference in pebble mass. When we trap a significant fraction of the solid mass in planetesimals, the pebble flux decreases and the planetary growth is limited.

APPENDIX C: DETAILS OF THE GAS ACCRETION MODEL

To calculate the self-consistent gas accretion rates we use a separate evolution model for the gas accretion (see Valletta & Helled 2020, for details). We use mesa version 10108 and the mesa sdk of the same version.

Gas accretion rates are calculated in the following way. Every time-step we compute the accretion radius, which is set by,

$$\begin{eqnarray} R_{\rm acc} = \frac{GM_{\rm p}}{c_{\mathrm{ s}}^{2}/k_{1} + G M_{\rm p}/\left(k_{2} R_{\mathrm{ H}} \right)}, \end{eqnarray}$$

(C1)

where G is the universal gravitational constant, M_p is the mass of the protoplanet, c_s is the sound speed in the disc at the location where the planet is forming, and R_H is the protoplanet’s Hill radius. k₁ and k₂ are constants to account for limited supply of gas at the protoplanet’s formation location due to disc perturbations. These are set to |$\frac{1}{2}$| and |$\frac{1}{4}$| respectively, which were found suitable values by Lissauer et al. (2009), albeit based on simulations of Jupiter’s formation. The values of k₁ and k₂ could be close to unity for less massive planets that do not perturb the disc; however, there can be other effects which reduce gas accretion (see Section 5.5). We then add gas to the envelope until the radius of the protoplanet equals the accretion radius. To construct the fit, we simulate the growth of 30 planets with varying disc and planet parameters, that all have final total masses between 1 and |$30\, \textrm {M}_{{\oplus }}$| if growth occurs via only pebble accretion with the constant model.

The gas accretion rate depends on three things. First, as the planetary mass increases the accretion radius is extended. Second, the luminosity of the growing core, which is set by the pebble accretion rate, inflates the envelope. Therefore, less gas is needed to fulfil the accretion criterion of R_acc = R_p. Third, the assumed grain opacities can notably limit the cooling efficiency of the accreted gas, which can further enhance the effect that the luminosity has on the gas accretion. Following Valencia et al. (2013) we combine Rossland opacities from Freedman et al. (2014; κ_mol), and grain opacities based on tables from Alexander & Ferguson (1994; κ_grain). So that the opacity is,

$$\begin{eqnarray} \kappa _{\textrm {tot}} = \kappa _{\textrm {mol}} + f_{\textrm {g}} \times \kappa _{\textrm {grain}} . \end{eqnarray}$$

(C2)

We consider cases where f_g = 1 and f_g = 0.1.

Since the original mesa data numerically oscillate in every step, we average the original data for every 100 steps and use the averaged data for fitting functions. Upper panels in Fig. C1 show the averaged gas accretion rate. When the envelope’s mass is small, the gas accretion rate depends on the core mass M_core. As the envelope’s mass increases, self gravity of the planetary envelope also regulates the gas accretion rate. From the mesa simulations, we find that the gas accretion rate starts to depend on the envelope mass M_env once the envelope-core mass fraction f_e/c = M_env/M_core exceeds ∼0.1. In order to follow each characteristic of gas accretion, we divide the data set into two sub-data sets using f_e/c; the data set with f_e/c < 0.1, and the data set with f_e/c > 0.5. The former data set are used for fitting equation (14), and the latter one are used for fitting equation (15), respectively. We fit the functions to each data set in logarithm space because the gas accretion rate changes orders of magnitude during the formation of Uranus and Neptune.

Figure C1.

Upper panels: Gas accretion rates obtained by mesa. Each line corresponds to one simulation, pebble accretion was modelled using the constant model. Left and right panels show the cases with f_g = 1.0 and f_g = 0.1, respectively. Lower panels show comparison of the obtained fitting function with the original data for 4 cases out of the 30 simulations, for each f_g. The solid lines show the fitting function equation (16). Circle plots are the averaged data of mesa’s simulations.

The lower panels in Fig. C1 show the comparison of mesa and the obtained fitting function equation (16). We do not fit the function when the gas accretion regimes shifts from equation (14) to equation (15; 0.1 < f_c/e < 0.5), but the summation of equations (14–15) [equation (16)] can reproduce mesa’s results well. We use equation (16) in all the simulations used in this study.

APPENDIX D: PLANETESIMAL ACCRETION MODEL

The mean collision rate in equation (17) can be written as (Inaba et al. 2001),

$$\begin{eqnarray} P_{\rm coll} = {\rm min} \left(P_{\rm med}, \left({P_{\rm high}}^{-2} +{P_{\rm low}}^{-2} \right)^{-1/2} \right), \end{eqnarray}$$

(D1)

where

$$\begin{eqnarray} P_{\rm high} &=& \frac{{\tilde{r}}^2}{2 \pi } \left(I_{\rm F} (\beta) +\frac{6 I_{\rm G} (\beta)}{\tilde{r} {\tilde{\mathrm{ e}}}^2 } \right), \end{eqnarray}$$

(D2)

$$\begin{eqnarray} P_{\rm med} &=& \frac{{\tilde{r}}^2}{4 \pi \tilde{i}} \left(17.3 +\frac{232}{\tilde{r}} \right), \end{eqnarray}$$

(D3)

$$\begin{eqnarray} P_{\rm low} &=& 11.3 {\tilde{r}}^{1/2}, \end{eqnarray}$$

(D4)

with

$$\begin{eqnarray} \tilde{e} &=& a_\mathrm{ p} \mathrm{ e}/R_{\rm H}, \end{eqnarray}$$

(D5)

$$\begin{eqnarray} \tilde{i} &=& a_\mathrm{ p} \mathrm{ i}/R_{\rm H}, \end{eqnarray}$$

(D6)

$$\begin{eqnarray} \tilde{r} &=& R_{\rm cap}/R_{\rm H}, \end{eqnarray}$$

(D7)

where e and i are mean eccentricity and inclination of planetesimals, a_p is the semimajor axis of the protoplanet, R_cap is the capture radius of the protoplanet. I_F and I_G are given by (Chambers 2006),

$$\begin{eqnarray} I_{\rm F} (\beta) &\simeq& \frac{1+0.95925 \beta +0.77251 \beta ^2}{\beta (0.13142 +0.12295 \beta)}, \end{eqnarray}$$

(D8)

$$\begin{eqnarray} I_{\rm G} (\beta) &\simeq& \frac{1 +0.3996 \beta }{\beta (0.0369 +0.048333 \beta +0.006874 \beta ^2)}, \end{eqnarray}$$

(D9)

where |$\beta =\tilde{i}/\tilde{e}$|⁠. The evolution of |$\tilde{e}$| and |$\tilde{i}$| is determined by the viscous stirring between planetesimals, the viscous stirring from the protoplanet, and the gas drag from the gaseous disc. The rates of the changes in the eccentricity and inclination of planetesimals are given by,

$$\begin{eqnarray} \frac{{\rm d} \mathrm{ e}^2}{{\rm d} t} &=& \left. \frac{{\rm d} \mathrm{ e}^2}{{\rm d} t} \right|_{\rm drag} +\left. \frac{{\rm d} \mathrm{ e}^2}{{\rm d} t} \right|_{\rm VS,M} +\left. \frac{{\rm d} \mathrm{ e}^2}{{\rm d} t} \right|_{\rm VS,m}, \end{eqnarray}$$

(D10)

$$\begin{eqnarray} \frac{{\rm d} \mathrm{ i}^2}{{\rm d} t} &=& \left. \frac{{\rm d} \mathrm{ i}^2}{{\rm d} t} \right|_{\rm drag} +\left. \frac{{\rm d} \mathrm{ i}^2}{{\rm d} t} \right|_{\rm VS,M} +\left. \frac{{\rm d} \mathrm{ i}^2}{{\rm d} t} \right|_{\rm VS,m}. \end{eqnarray}$$

(D11)

The gas damping rates are given by (Adachi, Hayashi & Nakazawa 1976; Inaba et al. 2001),

$$\begin{eqnarray} \frac{{\rm d} \mathrm{ e}^2}{{\rm d} t} &=& - \frac{2 \mathrm{ e}^2}{\tau _{\rm aero,0}} \left(\frac{9}{4} \eta ^2 +\frac{9}{4\pi } \zeta ^2 \mathrm{ e}^2 +\frac{1}{\pi } \mathrm{ i}^2 \right)^{1/2}, \end{eqnarray}$$

(D12)

$$\begin{eqnarray} \frac{{\rm d} \mathrm{ i}^2}{{\rm d} t} &=& - \frac{ \mathrm{ i}^2}{\tau _{\rm aero,0}} \left(\eta ^2 +\frac{1}{ \pi } \zeta ^2 \mathrm{ e}^2 +\frac{4}{\pi } \mathrm{ i}^2 \right)^{1/2}, \end{eqnarray}$$

(D13)

where ζ ∼ 1.211 and τ_{aero, 0} is given by,

$$\begin{eqnarray} \tau _{\rm aero,0} &=& \frac{2 m_{\rm pl}}{ C_{\rm d} \pi R_{\rm pl}^2 \rho _{\rm gas} v_{\rm K}}, \end{eqnarray}$$

(D14)

where m_pl is the mass of a planetesimal, C_d is the drag coefficient and is set to 1, R_pl is the radius of a planetesimal, ρ_gas is the density of the disc gas, and v_K is the Kepler velocity. We assume the vertical isothermal disc and use the mid-plane density for ρ_gas. The excitation rates of mean square orbital eccentricities and inclinations by the protoplanet are given by (Ohtsuki, Stewart & Ida 2002),

$$\begin{eqnarray} \left. \frac{{\rm d} \mathrm{ e}^2}{{\rm d} t} \right|_{\rm VS,M} &=& \left(\frac{M_{\rm p}}{3 b M_{*} P_{\rm orb}} \right) P_{\rm VS}, \end{eqnarray}$$

(D15)

$$\begin{eqnarray} \left. \frac{{\rm d} i^2}{{\rm d} t} \right|_{\rm VS,M} &=& \left(\frac{M_{\rm p}}{3 b M_{*} P_{\rm orb}} \right) Q_{\rm VS}, \end{eqnarray}$$

(D16)

where b is the full width of the feeding zone and is set to 10, M_* is the central star’s mass, and P_VS and Q_VS are given by,

$$\begin{eqnarray} P_{\rm VS} &=& \frac{73 {\tilde{e}}^2}{10 \Lambda ^2} \ln \left(1 +10\frac{\Lambda ^2}{{\tilde{e}}^2} \right) + \frac{72 I_{\rm PVS} (\beta)}{\pi \tilde{e}\tilde{i}} \ln \left(1 +\Lambda ^2 \right), \end{eqnarray}$$

(D17)

$$\begin{eqnarray} Q_{\rm VS} &=& \frac{4 {\tilde{i}}^2 +0.2 \tilde{i}{\tilde{e}}^3 }{10 \Lambda ^2 \tilde{e}} \ln \left(1 +10 \Lambda ^2 {\tilde{e}} \right) + \frac{72 I_{\rm QVS} (\beta)}{\pi \tilde{e}\tilde{i}} \ln \left(1 +\Lambda ^2 \right), \end{eqnarray}$$

(D18)

where |$\Lambda =\tilde{i} ({\tilde{e}}^2+{\tilde{i}}^2)/12$|⁠. For 0 < β ≤ 1, I_PVS and I_QVS can be approximated by (Chambers 2006),

$$\begin{eqnarray} I_{\rm PVS} (\beta) \simeq \frac{\beta -0.36251}{0.061547 +0.16112 \beta +0.054473 \beta ^2}, \end{eqnarray}$$

(D19)

$$\begin{eqnarray} I_{\rm QVS} (\beta) \simeq \frac{0.71946 -\beta }{0.21239 +0.49764 \beta +0.14369 \beta ^2}. \end{eqnarray}$$

(D20)

We also consider the excitation rates of mean square orbital eccentricities and inclinations by the mutual interactions of planetesimals are given by (Ohtsuki et al. 2002),

$$\begin{eqnarray} \left. \frac{{\rm d} \mathrm{ e}^2}{{\rm d} t} \right|_{\rm VS,m} &=& \frac{1}{6} \sqrt{\frac{G a_p}{M_{*}} } \Sigma _\mathrm{pl} h_\mathrm{m} P_{\rm VS}, \end{eqnarray}$$

(D21)

$$\begin{eqnarray} \left. \frac{{\rm d} \mathrm{ i}^2}{{\rm d} t} \right|_{\rm VS,m} &=& \frac{1}{6} \sqrt{\frac{G a_p}{M_{*}} } \Sigma _\mathrm{pl} h_\mathrm{m} Q_{\rm VS}, \end{eqnarray}$$

(D22)

with

$$\begin{eqnarray} h_\mathrm{m} = \left(\frac{2 m_\mathrm{pl}}{3 M_{*}} \right)^{1/3}. \end{eqnarray}$$

(D23)

APPENDIX E: SUCCESSFUL ANALOGUES WHEN THE ALLOWED MASS DIFFERENCE IS DOUBLED

In the main text, we identified successful Uranus and Neptune analogues by finding all simulations where the total mass of the inner and outer planet at the time of disc dissipation was within 1.5 of 14.5 and |$17.1\, \textrm {M}_{{\oplus }}$|⁠, respectively. Then we removed all simulations where the H–He mass fraction of any of the planets was above 20 per cent at disc dissipation. Furthermore, we also allowed for the planets to switch places after formation. In other words, we also searched for simulations where the mass of the inner planet was within |$1.5\, \textrm {M}_{{\oplus }}$| of the current mass of Neptune, and the mass of the outer planet was within |$1.5\, \textrm {M}_{{\oplus }}$| of the current mass of Uranus. Figs E1 and E2 show the number of successful analogues that were obtained when we increased the allowed mass difference to |$3\, \textrm {M}_{{\oplus }}$|⁠.

Figure E1.

Same as Fig. 5, but the allowed mass difference was increased to |$3\, \textrm {M}_{{\oplus }}$|⁠, resulting in more successful analogues.

Figure E2.

Same as Fig. 7, but the allowed mass difference was increased to |$3\, \textrm {M}_{{\oplus }}$|⁠, resulting in more successful analogues.

In Fig. E1 where the planetary embryos form and begin to accrete mass at the same time, increasing the allowed mass difference results in many more successful analogues while using the constant model. However, when using the more realistic evolving model, there is almost no change. There are two successful analogues found in the case with no planetesimal accretion, and none when there is planetesimal accretion. When the outer embryo is formed with a higher mass (Fig. E2), increasing the allowed mass difference more than doubles the amount of successful analogues found in both the constant model and the evolving model.