© ESO 2020

1 Introduction

Jupiter’s formation has been studied extensively since the very early works on planet formation (e.g. Bodenheimer & Pollack 1986; Pollack et al. 1996), and yet, there are still many open questions regarding its origin. In the standard model for giant planet formation, core accretion, the formation of a giant planet begins with the accretion of solids followed by the accretion of gas disc (see Helled et al. 2014 for review and references therein). The timescale for core formation depends on the size of the assumed solids. Classical models of core accretion assumed these solids to be large planetesimals with sizes of 100 km (Pollack et al. 1996; Alibert et al. 2005a), while more recent studies point out the importance of growing a core by cm-size pebbles (Lambrechts & Johansen 2014). In this work, we use the terms “solids” and “heavy elements” as synonyms, despite the physical state of such compounds after being accreted by the protoplanet. As the protoplanet reaches crossover mass (where the mass of the H–He envelope equals the mass of heavy elements), the gas accretion rate continuously increases and exceeds the heavy-element accretion rate. Soon after crossover time, the disc can no longer supply the gas required by the cooling and strong self-gravity of the protoplanet. At this point, the gas accretion rate is limited by the availability of the disc to supply gas, the protoplanet becomes detached from the disc, and the gas accretion rate is dictated by the disc’s hydrodynamics (Tanigawa & Watanabe 2002; Lissauer et al. 2009; Mordasini 2013).

It is still unknown whether giant planet formation is dominated by the accretion of small (pebbles) or large (planetesimals) solids. Typically, the accretion of planetesimal requires a longer core formation timescale due to the moderate accretion rates (Fortier et al. 2007, 2013). On the other hand, pebbles undergo orbital decay through the interaction with the disc of gas, and they can be easily and effectively accreted by a growing core due to efficient gas friction within the Hill’s sphere (e.g. Ormel 2017). This gave to the pebble accretion scenario a special attractiveness to solve the historical problem of planetesimal accretion being too slow to grow a critical core in timescales comparable to the disc’s lifetime (Lambrechts & Johansen 2014; Lambrechts et al. 2014).

When a planet grows by accreting pebbles, at certain point it reaches the so-called pebble isolation mass. This is the necessary mass to perturb the disc and create a pressure bump beyond the orbit of the protoplanet. Within a pressure bump, the sign of the pressure gradient is reversed, so the direction of the drift of pebbles changes accordingly. Thus, a pressure bump acts as a “trap” of pebbles, and pebbles cannot be accreted by the protoplanet any longer. Hence, once the protoplanet attains the pebble isolation mass, pebble accretion is halted. The pebble isolation mass is estimated to be ~ 20 M_⊕ and its exact value depends on the assumed disc properties (e.g. Lambrechts et al. 2014; Ataiee et al. 2018; Bitsch et al. 2018).

Recent analysis of meteoritic data of isotopic anomalies found in iron meteorites were used by Kruijer et al. (2017, hereafter K17) to constrain Jupiter’s formation chronology. According to K17, the isotopic anomalies indicate that carbonaceous and non-carbonaceous chondrites were accreted isolated from each other between times ~1 and ~3–4 Myr after formation of calcium-aluminium-rich inclusions (CAIs). K17 proposes that proto-Jupiter reached pebble isolation mass at time 1 Myr, and planetesimals were formed at both sides of proto-Jupiter, until time ~ 2 Myr for the non-carbonaceous chondrites and until time ~3–4 Myr for the carbonaceous chondrites. At that time, Jupiter became massive enough (~ 50 M_⊕) to scatter planetesimals, causing both reservoirs to reconnect, and to match in this way, the timescale inferred by the meteoritic record.

Alibert et al. (2018, hereafter A18), built upon the interpretation of K17, and investigated in detail the validity of this formation scenario, in which between times 1 and 3 Myr from the formation of CAIs, proto-Jupiter separates the solids. They found that in order to fulfil this time constraint, a solid accretion rate of ~10⁻⁶−10⁻⁵M_⊕ yr⁻¹ is needed to prevent the onset of runaway gas accretion before time ~3 Myr. Since the protoplanet reached the pebble isolation mass at time ~1 Myr, after this time pebbles could not be accreted by the planet, and hence, the heat source during that period must come from the accretion of planetesimals.

The exact formation process of Jupiter is not fully constrained and it is unknown whether the early stages of the planet’s formation are dominated by pebble accretion, planetesimal accretion, or both, as recently suggested by A18. In the hybrid formation scenario of A18, there are three formations stages. In a first phase, rapid pebble accretion supplies the major part of Jupiter’s core mass. In this stage planetesimal accretion is negligible, and large primordial planetesimals are thought to be excited by the growing planet reaching high collision velocities that lead to destructive collisions, and the creation of small, second-generation planetesimals. In a second phase, the accretion of (small) planetesimals provides enough energy to hinder runaway gas accretion for ~2 Myr. In this stage Jupiter is massive enough to prevent pebble accretion. Lastly, there is a third phase of runaway gas accretion, as usual in core accretion models (Bodenheimer & Pollack 1986; Pollack et al. 1996; Lissauer et al. 2009). In this scenario, phase-2 and phase-3 are somewhat similar to the ones in the traditional planetesimal accretion scenario (Pollack et al. 1996), although phase 2 is defined differently since it starts at the pebble isolation mass and not at the planetesimal isolation mass. The role of planetesimal and pebble accretion in giant planet formation, and therefore the characteristics of the different formation phases, is still being investigated.

Several important questions regarding Jupiter’s formation remain unanswered, for example: where did Jupiter’s formation start? How many heavy elements can be accreted by the growing planet, and where did they come from? What is the expected structure and global water abundance in Jupiter? Unfortunately, these questions cannot be answered easily. However, the combination of new data and new theoretical models allows us to better understand Jupiter’s origin. In this work, we investigate Jupiter’s growth history accounting for different planetesimal sizes and accretion rates. In particular, we focus on heavy-element accretion during the last stage of formation and connect it to Jupiter’s structure models as well as to the mass-metallicity relation of giant exoplanets.

2 Methods

We investigate Jupiter’s in situ formation in the framework of the core accretion model. As in A18, we assume that a certain core is already formed by pebble accretion at time = 1 Myr of the disc evolution. That is the approximate time at which proto-Jupiter reaches the pebble isolation mass and pebble accretion stops (K17, A18). Our simulations begin at that point, and we assume that from that time on the accreted heavy elements are in the form of large planetesimals. During the early stages of Jupiter’s growth, the protoplanet is embedded in the protoplanetary disc, and no clear boundary delimits the protoplanet from the disc. This phase is known as the attached phase. However, once the mass of the envelope surpasses the mass of the core, the gas accretion onto the protoplanet becomes so large, that the disc can no longer supply the amount of gas required by the contraction of the protoplanet. At this stage gas accretion is limited by the disc’s supply (the detached phase).

In addition to the attached and detached phases, we refer to the following three phases of Jupiter’s growth throughout the paper, in analogy to A18:

Phase 1corresponds to the growth of the core dominated by pebble accretion.

Phase 2starts when the pebble isolation mass is reached, and finishes at the crossover mass. During this phase the solid accretion is dominated by large planetesimals (10–100 km) (Sect. 2.1.1). The protoplanet accretes gas concurrently, at a rate that is dictated by the cooling of the envelope or its Kelvin–Helmholtz timescale (Sect. 2.3).

Phase 3 starts at crossover mass and finishes when the mass of Jupiter is attained. During this phase gas is accreted in a runaway fashion. Soon after crossover time, the detached phase (Mordasini 2013), or, equivalently, “the disc-limited accretion” phase (Shiraishi & Ida 2008) begins (see Sect. 2.4).

2.1 Planetesimal accretion

The initial size of planetesimals in the Solar Nebula is still a matter of debate (Kenyon & Bromley 2012; Weidenschilling 2011; Morbidelli et al. 2009). Simulations based on the streaming instability model, the most accepted scenario to explain planetesimal formation, predict initial sizes in the range of ~ 10–100 km (e.g. Simon et al. 2017). Therefore, in this work we use the more accepted planetesimal size of r_p = 100 km for the nominal model. This is in contrast to A18 where small planetesimals were required. We also show results for r_p = 10 km in Sect. 3.1.2.

2.1.1 Oligarchic growth

We implemented the planetesimal accretion rate given by Guilera et al. (2010) and Fortier et al. (2013), and as originally presented by Inaba et al. (2001); which is hereafter referred to as Oligarchic Growth (OG). This model includes the gravitational stirring on the planetesimals due to the presence of the protoplanet and the damping of eccentricities by gas drag. An equilibrium state between the stirring and the damping is assumed, which is appropriate for the case of large planetesimals as assumed here (Chambers 2006).

Under equilibrium, the root mean square eccentricity of planetesimals (e) is given by (Thommes et al. 2003): $$e=1.7{(\frac{m_{\text{p}}^{1/3}{\rho }_{\text{p}}^{2/3}}{b{C}_D{\rho }_{\text{gas}}{a}_{\text{p}}})}^{1/5}{(\frac{M_{\text{P}}}{M_{*}})}^{1/3},$$(1)

where m_p is the planetesimal’s mass, ρ_p its density, b a constant equal to 10, C_D the non-dimensional drag coefficient taken as 0.7, and ρ_gas the gas density. In the equilibrium approximation, i = e∕2.

The planetesimal accretion rate is given by (Chambers 2006): $${\dot{M}}_{\text{Z}}=\Omega \Sigma R_{\text{H}}^2{P}_{\text{coll}},$$(2)

where Ω is the Keplerian frequency, Σ is the surface density of planetesimals, R_H the Hill radius, and P_coll the collision probability between the embryo and the planetesimals.

The surface density of planetesimals evolves with time as a result of accretion onto the protoplanet and scattering. The protoplanet accretes planetesimals from a feeding zone of half-width = 5 R_H, and ejects planetesimals at a rate given by (Ida & Lin 2004): $$\frac{\text{accretion\hspace{0.17em}rate}}{\text{ejection\hspace{0.17em}rate}}={(\frac{v_{\text{esc},\text{disc}}}{v_{\text{surf},\text{planet}}})}^4,$$(3)

where $v_{\text{esc},\text{disc}}=\sqrt{2G{M}_{\odot }/a}$ is the escape velocity from the central star at the planet’s position, and $v_{\text{surf},\text{planet}}=\sqrt{G{M}_{\text{P}}/R_{\text{cap}}}$ is the planet’s characteristic surface speed; R_cap being the capture radius.

The collision probability onto the protoplanet, P_coll depends on the excitation between the embryo and planetesimals, and is given by: $$P_{\text{coll}}=\text{min}(P_{\text{med}},{(P_{\text{high}}^{-2}+P_{\text{low}}^{-2})}^{-1/2}),$$(4)

where $$\begin{array}{ll}{P}_{\text{high}}=\hfill & \frac{{(R_{\text{cap}}+r_{\text{p}})}^2}{2\pi R_{\text{H}}^2}(I_F(\beta )+\frac{6R_{\text{H}}{I}_{\text{G}}(\beta )}{(R_{\text{cap}}+r_{\text{p}}){\tilde{e}}^2}),\hfill \\ P_{\text{med}}=\hfill & \frac{{(R_{\text{cap}}+r_{\text{p}})}^2}{4\pi R_{\text{H}}^2\tilde{i}}(17.3+\frac{232R_{\text{H}}}{R_{\text{cap}}+r_{\text{p}}}),\hfill \\ P_{\text{low}}=\hfill & 11.3{(\frac{R_{\text{cap}}+r_{\text{p}}}{R_{\text{H}}})}^{1/2},\hfill \end{array}$$

where ẽ = ae∕R_H, ĩ =ai∕R_H, and β = ĩ∕ẽ. R_cap is the capture radius of the protoplanet, i.e, the effective cross section of collision when a core is surrounded by a gaseous atmosphere. The captureradius reduces to the core radius in the limit of no atmosphere. For the computation of the capture radius, we implement the prescription given by Inaba & Ikoma (2003).

The functions I_F(β) and I_G(β) are approximately given by (Chambers 2006): $$\begin{array}{ll}{I}_{\text{F}}(\beta )=\hfill & \frac{1+0.95925\beta +0.77251{\beta }^2}{\beta (0.13142+0.12295\beta )},\hfill \\ I_{\text{G}}(\beta )=\hfill & \frac{1+0.3996\beta }{\beta (0.0369+0.048333\beta +0.006874{\beta }^2)}.\hfill \end{array}$$

2.1.2 Detached phase

For our nominal model, we use Eq. (2) during all the growth of Jupiter. However, that accretion rate of planetesimals was calculated for Earth-mass embryos embedded in a gaseous disc. When the planet’s growth is dominated by the accretion of gas, Shiraishi & Ida (2008, hereafter SI08) showed, from analytical fits to N-body simulations, that the accretion rate of planetesimals is given by two equations, depending on the competition between the expansion of the feeding zone and the damping of eccentricity of planetesimals within it.

When the expansion of the feeding zone dominates, planetesimals enter continuously to the protoplanet’s feeding zone, so no planetesimal gap is created. Shiraishi & Ida (2008) express this condition in terms of η = v_H∕v_damp, where v_H is the rate of expansion of the feeding zone and v_damp the rate of damping of eccentricity by gas drag. The non-gap case fulfils η > 1, and in this case the accretion rate of solids is given by (Eqs. (17) and (24) of Shiraishi & Ida 2008): $${\dot{M}}_{\text{Z},\text{nogap}}={10}^{-6}{\rho }^{1/2}{(\frac{R_{\text{P}}}{R_{\oplus }})}^2{f}_{\text{d}}{(\frac{v_{\text{H}}}{v_{\text{scat}}})}^{0.8}{M}_{\oplus }{\text{yr}}^{-1},$$(10)

where ρ is the mean density of the protoplanet, R_P is the physical radius of theprotoplanet, f_d is the surface density of solids in Minimum Mass Solar Nebula (MMSN) units, and v_scat is the rate at which planetesimals are scattered from the feeding zone. The ratio of the speed of Hill’s radius expansion and planetesimal scattering follows (Shiraishi & Ida 2008, their Eq. (22)): $$\frac{v_{\text{H}}}{v_{\text{scat}}}=4.1{(\frac{a}{5\text{\hspace{0.17em}au}})}^{3/2}{(\frac{M_{\text{P}}}{M_{\oplus }})}^{-1/3}{(\frac{{\tau }_{\text{g}}}{{10}^4\text{\hspace{0.17em}yr}})}^{-1},$$(11)

where τ_g is the timescale to accrete gas, given by τ_g = M_P∕Ṁ_HHe.

When η < 1, a planetesimal gap is formed and the accretion rate of solids follows (Eqs. (17) and (25) of Shiraishi & Ida 2008): $${\dot{M}}_{\text{Z},\text{gap}}={10}^{-6}{\rho }^{1/2}{(\frac{R_{\text{P}}}{R_{\oplus }})}^2{f}_{\text{d}}{(\frac{v_{\text{H}}}{v_{\text{damp}}})}^{1.4}{M}_{\oplus }{\text{yr}}^{-1},$$(12)

where: $$\frac{v_{\text{H}}}{v_{\text{damp}}}=0.8{(\frac{a}{5\text{\hspace{0.17em}au}})}^{3/4}{(\frac{M_{\text{P}}}{M_{\oplus }})}^{-1/6}{(\frac{{\tau }_{\text{g}}}{{10}^4\text{\hspace{0.17em}yr}})}^{-1}{(\frac{{\tau }_{\text{damp}}}{{10}^4\text{\hspace{0.17em}yr}})}^{1/2},$$(13)

τ_damp is the timescale to damp the planetesimals’ eccentricity, and corresponds approximately to τ_damp ~ 10⁶ yr for large planetesimals as the ones considered in this work (Shiraishi & Ida 2008). Throughout this work, we refer to the implementation of Eqs. (10)–(13) as the SI08 prescription. It is important to remark that Eqs. (10) and (12) are the general form of the fit found by SI08, and not the particular case of ρ = 1 g cm⁻³ assumed in Eqs. (24) and (25) of SI08. We test the effect of these accretion rates during the runaway gas phase, when the planet is detached from the disc. The results are presented in Sect. 3.4.

The above solid accretion rates depend on the planetary radius as $~R_{\text{P}}^{1/2}$. 1D simulations have shown that during the detached phase the planetary radius contracts very fast to ~ 2 Jupiter radius, and remains fairly constant during the remaining of gas accretion (Lissauer et al. 2009; Mordasini 2013). Following those studies, we adopt R_P = 1.8 R_Jup when implementing Eqs. (10) and (12).

2.2 Disc model

Initially, giant planets grow while being embedded in a protoplanetary disc. The disc temperature and pressure affect the outer boundary of the protoplanet, and hence, a disc model is required to study the planetary growth in a consistent manner. We adopt the disc model from Bitsch et al. (2015, hereafter B15). This disc evolution model adopts the accretion rate Ṁ onto the central star derived from observations (Hartmann et al. 1998): $$\log \left(\frac{\dot{M}}{M_{\odot }{\text{yr}}^{-1}}\right)=-8.00-1.40\times \log \left(\frac{t+{10}^5\text{\hspace{0.17em}yr}}{{10}^6\text{\hspace{0.17em}yr}}\right),$$(14)

where t is the evolution time, and the error bars from the original expression in Hartmann et al. (1998) are neglected.

Ṁ is related to the disc’s viscosity and the gas surface density Σ_gas via: $$\dot{M}=3\pi \nu {\Sigma }_{\text{gas}}=3\pi \alpha H^2\Omega {\Sigma }_{\text{gas}},$$(15)

where H is the disc’s scale height, and Ω the Keplerian frequency. α represents the viscosity and only acts as a heating parameter. B15 adopted α = 0.0054 to perform the radiative-hydrodynamical simulations, so this is the value that must be assumed for consistency. After computing the temperatureprofile from the formulas presented in the appendix of B15, the disc scale height can be obtained from: $$T={(\frac{H}{r})}^2\frac{G{M}_{*}}{r}\frac{\mu }{R}.$$(16)

Finally, Σ_gas can be calculated from Eq. (15).

The initial surface density of planetesimals (at time = 1 Myr), Σ₁, is a free parameter of the model. It is important to note that this quantity is typically a free parameter in planet formation simulations, and is hidden in the assumed initial disc dust-to-gas ratio which is often set to ~0.01 guided by observations. However, studies of dust growth and drift (e.g, Birnstiel et al. 2012) suggest that the initial dust-to-gas ratio can vary significantly from the planetesimal-to-gas ratio in a given location of the disc (Dra̧żkowska et al. 2016). In addition, planetesimals can be formed in regions where dust is piled up, and the dust-to-gas ratio reaches values of ~1. In fact, the planetesimal-to-gas ratio could be even larger than 1 at t = 1 Myr of the disc’s evolution, a few 10⁵ yr after the onset of planetesimal formation (Dra̧żkowska et al. 2016, from their Figs. 2 and 4). Since various theoretical models and observations reveal an unprecedented complexity in the dust distribution within protoplanetary discs, the simplification of assuming an initial planetesimal-to-gas ratio equal to a fixed dust-to-gas ratio is no longer justified and is clearly an over-simplification.

2.3 Gas accretion: attached phase

We simulate the planetary growth by solving the standard planetary internal structure equations, assuming the luminosity results from the accretion of solids and envelope contraction (see Alibert et al. 2013; Venturini et al. 2016, for details): $$\begin{array}{lll}\frac{\text{d}{r}^3}{\text{d}m}\hfill & =\hfill & \frac{3}{4\pi \rho },\hfill \\ \frac{\text{d}P}{\text{d}m}\hfill & =\hfill & \frac{-G(m+M_{\text{core}})}{4\pi r^4},\hfill \\ \frac{\text{d}T}{\text{d}P}\hfill & =\hfill & \text{min}({\nabla }_{\text{conv}},{\nabla }_{\text{rad}}),\hfill \end{array}$$

where G is the gravitational constant, and r, P, T, ρ are the radius, the pressure, temperature and density inside the envelope, respectively.

The density is given as a function of T and P by the equation of state (EOS). For the nominal case, we assume that the solids reach the core, meaning that the envelope keeps a H–He composition. For this case we use the EOS of Saumon et al. (1995). Since we assume that large planetesimals are accreted, it is expected that despite of thermal ablation, the material is deposited deep in the envelope, especially the refractory materials at early formation stages (e.g. Valletta & Helled 2019). We also consider an enriched-envelope case in which we assume that all the volatiles remain in the envelope and mix uniformly throughout it (Venturini et al. 2016). Although this case is probably more realistic for pebble accretion, we still evaluate this possibility in 3.1.3. We suggest that future formation models should couple the heavy element deposition with the resulting compositional gradient and its effect on the internal structure.

The temperature gradient is given by either the radiative gradient (∇_rad): $${\nabla }_{\text{rad}}=\frac{3\kappa L}{64\pi \sigma G(m+M_{\text{core}})T^3}$$(20)

or the convective gradient, equal to the adiabatic one. In these formulas, σ is the Stefan–Boltzmann constant, κ is the opacity, and M_core is the mass of the core.

For the opacity, we take κ = max(f κ_dust, κ_gas), where the dust opacity is taken from Bell & Lin (1994), the gas from Freedman et al. (2014), and for the nominal model we take f = 1. For a solar composition this definition is equivalent to consider the Bell & Lin (1994) opacities throughout the envelope, but with the current definition we can scale down the dust opacity to consider grain growth. Moreover, the Freedman et al. (2014) gas opacities account for a large range of envelope metallicities, so they are better suited for the cases where envelope enrichment is considered during growth.

The internal structure equations are solved, using as boundary conditions the pressure and temperature in the protoplanetary disc at the position of the planetary embryo. The planetary radius as a combination of the Hill and Bondi radii (Lissauer et al. 2009): $$R_{\text{P}}=\frac{G{M}_{\text{P}}}{\left(C_{\text{S}}^2+4G{M}_{\text{P}}/R_{\text{H}}\right)}$$(21)

where $C_{\text{S}}^2$ is the square of the sound velocity in the protoplanetary disc at the planet’s location a, and the Hill radius $R_{\text{H}}=a{\left(\frac{M_{\text{P}}}{3M_{\odot }}\right)}^{1/3}$.

2.4 Gas accretion: detached phase

Once the gas accretion given by the gas disc dynamics is smaller than the gas required from the protoplanet by its own cooling, the protoplanet continues accreting gas at a rate determined by the disc flow (Ṁ_{gas, disc}). In other words, the gas accretion onto the protoplanet is given by: $${\dot{M}}_{\text{gas}}=\text{min\hspace{0.17em}}\{{\dot{M}}_{\text{gas},\text{KH}},{\dot{M}}_{\text{gas},\text{disc}}\}$$(22)

where Ṁ_{gas, KH} is the gas accretion obtained from computing the cooling of the protoplanet (KH standing for Kelvin–Helmholtz), as explained in Sect. 2.3, and Ṁ_{gas, disc} is the accretion rate of gas that the disc is able to supply to the protoplanet. For the latter, we take Eq. (3) from Lissauer et al. (2009, hereafter L09): $$\log (\frac{{\dot{M}}_{\text{gas},\text{disc}}}{{\Sigma }_{\text{gas}}{a}^2/P})=c_0+c_1\log (\frac{M_{\text{P}}}{M_{*}})+c_2{\log }^2(\frac{M_{\text{P}}}{M_{*}}),$$(23)

where P being the planet’s orbital period, c₀ = −18.67, c₁ = −8.97 and c₂ = −1.23. This expression was obtained as a fit to 3D hydrodynamical simulations for α =4 × 10⁻³. Following also L09, we multiply the gas accretion rate by a linear function equal to 1 at 0.85 Jupiter mass and equal to 0 at 1 Jupiter mass (1 M_J), in order to terminate gas accretion when reaching one Jupiter mass. It should be noted that in our simulations we do not consider the accretion shock during the last phase of rapid gas accretion which can lead to entropy gradients inside the planet and persist for up to 10 Myr after its formation (Cumming et al. 2018). While this process is highly important for the determination of Jupiter’s primordial entropy and heat transport mechanism at early ages, it is less likely to affect the heavy-element accretion during the last stage of formation. This, however, should be investigated in detail in future studies.

A summary of the main recipes and parameters used in this work is shown in Table 1.

3 Results

3.1 Formation at a = 5 au

To illustratethe output of the model, we first show, in Fig. 1, an example of planetary growth at a = 5 au that matches the meteoritic constraints. That is, Jupiter attains pebble isolation mass at time 1 Myr and M_P = 50 M_⊕ at a time of 3–4 Myr of disc evolution (K17). For this particular case, the initial core mass is M_core,1 = 14 M_⊕, and the initial planetesimal surface density is Σ₁ = 7 g cm⁻². The planetesimal size assumed is r_p = 100 km, and the planetesimals’ density, ρ_p = 1 g cm⁻³. The solid lines indicate the growth within the attached phase, and the dotted lines within the detached phase. The crossover mass is M_cross = 21.9 M_⊕ with the crossover time being t_cross = 3.7 Myr. The detached phase starts when the total planetary mass is M_P = 82 M_⊕, the mass of heavy elements is M_Z = 22.3 M_⊕, and the accretion rate of H–He is Ṁ_HHe = 9.7 × 10⁻³M_⊕ yr⁻¹. During the detached phase, additional 0.6 M_⊕ of heavy elements are accreted (OG prescription).

3.1.1 Meteoritic constraints and comparison with A18

Figure 2 summarises the results for all simulations at a = 5 au with planetesimal sizes of 100 km. We refer to these set of simulations as the “nominal case’. Each point in each panel is the growth of one planet, characterised by an initial core mass (M_core,1, the core mass att = 1 Myr), and planetesimal surface density at 1 Myr of disc evolution (Σ₁). The left panel plots the time required to reach M_scat, the mean mass the protoplanet must attain to scatter planetesimals into orbits inside proto-Jupiter, which range between 30 and 70 M_⊕ (Alibert et al. 2018). The green circles represent the models that fulfil the meteoritic constraint of reaching t_scat ≈ 3−4 Myr after CAI formation. The black circles indicate a too fast formation timescale, of t_scat < 2.5 Myr; and the opposite for grey ones, with t_scat > 4.5 Myr. Models in black in the upper left corner, above the red line, correspond to those whose M_core,1 is supercritical, and therefore, the growth cannot be computed. These cases would yield too short formation timescales.

Regarding the trends of the plot, for a fixed Σ₁, the larger the initial core mass, the less time is needed to reach the crossover mass, so the overall formation timescale shortens when we increase M_core,1. When we move in the horizontal direction, the larger the initial planetesimal surface density, the larger the accretion rate of solids, so the formation timescale also shortens in this direction. We find that when accreting planetesimals of 100 km size at a = 5 au, the pebble isolation mass must range between 6 and 16 M_⊕ to fulfil the K17 constraints.

At first sight, the values considered for Σ₁ might look too large. However, in terms of planetesimal-to-gas ratio, Σ₁ = 25 g cm⁻² corresponds to Σ₁∕Σ_gas = 0.231 in our model, value expected in regions where planetesimals are formed (Dra̧żkowska et al. 2016). The ranges of Σ₁ with its corresponding planetesimal-to-gas ratio and equivalent in units of MMSN are listed in Table 2. The right panel of the figure shows the mean accretion rate of solids during Phase-2 (between time 1 Myr and crossover time). We corroborate that the values range between Ṁ_Z ≈ 10⁻⁶−10⁻⁵M_⊕ yr⁻¹, as found by A18.

The timescale for Jupiter’s formation provided by the meteoritic constraints can be matched when the second stage of Jupiter’s growth is dominated by the accretion of 100 km planetesimals. This was not found by the simple estimates of accretion rates given in A18. The difference is that in this work, we follow the planet’s growth consistently, accounting for the enhancement factor of the accretion due to the presence of an atmosphere. Indeed, the capture radius can be larger than the core radius by a factor of ~ 20. This is shown in detail in Appendix A. A large capture radius increases the planetesimal accretion rate (Eqs. (5)–(7)), relaxing the implications on the minimum size of planetesimals that match the K17 scenario. In summary, we find that the “meteoritic time constraints” can be easily matched by accreting large planetesimals, which are presumably the most abundant type when planets are forming (e.g. Johansen et al. 2014).

Fig. 1 Example of planetary growth, where the initial surface density of planetesimals is Σ1 = 7 g cm−2 and the initial core mass is Mcore,1 = 14 M⊕. MZ is the mass of heavy elements, MHHe the mass of H and He, and MP the total mass of the planet. The solid lines indicate that the protoplanet is in the attached phase, whereas the dotted lines show the protoplanet in the detached phase. Bottom panel: simple zoom to appreciate better the growth of MZ.

Fig. 2 Grid of simulations for a = 5 au and rp = 100 km. Left: mean time to reach the scattering mass as a function of the initial surface density of planetesimals (Σ1) and the initial core mass (Mcore,1). The light green colour indicates the cases compatible with the meteoritic record. For grey circles the growth is too slow (longer than 4.5 Myr), and for black ones, too fast (shorter than 2.5 Myr). Points in the upper left corner above the red line are supercritical and their structure cannot be computed. Points in the lower left corner below the purple line correspond to cases whose growth is so slow that the mass of Jupiter is not reached within 10 Myr of disc evolution. Right: mean accretion rate of solids during Phase 2 as a function of the initial surface density of planetesimals. The initial core mass is shown in the colour-bar. Circles with black border represent cases where the meteoritic constraints are fulfilled (green band of left panel).

3.1.2 Heavy elements in Jupiter: comparison with Jupiter’s structure models

The left panel of Fig. 3 shows the total amount of heavy elements when the protoplanet reaches 1 Jupiter mass (M_Jup). We find that Jupiter could have accreted between 18 and 34 M_⊕ of heavy elements during its formation when implementing the OG prescription; and between 20 and 43 M_⊕ when implementing the prescription of SI08. These findings are in good agreement with estimates from recent structure models of Jupiter thatfit the new gravity data measured by the Juno spacecraft (Wahl et al. 2017) (green shaded area in the left panel). We note that the total heavy-element mass accreted until reaching Jupiter’s mass is strongly dependent on Σ₁.

The amount of heavy elements accreted during the final phase of Jupiter’s formation (the detached phase) is shown in right panel of Fig. 3. The range of heavy elements accreted in this fast stage varies between 0.3 and 14 M_⊕. It should be noted that our results on Jupiter’s enrichment during the different growth phases, are independent of the timescale constraint for Jupiter’s formation as presented by K17. At the same time, our estimates for the total heavy-element mass in the planet rely on prescriptions presented by other studies. Determining the heavy-element mass accreted during Jupiter’s formation is challenging and is still a topic of intensive investigation. A more robust estimate would be possible only when the different approaches and methods (i.e. quasi-hydrostatic growth, hydrodynamics, N-body interactions) are combined in one theoretical framework. Finally, the heavy-element mass estimates from Jupiter structure models also suffer from uncertainties. Our study demonstrates how different assumed parameters affect the inferred composition. This can then be used to better understand the origin of Jupiter’s enrichment and for the interpretation of the formation paths of enriched giant exoplanets, as we discuss below.

We next repeat the simulations assuming a planetesimal size of r_p = 10 km. The results are shown in Fig. 4. Compared to the runs of r_p = 100 km, the growth is more efficient, meaning with smaller Σ₁, the same growth timescale than with r_p = 100 km can be attained. This is a typical feature of oligarchic growth: smaller planetesimals are more easily accreted (e.g, Fortier et al. 2007; Guilera et al. 2010).

Interestingly, the width of the solutions that fit the scenario of K17 (green band of top left panel of Fig. 4) is narrower in this case. Regarding the solids accreted with r_p = 10 km until the planet reaches Jupiter mass, we note that the values are very similar to the case of r_p = 100 km (Fig. 4, right panel),but the range of Σ₁ was assumed smaller for this set of simulations.

The fullsummary of all heavy elements accreted in each phase for all the models is shown in Appendix D.

Fig. 3 Simulations for a = 5 au and rp = 100 km with the same range in (Mcore,1, Σ1) as in Fig. 2. Orange circles correspond to the solid accretion given by SI08, purple circles to Oligarchic Growth (OG). The darker colours represent all runs that reach Jupiter mass within 10 Myr of disc evolution. The lighter colours highlight the cases that match the meteoritic constraints of K17. Left: total amount of heavy elements accreted until reaching Jupiter mass as a function of the initial surface density of planetesimals (Σ1). The green shaded area indicates the values of MZ compatible with estimates from Juno (Wahl et al. 2017). Right: amount of heavy elements accreted during the detached phase with the different solid accretion prescriptions.

Fig. 4 Grid of simulations for a = 5 au and rp = 10 km. Left: mean time to reach the scattering mass as a function of the initial surface density of planetesimals (Σ1) and initial core mass (Mcore,1). The light green colour indicates the cases compatible with the meteoritic record. For grey circles the growth is too slow, and for black ones, too fast. Right: total amount of heavy elements accreted until reaching Jupiter mass as a function of the initial surface density of planetesimals (Σ1). Transparentcircles: all results. Solid circles: cases compatible with the meteoritic constraints.The green shaded area indicates the values compatible with estimates from Juno (Wahl et al. 2017).

3.1.3 Dependence on opacity and envelope composition

We next study the dependence of the accreted heavy-element mass on the assumed envelope’s opacity. The dust opacity could be reduced compared to ISM values due to grain growth and settling (e.g. Movshovitz et al. 2010). Therefore, we first analyse results for the case where we reduce the dust opacity by a factor of 10 (i.e. f =0.1).

Figure 5 shows the total amount of heavy elements accreted during the attached phase as a function of the initial surface density of planetesimals, for the simulations that fulfil the meteoritic constraints. The squares show the results for f =0.1 while the circles represent the nominal run (f = 1). We note that for all runs, the amount of solids accreted in the low opacity case is smaller compared to the nominal case, by a shift of ~ 4 M_⊕. The fact that the protoplanet accretes less solids when attaining the same total mass is expected: the lower the opacity, the faster the gas accretion (Ikoma et al. 2000; Movshovitz & Podolak 2008). In other words, if the opacity is lower, a smaller core can hold a moremassive atmosphere in hydrostatic equilibrium such that the total mass of the planet is the same as in the high opacity case. Therefore, at Jupiter mass, for the low opacity case, less solids are expected to be accreted.

A very similar behaviour is observed when envelope enrichment is included. We run simulations for the case where we assume that half of the solids (the ices) remain homogeneously mixed in the envelope, while the refractories are assumed to sink to the core (same set up as in Venturini et al. 2016; Venturini & Helled 2017). Since in this case the timescale of formation is shorter, we find thatwhen envelope enrichment is considered, the amount of solids accreted by Jupiter is reduced compared to the nominal case (triangles in Fig. 5).

Overall, the total amount of solids accreted in the three cases is in the range of 14–32 M_⊕, and the difference between the nominal case and the other two cases are small, of the order of 2–6 M_⊕. Therefore, we conclude that the effect of envelope composition does not affect significantly the estimates of the heavy-element mass accreted by Jupiter.

Interestingly, when envelope enrichment is taken into account, the relation between the amount of heavy elements accreted during the attached phase and the initial planetesimal surface density of solids is flatter than in the nominal case (triangles in Fig. 5). This weaker dependence between Σ₁ and the amount of accreted heavy elements in the case when envelope enrichment is accounted for, is important, because it would indicate a weak relation between star metallicity and planet metallicity. Such a weak relation seems to be present in giant exoplanets, as is shown in Fig. 8 of Thorngren et al. (2016).

Fig. 5 Amount of heavy elements accreted during the attached phase (including the initial core formed by pebble accretion) as a function of the initial surface density of planetesimals for different values of dust opacities (circles and squares) and considering or not envelope enrichment (triangles and circles). All the cases that reach Jupiter mass within 10 Myr of disc evolution are shown. All cases correspond to planetesimal sizes of rp = 100 km. Red circles indicate the nominal simulations (f = 1), squared yellow points represent the case with dust opacity reduced by a factor of 10 (f = 0.1), and blue triangles indicate the case where envelope enrichment is included.

3.2 Dependence on the formation location

The exact location at which Jupiter spent most of its formation is not known and is a subject of ongoing research (e.g. Alibert et al. 2005b; Walsh et al. 2011). Recent work by Pirani et al. (2019) suggests that Jupiter should have migrated inwards, from a = 9−18 au until its current position, to explain an asymmetry in the number of Jupiter’s trojans. On the other hand, Desch et al. (2018) claim that a much closer formation location is required (even inside the water iceline for some fraction of Jupiter’s growth) to solve the CAI storage problem. Their model predicts the correct composition and formation location of 11 chondrite and 5 achondrite types.

We start this section by testing two additional formation locations for Jupiter, the first inside the water iceline at a = 1 au; and the second beyond the current location of Jupiter, at a = 10 au.

Although our models do not include planetary migration, we have considered various parameters that affect the planetary growth such as formation locations, initial core masses and initial surface density of planetesimals. As a result, this study covers a wide range of the parameter space that could emerge during migration. While planetary migration can affect the estimates for the accreted heavy-element mass, the migration history of Jupiter is poorly constrained. In addition, the processes and timescales associated with planetary migration are still being investigated and at the moment, choosing a specific prescription for it would make our results even more model dependent. As a result, in this study we rather focuson in situ formation, which is complex enough, and hope to include planetary migration in future research.

We keep assuming large planetesimals (i.e. r_p = 100 km), but in the case of the formation inside the water iceline, we set the planetesimal mean density to be 3.2 g cm⁻³, corresponding to a rocky composition (Fortier et al. 2013). The results are summarised in Fig. 6. When the formation location is set inside the iceline, a clear difference with the growth at a = 5 au is that the cases that fulfil the meteoritic time constraints are confined to an initial core mass of 16–17 M_⊕ (see Appendix B for details). Also, when looking at all the simulations (not only those which fulfil the meteoritic constraints), we note that the amount of heavy elements accreted is not strongly dependent on the initial surface density of planetesimals, ranging between 13 M_⊕ and 21 M_⊕ when the mass of Jupiter is attained (bottom-left panel of Fig. 6), and reaching a maximum of 25 M_⊕ when the initial core mass is M_core,1 = 21 M_⊕. This means that a formation location inside the iceline is only consistent with the lower bound of the estimated heavy-element mass inferred from Jupiter’s structure models, and therefore resembles a less likely formation scenario.

For an assumed formation location of a = 10 au, it is found that the total heavy-element mass depends strongly on the prescription for solid accretion. The prescription of SI08 gives a larger Ṁ_Z for larger a (see Eqs. (11) and (13)). In this case the amount of heavy elements is in the range of 20–80 M_⊕, the upper value being beyond the estimated heavy-element mass of Jupiter for standard adiabatic models.

It should be noted that the reported accreted heavy-element masses correspond to the same assumed range of initial planetesimal surface density (5 ≲ Σ₁ ≲ 25 g cm⁻²) for all the simulations. In classical disc models, where the solid surface density decreases outwards, the assumed Σ₁ might seem unfeasible for certain locations, or might give the impression that completely different scenarios are being compared. As a result, we plot, in Fig. 7, the heavy-element mass attained at Jupiter mass for different formation locations, assuming the same Σ₁ in terms of thecorresponding MMSN values at the given location. For this figure we have added a case of Jupiter forming at a = 3 au, as assumed by Desch et al. (2018).

When Σ₁ is equivalent to an initial planetesimal surface density of only twice the MMSN (Fig. 7, left panel), at a = 10 au the growth is too slow and Jupiter’s mass cannot be reached within 10 Myr. At a = 1 au, the solid accretion rate is also low, so at Jupiter mass the heavy-element mass reaches a maximum of ~20 M_⊕. The cases of a = 3 and a = 5 au reach the lower bound of the heavy-element content as inferred by structure models.

Also when Σ₁ is equivalentto an initial surface density of planetesimals of five times the MMSN (Fig. 7, central panel), we confirm that a formation location of a = 1 au is very unlikely: only massive cores can reproduce the required enrichment of Jupiter, and for these cases the meteoritic time constraints of K17 are not met. Interestingly, we find that the case of a = 10 au is also very unlikely: only two initial core masses lead to the desired enrichment, and these cases also do not satisfy the meteoritic constraints. For a = 3 au and a = 5 au, a wide range of M_core,1 can reproduce Jupiter’s enrichment (in agreement with Shibata & Ikoma 2019), and for both set of simulations, four different initial core masses match the meteoritic record.

An intriguing aspect of all panels is the large differences between a = 1 au and a = 3 au. For a = 1 au, the longest phase of solid accretion happens after reaching the planetesimal isolation mass (see left panel of Fig. B.1), while for a ≥ 3 au, the planetesimal isolation mass is never reached, similarly to a = 5 au (also like in Fortier et al. 2007). This makes all cases beyond a = 3 au have a different planetesimal accretion rate, on average larger than at a = 1 au, and decreasing with larger orbital distance. The underlying reason for this behaviour is further discussed in Appendix B.

The right panel of Fig. 7 shows results assuming Σ₁ = 10 Σ_MMSN. In this extreme case in terms of available planetesimal mass, we see that a formation at a = 1 au would become possible. The same is true when a = 10 au, although without satisfying the meteoritic constraints. For a = 3 au all the assumed M_core,1 reproduce Jupiter’s enrichment, but no matching with K17 is found (Jupiter formation is too fast). For a = 5 au, both metallic and K17 constraints are still met.

Therefore, if we assume that Jupiter consists of 20–40 M_⊕ of heavy elements as suggested by standard structure models, it implies that after attaining the initial core (~5–15 M_⊕), Jupiter spent most of its time between ~1 and 10 au (being 1 and 10 au unlikely locations). Formation locations ranging between a = 3 and 5 au provide a very good match with Jupiter’s current bulk metallicity and with the meteoritic constraints of K17; for a wide range of the assumed Σ₁, ranging from two to ten times the MMSN. For Σ₁ corresponding to one MMSN, the metallicity of the forming Jupiter is too low and the meteoritic constraints cannot be matched. This suggests that in our scenario, Jupiter’s formation requires a “minimum planetesimal disc” corresponding to at least two times the MMSN.

It is however, hard to provide robust conclusions regarding Jupiter’s formation location and growth history since the type of the accreted solids as well as their distribution within the disc (and its time dependence) are not well constrained.

Finally, heavy elements can also be accreted by Jupiter after the gas disc has disappeared, and therefore our inferred metallicities should be taken as lower bound. Therefore it is important to acknowledge that in order to put tighter constraints on Jupiter’s origin additional information is required. Further hints could come in the form of cosmos-chemical and isotopic data, comparative planetary science (i.e. comparing Jupiter and Saturn), formation models that account for the formation of the entire planetary system, and information from giant exoplanet data.

Fig. 6 Grid of simulations for a = 1 au (left) and a = 10 au. rp = 100 km. Top: mean time to reach the scattering mass as a function of initial surface density of planetesimals and initial core mass. The light green colour indicates the cases compatible with the meteoritic record. For grey circles the growth is too slow, and for black ones, too fast. Bottom: total amount of heavy elements acquired by the planet when attaining Jupiter mass. Solid circles highlight the cases that match the meteoritic constraints. The accretion rate of solids is taken as that of SI08 in the case of orange circles, and of OG in the case of purple circles.

Fig. 7 Comparison of Jupiter’s bulk metallic content between different formation locations, assuming the initial planetesimalsurface density (Σ1) corresponds to 2, 5 and 10 MMSN at their given location. Σ1 for each run is shown in Table D.2. The initial core mass for all cases is the same as the one assumed in the nominal setup: 4 ≤ Mcore,1 ≤ 17 M⊕. The lighter colours remark the cases where the meteoritic constraints of K17 are met. The accretion of solids in the detached phase is given by SI08.

3.3 Connection to Jupiter’s internal structure

While linking Jupiter’s origin and its current-state interior is challenging, it is clearly desirable to compare the inferred heavy-element mass and its distribution from structure and formation models. Recent structure models of Jupiter that match the gravitational moments measured by the Juno spacecraft, suggest that Jupiter consists of 23.6–46.2 M_⊕ of heavy elements, with a minimum core of 6.21 M_⊕ (Wahl et al. 2017). Interestingly, it was also found that the only structure models that fit the new gravitational data are ones with inhomogenous distribution of heavy elements. The core in these adiabatic models can be either compact (central region of Z_core = 1) or diluted (Z_core < 1).

Alternatively, it is possible that the planetary interior is non-adiabatic and that the distribution of heavy elements changes gradually, with the heavy elements decreasing towards the atmosphere (e.g. Leconte & Chabrier 2012; Vazan et al. 2018; Debras & Chabrier 2019). Indeed, such gradual Z-profile could be expected from formation models by core accretion when the vaporisation of the incoming planetesimals is included (Lozovsky et al. 2017; Bodenheimer et al. 2018). Note however that Lozovsky et al. (2017) computed the heavy element mass deposition to pre-computed thermal/growth models, where the mass of all solids is added to the core. The work of Bodenheimer et al. (2018) is the first to compute self-consistently the mass deposition with the change of composition in the structure equations. Their models, however, are only developed for the formation of planets inside the iceline accreting pure silicate-planetesimals. Moreover, they follow the planetary growth up to only ~8 M_⊕. Thus, the role of compositional gradients in the giant planet formation and the expected distribution of heavy elements in young giant planets are still being investigated.

In our formation models, we assume that the protoplanet’s envelope is convective and it mixes the accreted heavy-elements quickly and efficiently. In all core accretion simulations, the first phase of growth is dominated by the accretion of heavy elements, resulting in a nearly pure-Z deep interior. Phase 2, of both heavy-element and gas accretion, creates an intermediate “layer” of lower metallicity, but still higher than the one expected during the third phase of runaway gas accretion, when most of Jupiter’s mass is accreted, and the metallicity decreases. We illustrate this behaviour in Fig. 8 where we show, for all the grid of simulations corresponding to the nominal case, the mean planetary metallicities resulting from the three mentioned phases of accretion. The core, presumably formed by pebbles, corresponds to the central region where the metallicity ranges between 0.91< Z₁ ≤ 1. The middle layer, which is formed during Phase 2, has metallicity values of 0.24≤ Z₂ ≤ 0.47. The outer layer, corresponding to Phase 3, when gas is accreted in a runaway fashion and planetesimal accretion continues, has a metallicity that ranges between 0.004 ≤ Z₃ ≤ 0.01 with the OG prescription and 0.009 ≤ Z₃ ≤ 0.06 with SI08 (see Fig. 8). The cases with the smallest planetesimal accretion during runaway (lowest Σ₁, see Fig. 8) can lead to sub-solar atmospheric metallicities.

A few of the formation paths we consider are consistent with the estimated metallicity of Jupiter’s atmosphere. The atmospheric metallicity inferred from the Galileo probe is Z_Gal = 0.0167 ± 0.006 (e.g. Debras & Chabrier 2019). The nominal models that yield an outer metallicity in that range corresponds to all those with 6 ≤ Σ₁ ≤ 14 g cm⁻² with the SI08 prescriptions. It is important to note that these cases also lead to a total heavy-element content compatible with what is inferred from structure models that fit Juno data (Wahl et al. 2017, see our Fig. 3). In particular, the models that fulfil the three constraints (Galileo measurement, structure models that fit Juno data, and K17) are the ones with a total heavy-element content of 24–33 M_⊕, of which 2–7 M_⊕ are accretedduring the detached phase (phase 3).

Fig. 8 Jupiter’s metallicity profile resulting from the 3 phases of accretion. All models attaining Jupiter mass within 10 Myr of disc evolution and corresponding to the nominal case with the SI08 prescription are shown. The values of atmosphericmetallicity compatible with the Galileo probe measurements are highlighted in green.

3.4 Mass-metallicity relation: link to giant exoplanets

The metallic content of giant exoplanets was first studied by Guillot et al. (2006). Since then, it has been suggested that there is a correlation between the amount of heavy elements and the total mass of giant exoplanets (Miller & Fortney 2011; Thorngren et al. 2016). The amount of heavy elements in a planet, however, cannot be directly measured and is inferred from models. The idea behind it is to take giant exoplanets with measured masses and radii, then together with the estimate of the stellar age one can run evolution models and search for the heavy-element mass (and therefore, bulk composition) that is required to fit the data. In particular, by computing the cooling of a sample of 47 giant exoplanets, Thorngren et al. (2016, hereafter T16) found that the amount of heavy elements in giant exoplanets follows the power law: $M_{\text{Z}}~M_{\text{P}}^{\Gamma }$, with Γ ≈ 0.6. Some studies have tried to reproduce such relation from formation models (Mordasini et al. 2014, 2016). In particular, Mordasini et al. (2014) find from population synthesis that Γ ≈ 0.12−0.32 (inferred from Γ = 1 + α from the range of α shown in their Table 4). Recently, Hasegawa et al. (2018) found from very simple analytical arguments that the T16 relation could be reproduced with planetesimal-gapped accretion. However, these simple estimates assumed a functional form of the SI08 prescription that is not the general one, but for a uniform planet density, as clarified in Hasegawa et al. (2018), Sect. 5.1 (see Sect. 2.1.2). As shown in Appendix C that the general expression of the gap-planetesimal accretion from SI08 leads to Γ = 0.4 when combined with gas accretion given by hydrodynamical considerations.

In Fig. 9, we show the M_Z−M_P relation we obtain from our simulations when adopting the SI08 prescription (top panel) and the OG prescription (bottom panel) during the detached phase. For the top panel, we note that at the beginning of the detached phase, the planetesimal accretion occurs without opening a planetesimal gap, and afterwards, for planet masses larger than ~150–200 M_⊕, a gap is formed. This happens due to the decrease of the gas accretion rate originated by the opening of a gap in the gas disc (Lissauer et al. 2009). The reduction of Ṁ_gas abates the expansion of the feeding zone, which diminishes the amount of planetesimals that can access it. The planet mass at which the transition between non-gapped and gapped planetesimal accretion occurs depends on the gas accretion rate. Lower gas accretion rates in the detached phase make the transition occur at lower planet masses.

To visualise the range of possible Γ for the M_Z−M_P, the figure shows different power-law relations for the different regimes. We note that the case of Γ =0.6 is extreme and can occur only for the gapped-planetesimal accretion in very metal-rich discs (large Σ₁), and for planets below Jupiter’s mass.

Overall,we find a good agreement between our numerical and analytical results (the latter are shown in Appendix C). For the gapped-planetesimal accretion, our numerical results yield 0.1 ≲ Γ ≲ 0.6, while our analytical estimate provides Γ = 0.4. For the non-gapped planetesimal accretion, our numerical results follow a clear power-law with 0.06≤ Γ ≤ 0.23, while our analytical estime gives Γ = 1∕6 (see Appendix C). These slopes are shown in Fig. 9 for clarity.

Our results should not be interpreted as a disagreement with what is inferred from observations. Following one particular growth case (one curve in Fig. 9) is not, strictly speaking, the correct way to compare with observations. The exoplanets in the T16 sample orbit stars with different properties, like different metallicities, masses and ages. If, for instance, there were some bias towards more massive planets detected around more metallic stars (which is apparently the case, according to Fig. 9 of T16), that would be a way to make the observed relation steeper according to our Fig. 9 (if a larger stellar metallicity implies a larger planetesimal surface density). In this regard, the SI08 prescription seems to be a promising recipe to reproduce the slope of the T16 mass-metallicity relation via population synthesis. This will be the topic of a follow-up work.

Another important aspect presented by T16 is the inferred absolute mass of heavy elements in the planetary sample. They find very large enrichments, with some exoplanets of Saturn mass reaching M_Z≈ 70 M_⊕, and of Jupiter mass reaching M_Z ≈ 100 M_⊕. From our formation models we find that this is very hard to reproduce.

Only when assuming extremely high surface density of planetesimals we managed to reach M_Z ~ 80 M_⊕ (Fig. 6, right panel). This case corresponds to assuming an initial planetesimal-to-gas ratio of ~ 0.6 at a = 10 au, much higher than the usual 0.01 value from classical models. This might suggest that giant planets with large metallic content could form at disc locations with extremely large planetesimal-to-gas ratios. Alternatively, large migration excursions could enhance the total amount of accreted solids (Shibata et al. 2019).

Fig. 9 Evolution of total amount of heavy elements as a function of planet mass during the detached phase, following SI08 (top) and OG (bottom). For clarity, only cases that fulfil the meteoritic time constraints are shown. The final mass for each track is the mass of Jupiter. The grey colour-bar indicates the initial planetesimal surface density (at time 1 Myr). Forthe top panel, the thick(thin) lines show cases with(without) planetesimal gap. For both panels, the green shaded area delimits the range of heavy elements in Jupiter inferred by Juno, and the colour-dashed curves show different $M_{\text{Z}}=\mathcal{C}{M}_{\text{P}}^{\Gamma }$ relations, with ($\mathcal{C},\Gamma$) = (15.5, 0.06) for the brown, ($\mathcal{C},\Gamma$) = (12.33, 1/6) for the red, ($\mathcal{C},\Gamma$) = (10.6, 0.23) for the purple, ($\mathcal{C},\Gamma$) = (12.5, 0.1) for the green, ($\mathcal{C},\Gamma$) = (3.52, 0.4) for the orange, and ($\mathcal{C},\Gamma$) = (1.48, 0.6) for the skyblue curve.

4 Discussion

4.1 Formation scenario: hybrid pebble-planetesimal model vs. pure pebble and pure planetesimal model

Our model assumes the scenario proposed by K17 (but our results are weakly dependent on this assumption, Sect. 3.1.2), where it is claimed that carbonaceous and non-carbonaceous chondrites were separated not only in orbital distance but also during a lapse time of at least 2 Myr. To fulfil this requirement, heat provided by planetesimal accretion is required, as shown in A18.

It is important to note that even if the K17 interpretation were taken with skepticism, when trying to understand the growth of Jupiter with a pure pebble accretion model, one runs into the problem of explaining the metal content of Jupiter inferred by Juno. This should be ~24–46 M_⊕. In a pure pebble accretion model, when the protoplanet reaches the pebble isolation mass, the accretion of heavy elements is halted, because the pebbles are trapped in the pressure bump generated in the vicinity of the protoplanet. This means that the amount of heavy elements accreted by Jupiter should correspond to the pebble isolation mass. The pebble isolation mass depends on the disc aspect ratio, viscosity, and radial pressure gradient (Ataiee et al. 2018; Bitsch et al. 2018). A pebble isolation mass larger than ~25 M_⊕ (approximate minimum amount of M_Z inferred by Juno) requires discs with very large aspect ratios (h∕H > 0.05) and/or very large viscosities (α > 10⁻³) (Ataiee et al. 2018; Bitsch et al. 2018). Thus, a pure pebble accretion model for Jupiter seems to run into difficulties when trying to explain the actual metal content of the planet.

It should be noted that pure planetesimal accretion models have tackled the problem of forming Jupiter since decades (Pollack et al. 1996; Alibert et al. 2005c; Hubickyj et al. 2005; Lissauer et al. 2009; Fortier et al. 2009; Guilera et al. 2011). Despite of the long core formation timescales associated with this scenario, especially when the gravitational interaction between the embryo and planetesimals is considered (oligarchic growth, Inaba et al. 2001; Fortier et al. 2007); these models, under certain assumptions, can lead to the formation of Jupiter. Typically, it is required to assume small planetesimals in order to reach crossover mass before disc dissipation (Fortier et al. 2013; Guilera et al. 2011). Shortening the formation timescale can be achieved when considering additional physical processes such as opacity reduction due to grain growth and settling (e.g. Movshovitz et al. 2010), or heavy-element pollution of the planetary envelope (Hori & Ikoma 2011; Venturini et al. 2016). The amount of heavy elements in Jupiter predicted by these different planetesimal-based models has a large range. For instance, models that invoke low dust opacities can form a Jupiter with a core of ~5 M_⊕ (Hubickyj et al. 2005), models that assume large dust opacities and large surface densities of solids can yield a core of up to ~40 M_⊕ (Guilera et al. 2010). Still, it s important to remark that no pure-planetesimal accretion model can form Jupiter with 100 km planetesimals when the interaction between embryo-planetesimals is correctly accounted for (Fortier et al. 2013; Guilera et al. 2014). This is overcome in the hybrid pebble-planetesimal accretion model.

Finally, it should also be noted that several studies suggested that the standard Phase 2 of giant planet formation, also in the planetesimal accretion scenario, can be extremely short or even non-existing (Zhou & Lin 2007; Shiraishi & Ida 2008). In that case, when the formation of a gap in the planetesimal disc is considered, the core mass is expected to be relatively small with a significant planetesimal accretion when the planet reaches a mass of ~ 100 M_⊕ (Shiraishi & Ida 2008). While this is an interesting option, we suggest that further calculations which account for the planetary internal structure and limit the planetesimal supply must be considered. In addition, it is yet to be shown that such a formation scenario with no Phase 2 is consistent with the cosmo-chemical constraints discussed in this study.

4.2 Water abundance and formation location

The Juno spacecraft is also expected to measure the water abundance in Jupiter using the microwave radiometer (MWR) which probes Jupiter’s deep atmosphere at radio wavelengths (Bolton et al. 2017). Data from MWR is still being collected, and constraints on the water abundance are expected in the near future. It should be noted, however, that Juno’s measurements are limited to the uppermost 100–200 bars, and might not represent the bulk (e.g. Helled & Lunine 2014).

In order to infer the global water abundance in Jupiter, one must assume that the heavy elements are homogeneous within the planet, an assumption that has recently been challenged by Jupiter’s structure models fitting Jupiter’s gravitational field (Wahl et al. 2017; Vazan et al. 2018; Debras & Chabrier 2019). In this case the MWR estimate is likely to provide a lower bound for the water abundance.

In Sect. 3.2 we analysed the possibility of Jupiter spending most of its formation time inside the water iceline. We found that in such case the heavy-element mass is in the lower bound of the range inferred from structure models. In terms of pebble isolation mass, at 1 au only a value of ~ 6 M_⊕ is expected (Lambrechts et al. 2014). This means that to fulfil the K17 scenario, a core of 16 M_⊕ should have migrated inwards, from beyond the iceline. In that case, half of it should be made of ices (e.g. Thiabaud et al. 2015), which would give Jupiter a maximum amount of ~8 M_⊕ of water, and probably very concentrated towards the centre of the planet (the accretion of dry solids in our simulations would be by planetesimals, in Phase 2 and Phase 3 of growth). Note, however, that this limit would increase if Jupiter accreted, post-disc dispersal, icy planetesimals/embryos.

From the simulations beyond the iceline (Figs. 2 and 9), we find that Jupiter could have a range of ~ 20−75 M_⊕ of heavy elements. The upper value is linked to the formation test case at a = 10 au (Sect. 3.2). In the nominal formation location of 5 au, the upper limit on Jupiter’s heavy elements is found to be M_Z= 43 M_⊕ (Fig. 3). According to condensation sequence studies, half of these heavy-elements could be volatiles (Thiabaud et al. 2015), which, if we neglect the presence of other volatiles than water, sets a maximum amount of water in the range of ~ 12–22 M_⊕. If beyond the iceline water resembles at least 50% in mass of the condensed material (Thiabaud et al. 2015), this leads to a minimum amount of water of ~6 M_⊕ for Jupiter.

The above considerations imply that if the maximum amount of water in Jupiter is below ~6 M_⊕, this would suggest that Jupiter spent a few million years within the iceline. In case the water content were confirmed to be above ~10 M_⊕, the conclusions are more blurry, because Jupiter could have spent some time within the iceline, but still have been polluted by ices after the disc dispersal. Still, if Jupiter had spent phase 2 and 3 inside the water iceline, it would have a larger rock-to-ice ratio compared to a full formation beyond the iceline. In the latter, the rock-to-ice ratio shouldbe close to 1. In the former (phase 2 and 3 within the iceline), Jupiter would have started with 8 M_⊕ of rocks and the same amount of ices, then accreted maximum 8 M_⊕ of rocks (to get a total of 24 M_⊕ of heavy elements, as expected by Juno). Thus, the rock-to-ice ratio would be closer to ~2.

The above discussion assumes the water iceline to be the location in the disc where planetesimals form with temperatures of ~170 K, and that the planetesimals do not change their composition during the disc evolution. In reality, the water snow line is expected to move as the disc evolves (e.g. Oka et al. 2011), but this would modify mainly the composition of pebbles (Morbidelli et al. 2016). In addition, the water abundance in the accreted solids is unknown. Therefore, linking Jupiter’s water abundance with its formation process is non-trivial: it depends on the formation location of the planet, its migration history, and the composition of the accreted planetesimals/pebbles. As a result, the inferred numbers should be taken as a guideline as they are derived under simplified assumptions. Finally, if certain evaporative mechanism operates efficiently, like heating from Al-26 (Lichtenberg et al. 2019), the above inferences would be only an upper bound on water content. Jupiter could be much drier than expected even if it spent phase 2 and 3 beyond the water iceline. In summary, constraining the planetary origin from water abundance estimates is a difficult task, and it suffers from various uncertainties and degeneracies in theoretical calculations.

Further constraints might come from linking cosmochemical data with astrophysical models. Indeed, Desch et al. (2018) show that in order to solve the CAI storage problem, the second phase of Jupiter’s growth had to be spent between a = 2.9 and a = 3.2 au. If Jupiter had been farther out, the amount of CAIs to form chondritic chondrites beyond Jupiter would have been too scarce. We found formation models of Jupiter that match both the estimated metallicity of the planet and the K17 constraints for 1 < a < 10 au, being a = 3 au a perfectly possible location for Jupiter to spend phase 2 (see Sect. 3.2). Many existing models explain different aspects of Jupiter’s observed properties, but clearly more work on linking meteoritic data, disc thermal evolution, simultaneous formation of Jupiter and Saturn with the correct metallic content, and special migration regimes (e.g, Masset & Snellgrove 2001) is required in order to develop a coherent picture of the history of our Solar Sytem.

5 Conclusions

We simulate Jupiter’s growth and followed the heavy-element accretion for various model assumptions. Many of our models can reproduce the amount of heavy elements in Jupiter as inferred from structure models (Wahl et al. 2017).

Similarly to Alibert et al. (2018), we show that the initial phase of core formation (phase 1) is likely to be dominated by pebble accretion, followed by a second stage (phase 2) of planetesimal accretion. By modelling the second stage coupled self-consistently with gas accretion up to Jupiter mass (phases 2 and 3 in our scenario), we find that the meteoritic time constraints derived by K17 can be successfully matched with the accretion of 100 km-size planetesimals. This is an important finding, because it is the first time that large planetesimals are shown to play a crucial role in Jupiter’s growth and enrichment when the excitation of planetesimals by the embryo is consistently accounted for (oligarchic growth regime). It is also important to note that the total heavy-element mass in Jupiter inferred from our simulations is independent on the assumption that Jupiter is responsible for reconnecting the carbonaceous and non-carbonaceous reservoirs (K17, A18).

We suggest that a formation location inside the water iceline after reaching pebble isolation mass is less likely for Jupiter, since the amount of heavy elements that can be accreted in the inner regions of the disc is much lower than at 3–5 au, making the total heavy-element mass accreted by Jupiter in the lower-bound of the heavy-element mass estimates from structure models. Our calculations favour Jupiter’s formation at 1 < a < 10 au after pebble isolation mass is reached.

We find that our hybrid pebble-planetesimal model can produce a Jupiter that matches the estimates for the planet’s bulk metallicity and the meteoritic time constraints derived by K17. In addition, when examining the average metallicity for a given formation phase, and assuming that the layer corresponding to this formation phase is homogeneously mixed, we can reproduce Jupiter’s atmospheric metallicity as inferred by the Galileo probe. For these models, the heavy element-mass accreted during the detached phase is found to be M_Z ≈ 2−7 M_⊕.

We also find that Jupiter can accrete ~1–15 M_⊕ of heavy elements during its final formation phase. The exact value depends on the initial surface density of planetesimals and their dynamical properties, which are highly unconstrained. We find that if the initial planetesimal surface density corresponds to approximately two times the MMSN, the accreted heavy-element mass during Phase 3 would be ~1 M_⊕, and of ~6–15 M_⊕ for ten times the MMSN. This is an important finding since it has direct implications of our understanding of the enrichment of giant planets in general. If a non-negligible amount of heavy elements can be accreted in this last stage of formation, it would lead to a higher metallicity in the outer planetary envelope. For the case of Jupiter, we find that such “late enrichment” can yield an envelope metallicity of ~0.5–3 times solar. This results suggest that the diversity in planetary composition is a natural outcome of planet formation models, with the exact composition depending on the exact formation location and growth history.

We have, additionally, investigated the link to the metallicity of giant exoplanets. We find that the very high bulk metallicities inferred for many exoplanets in the Thorngren et al. (2016) sample are hard to reproduce. It is found that an extremely high initial planetesimal surface density, corresponding to an initial planetesimal-to-gas ratio of ~0.6, could yield M_Z ~ 80 M_⊕ at Jupiter mass and a = 10 au. This could suggest a migration history for the metal-rich warm giant exoplanets, as recently proposed by Shibata et al. (2019).

Finally, we examined the slope of the mass-metallicity relation of giant exoplanets, and found analytically that $M_{\text{Z}}~M_{\text{P}}^{\Gamma }$, with Γ ≈ 1∕6 when no gap in the planetesimal disc is created, and with Γ ≈ 0.4 with planetesimal gap. These values are in general agreement with our numerical results, but are lower than the 0.6 value found by T16. The solid accretion rate provided by Shiraishi & Ida (2008) could lead to a better match with T16 when implemented in population synthesis models, and we hope to address this point in future work.

Acknowledgements

We thank Sho Shibata and Masahiro Ikoma for enlightening discussions, and Octavio Guilera for relevant suggestions. We also thank the anonymous referee for the insightful comments, which helped us to improve the manuscript. This work has been carried out, in part, within the framework of the National Centre for Competence in Research PlanetS, supported by the Swiss National Foundation. R.H. acknowledges support from SNSF grant 200021_169054.

Appendix A Details on planetesimal accretion: the importance of the capture radius and the initial large core mass

A long standing problem in planet formation has been the slow planetary growth when the solid accretion is given solely by large planetesimals. Indeed, simulations show that the formation timescale of gas giants typically exceeds discs’ lifetimes (Pollack et al. 1996; Fortier et al. 2007). Two broad solutions for this “timescale problem” were proposed:

• (i)

  Planets grow mainly by accretion of small planetesimals (100 m–1 km) in discs with relatively large metallicities (Fortier et al. 2009, 2013; Guilera et al. 2011; Alibert et al. 2013).

• (ii)

  Planets grow predominantly by cm-size pebbles (Lambrechts et al. 2014).

One problem with the first scenario is that planetesimals seem to be formed large (~100 km) by streaming instability (e.g, Simon et al. 2017). In addition, when collisions are included, planetesimals in the size range of ~km do not live long (Kobayashi et al. 2010). On the other hand, the pure pebble accretion scenario runs into other type of difficulties. First, pebbles drift fast towards the central star, and protoplanetary discs have limited sizes of ~10–200 au (Andrews et al. 2010), so the prevalence of pebbles in the disc during the ~2–10 Myr of discs’s lifetime is not guaranteed. Second, as the disc evolves and cools, icy pebbles can reach the inner regions of the system, polluting it with volatiles (Morbidelli et al. 2016). Recent interpretation of the Kepler data suggests that super-Earths are dry (Owen & Wu 2017), which is at odds with a pure pebble accretion scenario. In addition, pebble accretion is either too efficient (discs with low viscosity) or too inefficient (discs with high viscosity), which leads to the formation of either too many gas giants or too manyembryos (Lin et al. 2018; Venturini & Helled 2017). Indeed, intermediate-mass planets (Neptunes), which are very abundant (Batalha et al. 2013; Suzuki et al. 2018), are hard to form only with pebbles (Venturini & Helled 2017).

Fig. A.1 Growth of a planet corresponding to the nominal setup, with Σ1 = 15 g cm−2. The green curves assume an initial core mass of Mcore,1= 10 M⊕ ; and the orange curves an initial core mass of Mcore,1= 0.01 M⊕. Solid lines indicate the growth of the core, dashed lines the growth of the envelope. The envelope mass of the latter case is always smaller than 0.001 M⊕.

A hybrid pebble-planetesimal accretion scenario for giant planet formation does not only provide a solution for the isotopic anomalies of carbonaceous and non-carbonaceous meteorites (Alibert et al. 2018), but it also has the potential to overcome the aforementioned problems. An initial dominant phase of pebble accretion is plausible, given that discs are expected to be very abundant in pebbles at the beginning of their lifetimes (e.g. Lambrechts & Johansen 2014).

Fig. A.2 Upper panel: growth of a planet for Mcore,1 = 10 M⊕ and Σ1 = 15 g cm−2, with (green) and without (purple) the effect of the capture radius. The solid lines show the growth of the core and the dashed ones the growth of the envelope. Bottom panel: evolution of capture radius in units of core radius for the growth cases shown in the top panel. In the purple case the capture radius is always fixed to the core radius. The green curve drops to a radius of 1.8 RJup when the detached phase starts.

In the light of this discussion, it is important to clarify why in the new pebble-planetesimal accretion scenario, the accretion of large planetesimals is consistent with the formation of a gas giant. The key for 100 km-planetesimal accretion to be successful is twofold:

• 1.

  The initial core mass cannot be the usual lunar mass (value that came from the original pure N-body simulations of oligarchic growth, Kokubo & Ida 1998). It has to be at least 4 M_⊕. This core can be formed easily by pebble accretion.

• 2.

  The atmospheric enhancement factor for the accretion rate (the capture radius) has to be considered (Inaba & Ikoma 2003; Fortier et al. 2009).

We illustrate the importance of these two effects in Figs. A.1 and A.2. The former shows the planetary growth with the nominal set up, and initial planetesimal surface density of Σ₁ = 15 g cm⁻². One growth-case assumes an initial core mass of M_core,1 = 10 M_⊕ and the other, of M_core,1 = 0.01 M_⊕. When the core starts its growth with a lunar mass, at t = 10 Myr it reaches a core mass of only M_core = 0.24 M_⊕, whereas in the case where M_core,1= 10 M_⊕, crossover mass (25 M_⊕ in this case) is attained at t = 4 Myr. This means that for the case that starts with a lunar embryo, the average accretion rate of solids is 0.23 M_⊕ per 10 Myr (or 2.3 × 10⁻⁸ M_⊕ yr⁻¹); while in the nominal case the mean accretion rate of solids is 15 M_⊕ per 4 Myr (or 3.75 × 10⁻⁶ M_⊕ yr⁻¹). The planetesimal accretion rate during the OG goes with $R_{\text{H}}^2~M_{\text{P}}^{2/3}$, which gives a ratio between the two starting embryos of ${(10/0.01)}^{2/3}$ = 100. This is practically the same difference in ratio that we obtain between the two mean planetesimal accretion rates.

The effect of the capture radius on the core’s growth is demonstrated on Fig. A.2, where, as before, Σ₁ = 15 g cm⁻² and M_core,1 = 10 M_⊕. In one simulation the capture radius is computed and included in the cross section of accretion, and in the other case, the capture radius is set to the core radius for the whole simulation (i.e, the enhancement factor of accretion is neglected). The effect on the core growth is shown on the top panel of Fig. A.2. The bottom panel shows the value of the capture radius in units of core radius for both setups. We can see that without accounting for the capture radius, the core grows too slowly to attain crossover mass before the dissipation of the disc.

Appendix B Difference in solid accretion between growing Jupiter inside or beyond the iceline

We mentioned in Sect. 3.2 that the formation scenario at a = 1 au is unfavourable for Jupiter’s formation due to the relatively low bulk metallicity of Jupiter. Although this formation scenario is not very likely, it is interesting to understand the different behaviours between Figs. 2 and 6. We therefore compare the growth of the heavy elements and total planet mass for the cases of 1 and 5 au in Fig. B.1. For clarity, only cases where M_core,1 = 16 M_⊕ (a = 1 au) and M_core,1 = 13 M_⊕ (a = 5 au) are shown. The different influence of Σ₁ at both locations is very clear. The green-shaded area shows the cases that match the meteoritic time constraint. For an equal M_core,1, all the values of Σ₁ produce growth paths that match the meteoritic time constraints for a = 1 au, whereas for a = 5 au, this depends on the value of Σ₁.

The growth of the heavy-element mass is very different in both cases. In the former, the core grows quickly at the beginning, consuming all the planetesimals available in its feeding zone in the first ~300 thousand years. The planetesimal isolation mass is a strong function of the semi-major axis, with a milder dependance on the initial planetesimal surface density ($M_{\text{iso}}~a^3{\Sigma }_1^{3/2}$, Pollack et al. 1996). This implies that at a = 5 au, the planetesimalisolation mass is ~100 larger than at a = 1 au (Σ₁ being constant). Consequently, the case at a = 5 au never reaches the planetesimal isolation mass: there are always enough planetesimals to be accreted in the feeding zone. Indeed, we checked that at a = 5 au, the amount of planetesimals accreted at each time step until crossover mass never exceeds 1% of the total amount of planetesimals available in the feeding zone at that time. The duration of phase 2¹ is less dependent on Σ₁ than phase 1. The main factor influencing the duration of phase 2 is the composition of the envelope, which is the same for all the cases plotted here. It is exactly the fact that the planetesimal isolation mass is so much smaller at a = 1 au than at a = 5 au what makes it difficult for a giant planet to accrete substantial amounts of solids in close-in orbits: there is simply not enough material.

Fig. B.1 Growth of planet at a = 1 au (left) and a = 5 au (right).For clarity, only cases with an initial core of 16 M⊕ and 13 M⊕ are shown (left and right panels, respectively). Top panels: evolution of the mass of heavy elements. Bottom panels: evolution of the total planet mass. Colour-bar: initial surface density of solids. Green shaded area: region where the meteoritic time constraints are met (up to MP ~ 50 M⊕).

Appendix C An Analytical estimate of the M_P –M_Z relation from formation

In this section we derive analytically the M_P –M_Z relation assuming the accretion rate of solids is given by Shiraishi & Ida (2008, hereafter SI08) and that the gas accretion is dictated by the supply of the disc (Tanigawa & Watanabe 2002). We follow the same approach as Hasegawa et al. (2018, hereafter H18), but using the general form of the planetesimal accretion from SI08 (Eqs. (10) and (12)), and assuming that the protoplanet’s radius has shrunk to a fixed value, a result that is found by several studies during the detached phase (Lissauer et al. 2009; Mordasini 2013). It should be noted that H18 assume that the density of the protoplanet remains uniform during this final stage of accretion, i.e, $R_{\text{P}}~M_{\text{P}}^{1/3}$, which is not what is found when computing the actual radius during the detached phase.

We start by analysing the gapped-planetesimal case. From Eq. (12) we have the following functional dependance, $${\dot{M}}_{\text{Z},\text{gap}}~{\rho }^{1/2}{R}_{\text{P}}^2{(\frac{v_{\text{H}}}{v_{\text{damp}}})}^{\alpha }.$$(C.1)

Considering that ${\tau }_{\text{g}}~M_{\text{P}}^D$ and using Eq. (13) only with the mass dependence (as in H18), the mass dependence of Eq. (C.1) reads, $${\dot{M}}_{\text{Z},\text{gap}}~M_{\text{P}}^{1/2-\alpha (1/6+D)}.$$(C.2)

The metal content (M_Z) as function of planet mass (M_P) can be obtained from integrating dM_Z/dM_P, which is given by: $$\frac{\text{d}{M}_{\text{Z}}}{\text{d}{M}_{\text{P}}}={\dot{M}}_{\text{Z}}\frac{\text{d}t}{\text{d}{M}_{\text{P}}}~{\dot{M}}_{\text{Z}}\frac{{\tau }_{\text{g}}}{M_{\text{P}}}.$$(C.3)

Using the above expressions, we obtain, for gapped-planetesimal accretion, dM_Z/dM_P = $M_{\text{P}}^{\gamma }$, with $$\gamma =-\frac{1}{2}-\frac{\alpha }{6}+D(1-\alpha ).$$(C.4)

Hence, when integrating dM_Z/dM_P in the planet mass, we get $M_{\text{Z}}=M_{\text{P}}^{\Gamma }$, where Γ = 1 + γ, i.e: $$\Gamma =\frac{1}{2}-\frac{\alpha }{6}+D(1-\alpha ).$$(C.5)

Plugging-in the values for gap-planetesimal accretion (α = 1.4) and for the timescale to accrete gas given by the supply of the disc (see H18) (D = −1∕3), we obtain Γ = 0.4. This exponent of the M_P –M_Z is not far from the one derived by T16 (Γ = 0.61 ± 0.08) and is actually more consistent with our numerical results (top panel of Fig. 9).

Analogously, for the non-gapped planetesimal case, the functional dependence of the accretion rate of solids goes as (Eq. (10)): $${\dot{M}}_{\text{Z},\text{nogap}}~{\rho }^{1/2}{R}_{\text{P}}^2{(\frac{v_{\text{H}}}{v_{\text{scat}}})}^{\alpha },$$(C.6)

where α = 0.8

Applying again Eq. (C.3) and assuming the accretion rate of gas is limited by the disc supply, we reach, by using Eqs. (10) and (11): $$\Gamma =\frac{1}{2}-\frac{\alpha }{3}+D(1-\alpha ).$$(C.7)

Plugging-in the corresponding values of α = 0.8 and D = −1/3, we obtain Γ = 1∕6.

Appendix D Heavy element accretion in each phase of Jupiter’s growth

Table D.1 lists that heavy-element accretion mass during each phase of Jupiter’s growth as inferred from all the simulations presented in this study. The table is constructed as follows. For each run, the maximum and minimum ranges of the assumed Σ₁ are shown, and for each Σ₁, the maximum and minimum values of M_core,1 that yield to a Jupiter mass planet within 10 Myr of disc evolution. Thus, for each run, typically four extreme cases are shown, and all the corresponding intermediate values of (Σ₁, M_core,1) yield an accreted amount of solids that lie in between the shown boundaries. For one case there is only one value of M_core,1 for a given Σ₁. For that case, a lower M_core,1 resulted in a planet that has not reached a Jupiter mass within 10 Myr, and a larger value of M_core,1 yielded a supercritical planet. By comparing the amount of solids accreted in the different phases, one can see that the ones accreted during the attached Phase 3 have always very low values. This is because the detached phase starts typically soon after the crossover time. For a few cases in the run of envelope enrichment, the crossover mass is reached just after the detached phase starts, meaning that there is no attached- Phase 3.

We note that for a given Σ₁, typically the amount of heavy elements accreted during Phase 2 is larger for the smaller M_core,1. This is because the time duration of Phase 2 is considerably longer for the smaller M_core,1. For instance, for the nominal case and Σ₁ = 25 g cm⁻², Phase 2 for M_core,1 = 4 M_⊕ lasts 5.3 Myr and for M_core,1 = 17 M_⊕, 0.6 Myr. This is the reason that the averaged heavy-element accretion rate is always higher for the larger core mass, as indicated in Fig. 2.

Table D.2 provides the same information as Table D.1 but for the runs shown in Fig. 7; where the initial surface density of planetesimals is fixed in terms of the MMSN values (Hayashi 1981): ${\Sigma }_{1,\text{MMSN}}=30.6f_{\text{ice}}{(a/\text{au})}^{-1.5}$, where f_ice = 0.5 inside the water iceline (at a = 1 au in our simulations) and 1 for larger semi-major axes.

References

All Tables

All Figures

Fig. 1 Example of planetary growth, where the initial surface density of planetesimals is Σ1 = 7 g cm−2 and the initial core mass is Mcore,1 = 14 M⊕. MZ is the mass of heavy elements, MHHe the mass of H and He, and MP the total mass of the planet. The solid lines indicate that the protoplanet is in the attached phase, whereas the dotted lines show the protoplanet in the detached phase. Bottom panel: simple zoom to appreciate better the growth of MZ.

In the text

Fig. 2 Grid of simulations for a = 5 au and rp = 100 km. Left: mean time to reach the scattering mass as a function of the initial surface density of planetesimals (Σ1) and the initial core mass (Mcore,1). The light green colour indicates the cases compatible with the meteoritic record. For grey circles the growth is too slow (longer than 4.5 Myr), and for black ones, too fast (shorter than 2.5 Myr). Points in the upper left corner above the red line are supercritical and their structure cannot be computed. Points in the lower left corner below the purple line correspond to cases whose growth is so slow that the mass of Jupiter is not reached within 10 Myr of disc evolution. Right: mean accretion rate of solids during Phase 2 as a function of the initial surface density of planetesimals. The initial core mass is shown in the colour-bar. Circles with black border represent cases where the meteoritic constraints are fulfilled (green band of left panel).

In the text

Fig. 3 Simulations for a = 5 au and rp = 100 km with the same range in (Mcore,1, Σ1) as in Fig. 2. Orange circles correspond to the solid accretion given by SI08, purple circles to Oligarchic Growth (OG). The darker colours represent all runs that reach Jupiter mass within 10 Myr of disc evolution. The lighter colours highlight the cases that match the meteoritic constraints of K17. Left: total amount of heavy elements accreted until reaching Jupiter mass as a function of the initial surface density of planetesimals (Σ1). The green shaded area indicates the values of MZ compatible with estimates from Juno (Wahl et al. 2017). Right: amount of heavy elements accreted during the detached phase with the different solid accretion prescriptions.

In the text

Fig. 4 Grid of simulations for a = 5 au and rp = 10 km. Left: mean time to reach the scattering mass as a function of the initial surface density of planetesimals (Σ1) and initial core mass (Mcore,1). The light green colour indicates the cases compatible with the meteoritic record. For grey circles the growth is too slow, and for black ones, too fast. Right: total amount of heavy elements accreted until reaching Jupiter mass as a function of the initial surface density of planetesimals (Σ1). Transparentcircles: all results. Solid circles: cases compatible with the meteoritic constraints.The green shaded area indicates the values compatible with estimates from Juno (Wahl et al. 2017).

In the text

Fig. 5 Amount of heavy elements accreted during the attached phase (including the initial core formed by pebble accretion) as a function of the initial surface density of planetesimals for different values of dust opacities (circles and squares) and considering or not envelope enrichment (triangles and circles). All the cases that reach Jupiter mass within 10 Myr of disc evolution are shown. All cases correspond to planetesimal sizes of rp = 100 km. Red circles indicate the nominal simulations (f = 1), squared yellow points represent the case with dust opacity reduced by a factor of 10 (f = 0.1), and blue triangles indicate the case where envelope enrichment is included.

In the text

Fig. 6 Grid of simulations for a = 1 au (left) and a = 10 au. rp = 100 km. Top: mean time to reach the scattering mass as a function of initial surface density of planetesimals and initial core mass. The light green colour indicates the cases compatible with the meteoritic record. For grey circles the growth is too slow, and for black ones, too fast. Bottom: total amount of heavy elements acquired by the planet when attaining Jupiter mass. Solid circles highlight the cases that match the meteoritic constraints. The accretion rate of solids is taken as that of SI08 in the case of orange circles, and of OG in the case of purple circles.

In the text

Fig. 7 Comparison of Jupiter’s bulk metallic content between different formation locations, assuming the initial planetesimalsurface density (Σ1) corresponds to 2, 5 and 10 MMSN at their given location. Σ1 for each run is shown in Table D.2. The initial core mass for all cases is the same as the one assumed in the nominal setup: 4 ≤ Mcore,1 ≤ 17 M⊕. The lighter colours remark the cases where the meteoritic constraints of K17 are met. The accretion of solids in the detached phase is given by SI08.

In the text

	Fig. 8 Jupiter’s metallicity profile resulting from the 3 phases of accretion. All models attaining Jupiter mass within 10 Myr of disc evolution and corresponding to the nominal case with the SI08 prescription are shown. The values of atmosphericmetallicity compatible with the Galileo probe measurements are highlighted in green.
In the text

Fig. 9 Evolution of total amount of heavy elements as a function of planet mass during the detached phase, following SI08 (top) and OG (bottom). For clarity, only cases that fulfil the meteoritic time constraints are shown. The final mass for each track is the mass of Jupiter. The grey colour-bar indicates the initial planetesimal surface density (at time 1 Myr). Forthe top panel, the thick(thin) lines show cases with(without) planetesimal gap. For both panels, the green shaded area delimits the range of heavy elements in Jupiter inferred by Juno, and the colour-dashed curves show different $M_{\text{Z}}=\mathcal{C}{M}_{\text{P}}^{\Gamma }$ relations, with ($\mathcal{C},\Gamma$) = (15.5, 0.06) for the brown, ($\mathcal{C},\Gamma$) = (12.33, 1/6) for the red, ($\mathcal{C},\Gamma$) = (10.6, 0.23) for the purple, ($\mathcal{C},\Gamma$) = (12.5, 0.1) for the green, ($\mathcal{C},\Gamma$) = (3.52, 0.4) for the orange, and ($\mathcal{C},\Gamma$) = (1.48, 0.6) for the skyblue curve.

In the text

	Fig. A.1 Growth of a planet corresponding to the nominal setup, with Σ1 = 15 g cm−2. The green curves assume an initial core mass of Mcore,1= 10 M⊕ ; and the orange curves an initial core mass of Mcore,1= 0.01 M⊕. Solid lines indicate the growth of the core, dashed lines the growth of the envelope. The envelope mass of the latter case is always smaller than 0.001 M⊕.
In the text

Fig. A.2 Upper panel: growth of a planet for Mcore,1 = 10 M⊕ and Σ1 = 15 g cm−2, with (green) and without (purple) the effect of the capture radius. The solid lines show the growth of the core and the dashed ones the growth of the envelope. Bottom panel: evolution of capture radius in units of core radius for the growth cases shown in the top panel. In the purple case the capture radius is always fixed to the core radius. The green curve drops to a radius of 1.8 RJup when the detached phase starts.

In the text

Fig. B.1 Growth of planet at a = 1 au (left) and a = 5 au (right).For clarity, only cases with an initial core of 16 M⊕ and 13 M⊕ are shown (left and right panels, respectively). Top panels: evolution of the mass of heavy elements. Bottom panels: evolution of the total planet mass. Colour-bar: initial surface density of solids. Green shaded area: region where the meteoritic time constraints are met (up to MP ~ 50 M⊕).

In the text