© The Authors 2023

1 Introduction

The discovery and characterization of exoplanets over the course of the last 30 yr has brought about a seismic shift in our comprehension of planetary formation and evolution. Nonetheless, the first objects to be discovered in large numbers – hot Jupiters – continue to stand out as an enigmatic class of astrophysical bodies. These giants are close to their host stars (with orbital periods of less than 10 days) and, despite their unexpected nature, have attracted extensive scrutiny due to their distinctive features and relative ease of observation. In particular, precise mass and radius determinations are substantially more common within the presently known census of hot Jupiters than for other types of extrasolar planets.

Conventional giant planet structure theory holds that the mass-radius relation for degenerate Jovian planets is approximately flat (i.e., mass-independent; Stevenson 1982), meaning that their size is largely dictated by their composition. Within this framework, Jupiter’s radius in the first approximation corresponds to a roughly solar mixture of hydrogen and helium, meaning that any smaller radius is indicative of a substantially super-solar overall metallicity. The estimate of the total mass of heavy elements, here denoted as M_h, can be sharpened further through detailed modeling (which accounts for the corrections due to age, total mass, etc.). Guillot et al. (2006) were the first to carry out this analysis for a pool of nine well-characterized hot Jupiters. Intriguingly, they found that some objects have M_h on the order of ~ 100 M_⊕, where M_⊕ denotes the mass of Earth, and pointed out an apparent correlation between the planet’s M_h and the metallicity of the central star. This investigation was later extended and confirmed in Guillot (2008), Laughlin et al. (2011), and Moutou et al. (2013).

An important complication that arises within such analyses is that hot Jupiters experience substantial radius inflation, and, as such, their interiors are not in a fully degenerate state. Though a number of physical mechanisms – including tidal damping (Bodenheimer et al. 2003; Levrard et al. 2007), breaking gravity waves (Guillot & Showman 2002), impeded cooling due to enhanced atmospheric opacity (Burrows et al. 2007), double-diffusive convection (Chabrier & Baraffe 2007), turbulent burial of atmospheric entropy (Youdin & Mitchell 2010), and Ohmic dissipation (Batygin & Stevenson 2010) – have been proposed to explain this anomalous heating within the planets’ envelopes, statistical analyses (Thorngren & Fortney 2018) have shown that the strong dependence of the degree of inflation on stellar irradiation predicted by the Ohmic dissipation mechanism is indeed reflected in the data (see also Knierim et al. 2022). Consequently, the application of conventional giant planet evolution models to strongly irradiated planets can yield negative heavy element masses. In turn, this implies that the values of M_h reported in the aforementioned studies are lower bounds.

To circumvent this problem, Thorngren et al. (2016) considered a subset of 47 giant planets that are not strongly irradiated by their central star. Within this subset of objects, anomalous heating of the interior could be reasonably assumed to be negligible, meaning that the computed values of M_h likely represent the actual masses in heavy elements and not their lower bounds (indeed, no negative values of M_h appear in the Thorngren et al. calculations). We note that, strictly speaking, many of these planets fall into the “warm Jupiter” category because they have periods that exceed the nominal 10 day boundary, but for simplicity we still refer to them as hot Jupiters.

The main result of Thorngren et al. is reproduced in Fig. 1, which shows that many hot Jupiters are more enriched in heavy elements than Jupiter or Saturn. Some Jupiter-mass planets exceed 100 Earth masses in heavy elements. On average, the mass in heavy elements, M_h, is correlated to the total planet mass, M_p, as $$M_{\text{h}}=57.9\pm 7.0M_{\oplus }{\left(\frac{M_{\text{p}}}{M_{\text{jup}}}\right)}^{\beta },$$(1)

where M_jup is the mass of Jupiter. The exponent β of the M_h(M_p) correlation is 0.61 ± 0.08. Thorngren et al. also confirmed the correlation between M_h/M_p (the planet metallicity) and the stellar metallicity¹. When this correlation is accounted for, the scatter of the data around the correlation law given by Eq. (1) is significantly reduced.

These results are surprising. According to the core-accretion theory of giant planet formation, giant planets are nucleated by the gradual accumulation of solid material into a ~15 M_⊕ core, which then accretes a massive envelope of gas and small dust with approximately stellar metallicity. This process is expected to result in a range of solid-to-gas ratios for giant planets that is appreciably super-stellar but nonetheless much smaller than that given by Eq. (1). For instance, a Jupiter-mass planet would be expected to have ~18 M_⊕ of heavy elements. This mismatch between expectations and observations suggests that the process of hot-Jupiter formation may be more complex than originally thought.

Thorngren et al. proposed an explanation for the surprising heavy element enrichment observed in hot Jupiters. They conjectured that planets accrete all planetesimals located within their feeding zone, an annulus with a radial width proportional to the planet’s Hill radius, R_H = a(M_p/3M_star)^1/3, where a is the semi-major axis of the accreting planet. The exponent 1/3 is smaller than that of Eq. (1) but still not grossly inconsistent with the data. Venturini & Helled (2020) proposed a similar explanation for the enrichment in heavy elements in Jupiter and, through a more sophisticated planetesimal accretion model, predicted an exponent of 2/5 when the formation of a planetesimal gap is considered, which is a bit closer to the measured value of β than the estimate of Thorngren et al. The combination of accretion of gas with a 1% metallicity and the accretion of planetesimals with the $M_{\text{p}}^{2/5}$ relationship gives the magenta curve in Fig. 1, which explains some of the planets, but clearly not the majority of them. Shibata et al. (2020, 2022) show that planet migration can enhance the efficiency of planetesimal accretion, due to a combination of resonant shepherding and gas drag. This is particularly efficient when the planet is migrating through specific locations of the disk (depending on the parameters). This may potentially make hot Jupiters more metal-rich than Jupiter itself.

A radically different mechanism for the heavy element enrichment of giant planets has been proposed by Schneider & Bitsch (2021), elaborating on an original idea of Guillot & Hueso (2006). In their model, inward drifting dust particles (also known as pebbles) evaporate their volatile elements, each at a specific distance (the sublimation line of the corresponding volatile species). Therefore, the gas in the inner part of the disk gets enriched in volatile element vapor, and by a substantial amount for some disk parameters. Because a planet accretes H, He, and heavier element vapors indiscriminately, it can be enriched in heavy volatile elements through this process.

It is worth noting that this model is not restricted to hot Jupiters; it also applies to any planet accreting a substantial fraction of its gaseous envelope in the inner part of the protoplanetary disk, where volatile elements are in a vapor state. Nevertheless, this model falls a bit short of reproducing correlation (1). The M_h are typically below the value predicted by correlation (1) up to planets with M_p = 2 M_Jup and can reach the observed mean values only for more massive planets formed in the most metallic disks. Planets with M_h > 100 M_⊕ are not expected in the Schneider and Bitsch model unless their total mass is larger than ~3 M_Jup.

Here we propose that the heavy element enrichment of hot Jupiters can be explained by the accretion of very dust-polluted gas in the inner part of the disk (essentially in situ; Bodenheimer et al. 2000; Batygin et al. 2016). This idea does not contradict the picture proposed by Schneider and Bitsch but introduces a previously overlooked effect. It can explain the enrichment in refractory elements, whereas Schneider and Bitsch predict an enrichment only in volatile elements.

It is generally expected that the gas accreted by a giant planet is metal-poor. This is because dust coagulation converts most of the solid mass into pebble-size objects, which are moderately coupled to the gas; pebbles cannot be accreted by planets exceeding ~ 20–40 Μ_⊕ (the exact value depending on the disk’s viscosity and scale height; Lambrechts et al. 2014; Bitsch et al. 2018) because they remain trapped in the pressure bump produced at the outer edge of the planet-induced gap.

In this manuscript we show that, while this is true for giant planets in the central and outer parts of the disk, at small orbital radii pebbles readily permeate the gap of close-in giant planets due to two important factors: (1) the gas flows through the planet’s gap in the outside-in direction if the planet is sufficiently close to the star, whereas it flows through the gap in the inside-out direction if the planet is farther out (Dürmann & Kley 2015), and (2) pebbles have very low Stokes numbers in the inner part of the disk (Batygin & Morbidelli 2022) and therefore can be entrained into the gap by the radial flow of the gas despite the existence of a pressure bump at the gap’s outer edge. Pebbles are expected to fragment while flowing into the gap, thus maintaining a low Stokes number even inside the gap. The strong coupling with the gas then makes the accretion of solids possible only in conjunction with the accretion of gas. Nevertheless, this can considerably increase the planet’s heavy element budget because of the very high dust/gas density ratio that can be achieved in the inner disk due to the rapid drop in the dust/gas radial velocity ratio.

For simplicity, we elaborate on this process in Sect. 2 by adopting a viscous α-disk model, but we then discuss in Sect. 3 how the process is affected if the radial transport of gas is mostly due to angular momentum removal in disk winds and/or occurs only near the surface of the disk. Thus, we present a comprehensive revision of the problem of the interaction of a gap-opening planet with the flow of gas and dust or pebbles, hopefully correcting some misconceptions that are often found in the literature.

This manuscript ends with a discussion on super-Earths and the case of Jupiter and Saturn in Sect. 4, and we summarize our conclusions in Sect. 5.

Fig. 1 Reproduction of Fig. 7 from Thorngren et al. (2016), showing the mass in heavy elements of weakly irradiated transiting giant planets as a function of the total planet mass. Relationship (1) is also illustrated. The magenta curve shows the expected enrichment if one considers the accretion of gas with 1% metallicity and the accretion of planetesimals following Venturini & Helled (2020). The green curve represents the case with a gas metallicity of 30%, the most favorable case expected from Sect. 2.

2 Pebble dynamics in the vicinity of a giant planet in an α-disk model

The dynamics of a pebble in a protoplanetary disk is dictated by its coupling with the gas. Due to gas drag, the pebble’s velocity, v, relative to the gas velocity, u, is damped following the equation $$\frac{\text{d}\upsilon }{\text{d}t}=-\frac{1}{t_{\text{f}}}\left(\upsilon -u\right),$$(2)

where t_f is called the friction timescale. The Stokes number, S_t, of a pebble is the value of its friction timescale in units of the local orbital timescale: $$S_t=t_{\text{f}}\text{Ω},$$(3)

where Ω is the local orbital frequency.

The radial velocity of a pebble is given by $${\upsilon }_{\text{r}}=\frac{u_{\text{r}}}{1+S_t^2}-2S_t\left({\upsilon }_{\theta }-u_{\theta }\right)$$(4)

(Nakagawa et al. 1986; Takeuchi & Lin 2002), where the r and θ subscripts denote the radial and azimuthal components of the velocities. This formula does not account for the back-reaction of dust over gas (Dipierro et al. 2018). For a pebble on a circular orbit, the difference in azimuthal velocities ν_θ − u_θ is a fraction η (depending on r) of the Keplerian velocity v_K, where $$\eta \left(r\right)=-\frac{1}{2}{\left(\frac{H}{r}\right)}^2\frac{\partial \log P}{\partial \log r},$$(5)

Η is the scale height of the disk, and $P={\left(H\text{Ω}\right)}^2\text{Σ}/\left(\sqrt{2\pi }H\right)$ is the internal pressure of the gas. Assuming H ∝ r (i.e., neglecting disk flaring), Eq. (5) can be approximated as $$\eta \left(r\right)=-\frac{1}{2}{\left(\frac{H}{r}\right)}^2\left[\frac{r}{\text{Σ}}\frac{\partial \text{Σ}}{\partial r}-2\right],$$(6)

which is the equation we use in the rest of this paper.

Equation (6) reveals that, where Σ monotonically decays with r, η > 0. In this case, the azimuthal drag of the gas on the pebble leads to the star-ward radial drift of the pebble. However, where Σ has a sufficiently positive radial gradient, η < 0 and the drag is reversed. The location where η = 0 is called a pressure bump.

Giant planets open deep gaps in the gas distribution of the disk. Along the outer edge of the gap δΣ/∂r is positive and large, so that η < 0, whereas far from the planet’s orbit η > 0. So, a pressure bump is established whenever a giant planet forms. If the radial velocity of the gas u_r > 0 then pebbles cannot drift into the gap, whatever their Stokes number. If instead u_r < 0 only particles with S_t > min[u_r(r)/2η(r)v_K] do not penetrate the gap, where the minimum is computed for r ranging from one to multiple Hill radii beyond the planet’s orbit.

In most hydrodynamical studies on gap opening by giant planets published in the literature, the giant planet is kept on a fixed orbit. The gas is then observed to flow through the gap, from the outer part of the disk to the inner part (Bryden et al. 1999; Kley 1999; Lubow et al. 1999; Weber et al. 2018). Thus, Weber et al. computed that only particles with S_t > 10⁻³ do not penetrate a gap opened by a Jupiter-mass planet, for the disk parameters used in their nominal simulations. Consequently, it is often considered in the literature that the so-called planet barrier against the radial drift of pebbles is effective only for pebbles with a Stokes number larger than this order of magnitude (e.g., Stammler et al. 2023).

However, the situation is radically different if the planet is allowed to migrate, instead of being kept artificially on a fixed orbit. It has been shown (Duffell et al. 2014; Dürmann & Kley 2015; Robert et al. 2018) that, although the migration speed of giant planets is proportional to disk viscosity, giant planets migrate faster than the unperturbed radial velocity of the gas (which, in a viscous, disk is u_r = −3/2(v/r), where ν is the disk’s viscosity, expressed as ν = αΗ²Ω in the so-called α-disks) when $\text{Σ}\left(r_{\text{p}}\right)r_{\text{p}}^2/M_{\text{p}}>0.2$, where Σ(r_p) is the gas’s unperturbed surface density at the radial distance of the planet r_p. In this case, the gas flows through the gap in the inside-out direction (even if the radial motion of the gas in an absolute reference frame can remain negative). Instead, the planet migrates significantly slower than the radial motion of the gas, allowing gas to pass through the gap in the outside-in direction, if $$\frac{\text{Σ}\left(r_{\text{p}}\right)r_{\text{p}}^2}{M_{\text{p}}}\ll 0.2$$(7)

(Dürmann & Kley 2015). Because of the $r_{\text{p}}^2$ dependence of this formula and the typical weak radial decay of Σ (usually proportional to $1/\sqrt{r}$) this happens only when the planet is in the very inner part of the disk. For a typical disk with surface density $175\text{\hspace{0.17em}g\hspace{0.17em}}{\text{cm}}^{-1^2}/\sqrt{\left(r/5.2\text{\hspace{0.17em}au}\right)}$, (delivering a stellar accretion rate of ~4 × 10⁻⁸ M_⊙ yr⁻¹ for α = 3 × 10⁻³) and a Jupiter-mass planet, condition (7) translates to r_p ≪ 2.5 au (Robert et al. 2018), which is well satisfied by the planets studied in Thorngren et al. (2016).

Obviously, in a reference frame co-moving radially with the planet, u_r > 0 in the first case (inside-out flow through the gap) and u_r < 0 in the second case. This means that the planet barrier is effective for particles of all sizes and Stokes numbers when a giant planet migrates in the central or outer parts of the disk², whereas the barrier starts to be leaky when the planet reaches the vicinity of the star.

In the following, we present some quantitative estimates of the Stokes number of pebbles that can penetrate a gap opened by a giant planet, once the latter is close enough to the star for condition (7) to be true.

2.1 Searching for a dust trap: Are particles characterized by a size or a Stokes number?

Before computing the critical Stokes number below which pebbles can penetrate the gap, we need to determine whether particles should be characterized by a given size or a given Stokes number, whatever their position in the disk. In fact, if particles are characterized by a given S_t, we can apply Eq. (4), computing u_r and η as a function of r and setting $S_t^{\text{crit}}=\min \left[u_{\text{r}}\left(r\right)/2\eta \left(r\right){\upsilon }_K\right]$. If instead particles are characterized by a given size, the Stokes number depends on the gas surface density Σ as $S_t=\sqrt{2\pi }\frac{{\rho }_{\text{p}}R}{\text{Σ}}$, where ρ_p is the bulk density of the pebble of radius R; thus, we need to set $S{}_t\left(r\right)=S_t^0\text{Σ}\left(r_0\right)/\text{Σ}\left(r\right)$, where $S_t^0$ is the Stokes number of the pebble at a reference location r₀ and Σ(r) is the actual density of the disk due to the presence of the gap, and solve Eq. (4) with respect to $S_t^0$.

Particles continuously collide with each other, and break or coagulate depending of their collision speed. If the size of particles is limited by the fragmentation barrier, they achieve at collisional equilibrium a Stokes number given by (Birnstiel et al. 2009) $$S_t^{\text{coll,eq}}={\upsilon }_{\text{frag}}^2/\left(3\alpha c_{\text{s}}^2\right),$$(8)

where u_frag is the velocity threshold for fragmentation and c_s is the sound speed. Equation (8) is independent of the gas surface density, Σ, if not indirectly through c_s ∝ Σ^1/8, a weak dependence that we neglect in the following for simplicity³. Equation (8) applies if the timescale over which a particle experiences a change in Σ is longer than the collision timescale with other particles.

In the search for a dust trap at some location along the Σ-gradient characterizing the gap, it is correct to assume that the particle’s Stokes number is that given by Eq. (8). In fact, when a particle is trapped, it has the time to experience collisions and grind down until its Stokes decreases to $S_t^{\text{coll.eq}}$. Thus, the trapping can be permanent only if the condition v_r = 0 in Eq. (4) holds for $S_t=S_t^{\text{coll.eq}}$. Searching for a dust-trap location assuming a fixed dust size is not a valid approximation.

2.2 Evaluating the velocity of the gas

Having determined that S_t is roughly constant and that the critical Stokes number below which particles can penetrate into the gap is $S_t^{\text{crit}}=\min \left[u_{\text{r}}\left(r\right)/2\eta \left(r\right){\upsilon }_K\right]$, we now proceed to evaluate η(r) and u_r(r).

We started from the gap profile formula provided in Crida et al. (2006), which gives an expression for η as a function of Δ = (r − r_p)/R_H, α and q = M_p/M_star: $$\eta =-\frac{1}{2}\left[\frac{0.4q^2{r}_{\text{p}}^4{\left(\frac{1}{R_{\text{h}}\text{Δ}}\right)}^4}{\frac{3}{2}\alpha +\frac{R_{\text{h}}/r_{\text{p}}}{8\text{Δ}}+200\frac{R_{\text{h}}/r_{\text{p}}}{{\text{Δ}}^{10}}}-2{\left(\frac{H}{r}\right)}^2\right],$$(9)

where, with respect to formula (14) in Crida et al., we have retained only the term corresponding to (6), assumed r ~ r_p and retained only the dependence on Δ.

Formula (9) gives a minimum of η of −0.024 at Δ = 2.2 for q = 1 × 10⁻³, α = 3 × 10⁻³ and H/r = 0.05, in good agreement with the nominal hydrodynamical simulation of Weber et al. (2018, −0.03; see their Fig. 2). Instead, for q = 7.6 × 10⁻⁵, α = 10⁻³ and H/r = 0.05, formula (9) gives η = −0.009, whereas Bitsch et al. (2018) found η = 0. It is well known that the model of Crida et al. is not very accurate for planets of moderate mass. Thus, we introduce an empirical correction to formula (9) by dividing the first term in the [.] by 3.5(7.6 × 10⁻⁵/q + 0.2).

Concerning the radial velocity of the gas u_r we remark that, due to the conservation of radial mass flow, one has $u_{\text{r}}=u_{\text{r}}^0{\text{Σ}}^0\left(r\right)/\text{Σ}\left(r\right)$, where $u_{\text{r}}^0=-3/2\alpha {\left(H/r\right)}^2{\upsilon }_K$ is the gas radial velocity in the unperturbed α-disk. We note that this implicitly assumes that the radial migration of the planet is much slower than −3/2α(Η/r)²υ_K, that is to say, that condition (7) holds true. Otherwise, we should use the relative velocity $u_{\text{r}}^0=-3/2\alpha {\left(H/r\right)}^2{\upsilon }_K-{\upsilon }_{\text{r}}^P$, where ${\upsilon }_{\text{r}}^P$ is the radial migration velocity of the planet. Again, if $u_{\text{r}}^0>0$ there is no possibility for particles with any S_t to penetrate into the gap.

To evaluate Σ(r), and then u_r, we turned again to the analytic gap model of Crida et al. (2006). From Eq. (6) one has $$\text{Σ}\left(r\right)={\text{Σ}}^0\left(r\right)+{\displaystyle {\int }_r^{\infty }2\frac{\text{Σ}\left(r\right)}{r}\left(1-\frac{\eta \left(r\right)}{{\left(H/r\right)}^2}\right)\text{\hspace{0.17em}d}r}.$$(10)

Equation (10) is implicit, but it can be solved iteratively, starting from r = ∞ and setting Σ(∞) = Σ⁰(∞). Of course, a physical approximation of infinity is 10R_H or so, which is well beyond the planet’s gap.

Figure 2 shows S_t(r), solution of v_r = 0 in (4) in the nominal case considered by Weber et al. (2018). The critical Stokes number $S_t^{\text{crit}}$ is the minimum of S_t (r) and is ~ 10⁻³, in good agreement with Weber et al. numerical results.

Figure 3 shows a map of the critical Stokes number as a function of q and α for H/r = 5% (top) and 3% (bottom). We remark that the two panels are quite similar, revealing a weak dependence on H/r. It is also worth noticing that, for a given value of α, the critical Stokes number does not monotonically decrease with increasing planet’s mass, as one could naively expect. This is due to the fact that, although the pressure bump becomes stronger with increasing q, −u_r(r) increases at any value of r because the gap becomes wider and deeper. Consequently, the value $S_t^{\text{crit}}$ is achieved farther away from the planet and can result bigger than that computed for a smaller planet. According to Fig. 3, for H/r = 0.05 the planets that oppose the most severe barriers to pebble drift are those of approximately Saturn’s mass (log q ~ − 3.5).

It is clear from Fig. 3 that, even for the least viscous (α ~ 10⁻⁴) and shallow (H/r = 0.03) disks, giant planets’ gaps are leaky for pebbles with S_t ≲ 10⁻⁴. This value is small, but it is the characteristic Stokes number of pebbles that reach the fragmentation limit (8) in the inner part of the protoplanetary disk, where condition (7) is fulfilled (Fig. 4). In particular, when planets approach the very inner part of the disk where the magneto-rotational instability (MRI) is active, the α parameter is expected to increase from 10⁻⁴ toward several times 10⁻³ (Flock et al. 2016). If α ~ 10⁻³, Fig. 4 shows that pebbles have S_t < 10⁻⁴, which is significantly smaller than $S_t^{\text{crit}}$ (Fig. 3).

Thus, we conclude that when giant planets migrate near the central star, eventually pebble isolation disappears and planets can potentially feed from the full radial flow of solid (refractory) material. The efficiency of the accretion process is discussed next.

Fig. 2 Value of the Stokes number of particles with zero radial velocity, as a function of the radial distance beyond the planet, normalized to the planet’s RH. The parameters used for the calculations are the same as those of the nominal simulation in Weber et al. (2018): q = 10−3, α = 3 × 10−3, and H/r = 0.05, so that the curve depicted here can be compared with the red curve in the bottom panel of Fig. 4 of Weber et al.. The critical Stokes number is the minimal value along the curve.

2.3 Once in the gap: Accretion on the planet

Particles with $S_t<S_t^{\text{crit}}$ penetrate into the gap. For these particles, the evolution of the Stokes number is not obvious. In fact, it is important to remember that Eq. (8) only applies if the timescale over which a particle experiences a change in Σ is longer than the timescale of collisions with other particles. Because the particles radial speed increases as they go into the gap as Σ⁰(r)/Σ(r), where Σ⁰(r) is the unperturbed surface density and Σ(r) the actual density, it is possible that the migration timescale becomes shorter, so that particles preserve their sizes and increase in Stokes number. By comparing the timescale of the unperturbed radial motion of particles coupled to the gas (i.e., T_drift = r/u_r = 2/[3α(Η/r)²Ω]) with the particle collision timescale, T_coll = 1/[f_dustΩ], and adopting nominal parameters (α = 10⁻³, (H/r) = 0.05, and f_dust = 0.01), we find that the drift timescale becomes shorter than the collision timescale only in gaps that are at least three orders of magnitude deep, namely for planets significantly more massive than Jupiter.

If the gas is very dust-rich, as argued below (Sect. 2.4), the particle collision timescale becomes even shorter. However, for large values of f_dust the stirring of dust by turbulent diffusion in the gas is reduced by their collective inertia (see Schneider & Wurm 2019 for laboratory experiments and Sect. 4 of Batygin & Morbidelli 2020 for an analytic derivation). This lowers the impact velocity and raises the value of the Stokes number at the fragmentation threshold.

Nevertheless, we consider it unlikely that the Stokes number of particles can increase by orders of magnitude. Thus, particles should remain very coupled to the gas (e.g., S_t ~ 10⁻⁴−10⁻³). The usual formulæ for pebble accretion cannot be applied in our case because they are derived for small planets that do not accrete (nor perturb) the surrounding gas. Instead, the dynamics of gas in the vicinity of a giant planet is very perturbed. Hydrodynamical simulations (e.g., Lambrechts et al. 2019) show that part of the gas entering the planet’s Hill sphere is accreted into a bound atmosphere, while some merely passes through the Hill sphere with a residence timescale typically shorter than the Keplerian orbital period around the planet itself. Given the low Stokes number, particles coming with the accreted flow of gas will also be accreted in the envelope, whereas those carried by the unbound flow will not have enough time to decouple from the background flow and will be eventually transported away. Although solid particles can be accreted in this regime, they are accreted along with the gas. For this reason, it would be misleading refer to this process as pebble-accretion, since this term specifically refers to the selective accretion of dust over gas.

The accretion of dust together with gas may suggest that the overall planet metallicity cannot increase in this process. However, we show below that the gas in the inner part of a planetary disk can be very dust-rich, with ƒ_dust that can approach unity.

Fig. 3 Critical Stokes number as a function of planetary mass, q (normalized to the stellar mass) and viscosity parameter, α, for H/r = 0.05 (top) and H/r = 0.03 (bottom). The dark red color in the top-left corner of the top panel denotes a region of parameter space where particles with any Stokes number penetrate into the gap.

Fig. 4 Reedition of Fig. of Batygin & Morbidelli (2022), re-centered in the region of interest, from 0.01 to 10 au.

2.4 Metallicity of the inner part of the disk

In steady state, if dust drifts radially in the disk at the speed υ_r and gas at the speed u_r, the dust/gas ratio ƒ_dust is given by $$f_{\text{dust}}\left(r\right)=f_{\text{dust}}\left(r_0\right)\frac{u_{\text{r}}}{{\upsilon }_{\text{r}}}\left(r\right)\frac{{\upsilon }_{\text{r}}}{u_{\text{r}}}\left(r_0\right),$$(11)

where r₀ is a reference distance in the disk. Correspondingly, from Eq. (4) wehave $$\frac{{\upsilon }_{\text{r}}}{u_{\text{r}}}=1+2S{}_t\eta \frac{{\upsilon }_K}{u_{\text{r}}}=1+\frac{4}{3}S{}_t\eta {\alpha }^{-1}{\left(\frac{H}{r}\right)}^{-2},$$(12)

where for u_r we have assumed the usual formula for an unflared α-disk.

As a standard assumption, we use H/r = 0.05 and η = 0.003 (Bitsch et al. 2015). Following Flock et al. (2016), we also assume $$\alpha \left(r\right)=\frac{{10}^{-2}-{10}^{-4}}{2}\left[1-\tanh \left(10-\frac{1}{r}\right)\right]+{10}^{-4},$$(13)

where we have adopted a 1/r dependence of the temperature and a critical temperature of T_MRI ~ 1000 K to activate the MRI at 0.1 au. Finally, for S_t we used the value derived from a self-consistent description of an α-disk in Batygin & Morbidelli (2022): $$S_t\left(r\right)={10}^{-4}{r}^{9/10}{\left(\frac{{10}^{-3}}{\alpha }\right)}^{4/5}.$$(14)

Using Eqs. (13) and (14), the resulting dust/gas radial velocity ratio given by Eq. (12) is illustrated in Fig. 5. We recall from Eq. (11) that the dust/gas ratio (i.e., the gas metallicity) is inversely proportional to Eq. (12). Thus, in the inner disk, the metallicity can be ~2.5 times higher than at 1 au and ~6 times higher than at 3 au.

To use this result, we need to anchor the dust/gas ratio somewhere in the disk. In a steady-state scenario, where the dust radial velocity is dominated by the −2S_tηυ_K term in Eq. (4), the dust to gas ratio is (Ida et al. 2016, formula 46) $$f_{\text{dust}}=2\times {10}^{-3}{\left(\frac{S_t}{0.1}\right)}^{-1}{\left(\frac{t}{1\text{\hspace{0.17em}My}}\right)}^{-1/3},$$(15)

where t is the age of the disk. We recognize the inverse dependence on S_t (itself proportional to r^9/10; see Eq. (14)) beyond ~0.3 au in Fig. 5. Formula (15) is valid only until the pebble formation front, r_pƒ ~ 50 au(t/1 My)^7/3 (Ida et al. 2016), reaches the outer edge of the disk. After this event, the dust flux drops and so drops also the dust/gas ratio. Even before that time, there may be obstacles to dust drift, due to gaps opened by distant giant planets (Weber et al. 2018) or the formation of dust-trapping rings due to nonideal magnetohydrodynamic effects (Riols et al. 2020). Until all these limitations become real, formula 15 predicts ƒ_dust = 0.2 at 1 au, where S_t ~ 10⁻³ and t ~ 1 Myr (i.e., enhanced by a factor of 50 with respect to solar metallicity). However, the enhancement has to be limited by planetesimal formation. For S_t > 10⁻³ the streaming instability converts dust into planetesimals when the dust/gas volume density ratio on the midplane is of order unity (Li & Youdin 2021). Depending on the disk viscosity and its ability to stir the vertical distribution of dust (usually encapsulated in the so-called Schmidt number S_c) the vertical scale height of the dust layer can be ~1/10 to ~1/3 that of the gas disk (i.e., for α = 10⁻⁴ and S_c = 10 ans 1, respectively); consequently, when ƒ_dust exceeds ~0.1–0.3 planetesimal formation is expected to convert the dust excess into macroscopic bodies. For this reason, we assume that ƒ_dust at 1 au cannot exceed these values. Therefore, we predict that, in the best case scenario, the ƒ_dust in the region of warm and hot Jupiters, where planetesimal formation is inhibited by the Stokes number being ≪10⁻³, can be 0.3 to 1, that is, enriched in metallicity by a factor of 30 to 100. The dust/gas ratio nevertheless remains smaller than one, which justifies the use of Eq. (4) for the dust velocity, even if it neglects the back-reaction of dust on gas.

The green curve in Fig. 1 shows the mass in heavy elements expected for giant planet accreting most of their atmosphere in the inner disk, from a gas enhanced in metallicity by a factor of 30. We note that we do not expect all planets to lie on this curve; the curve illustrates how metal-rich hot Jupiters can become under the most favorable conditions. Correspondingly, a large scatter of values are expected, depending on the actual metallicity of gas and the fraction of the planet’s atmosphere accreted in the inner disk.

Fig. 5 Dust/gas radial velocity ratio, υr/ur, as given by Eq. (12) as a function of radial distance from the star. The validity of the curve extends to the location of the snow line (2–15 au, depending on the stellar accretion rate for the assumed α = 10−4; see Batygin & Morbidelli 2022), where St (and hence υr/ur) jumps by a factor equal to the square of the ratio of fragmentation velocity thresholds for icy and rocky particles.

3 Wind-driven disks and layered accretion

One-dimensional α-disk models – of the type we have employed thus far – are routinely adopted for their simplicity and the ability to obtain analytical estimates. Their realism, however, is diminished in part by the fact that they do not account for other modes of angular momentum transfer within the system. In this section we address this drawback by describing what we expect to happen if the flow of gas toward the central star ensues due to angular momentum removal in disk winds, and if the transport of gas occurs predominantly on the surface layer of the disk. This is typical of disks where Ohmic diffusion dominates (Lesur 2021), namely in the inner disk region (0.1 < r < 1 au, Armitage 2011), but can also occur in viscous disks, due to the meridional circulation unveiled in Takeuchi & Lin (2002) or if there is a dead zone near the midplane (Bitsch et al. 2014).

We started with wind-driven disks. A main difference with the a-disk considered in the previous section is the radial dependence of the gas radial velocity u_r. As we have seen above, in an α-disk u_r scales with υ_K. Thus, the balance between radial drag and pressure bump (see Eq. (4)) is independent of the distance from the star and only depends on the particles’ Stokes number. Particles penetrate into the gap in the inner part of the disk because their Stokes number is smaller there. In a wind-driven disk, instead, u_r scales as 1/r, which is to say it increases faster than υ_K as the gas approaches the star (Lega et al. 2022). This is because the wind removes angular momentum only from a sufficiently ionized layer and ionization depends on the amount of gas encountered by radiation as it penetrates from the surface of the disk toward the mid-plane. Thus, the vertically integrated column density of gas of the “active” layer of the disk, where the gas flows toward the central star, is independent of r. Then, conservation of mass-flow implies u_r ∝ 1/r. In this case, it is even easier for particles to penetrate into a planet-carved gap in the inner disk, because the radial entrainment due to the gas radial velocity is stronger. For the same reason, the condition for outside-in flow of gas across the planet gap (a necessary condition for particles of any size to reach the planet) is no longer given by Eq. (7) and in general can be fulfilled farther out in the disk (Lega et al. 2022). To this end, Aoyama & Bai (2023) reported an acceleration of the gas radial velocity near the gap due to a concentration of magnetic field lines there. In other words, the scenario envisioned above appears even more favorable for efficient accretion of dust by hot Jupiters if the disk's evolution is driven by magnetized winds.

The fact that the flow of gas toward the star occurs near the surface in wind-driven disks (and also in some two-dimensional α-disks models, with or without dead zone) is not an obstacle to dust accretion. The small Stokes number of particles previously considered (S_t ~ 10⁻⁴) ensures that the dust is well coupled to the gas and uniformly distributed in the vertical direction of the disk, unless the turbulent vertical stirring is pathologically low (α/S_c ≪ 10⁻⁴), which is unlikely, particularly in the inner disk. Under these conditions, the flow of dust from a given disk layer into the gap prompts the vertical redistribution of the dust. This prevents the dust from accumulating indefinitely on the midplane, even if the pressure bump operates there.

This argument can be formalized. To fix ideas, we can imagine a disk of gas where u_r = 0 everywhere but in a near-surface layer, where u_r is large enough to transport the dust in that layer into the gap. The vertically integrated radial mass-flow of dust toward the pressure bump is $${\dot{M}}_r=4\pi r\eta {\upsilon }_K{S}_t{\text{Σd}}_{},$$(16)

where η ~ 3 × 10⁻³ is the value of η beyond the pressure bump and Σ_d = ƒ_dust∑_g is the surface density of the dust. The vertical flow of dust to restore a vertically uniform distribution of the dust/gas ratio is $${\dot{M}}_z=2\pi r\text{Δ}rD{\rho }_{\text{g}}\frac{\partial }{{\partial }_z}\left(\frac{\rho \text{d}}{\rho \text{g}}\right),$$(17)

where ρ_d(z) and ρ_g(z) are the volume densities of dust and gas at height z and ∆r is the width of the ring where the dust tends to be concentrated due to the pressure bump on the mid-plane. The latter is $\text{Δ}r=r\text{Δ}{w}_0\sqrt{\alpha /S_t}$ (Dullemond et al. 2018), where ∆w₀ is the radial width of the pressure bump in normalized units. In Eq. (17) D = αH²Ω is the diffusion coefficient. Because on the surface layer of the disk the dust flows into the gap, ρ_d is reset to ƒ_dustρ_g there, and therefore we can approximate ∂/∂z (ρ_d/ρ_g) ~ 1/H[(ρ_d/ρ_g)_z=0 − ƒ_dust].

In a steady state, Eqs. (16) and (17) have to be equal. Approximating ρ_g ~ Σ_g/(2H), this gives $${\left(\frac{{\rho }_{\text{d}}}{{\rho }_{\text{g}}}\right)}_{z=0}=\left[1+4{\left(\frac{S_t}{\alpha }\right)}^{3/2}\frac{\eta }{\text{Δ}{w}_0}\right]f_{\text{dust}}.$$(18)

The term S_t/α is of order unity, as required for a uniform vertical dust distribution. The ratio η/∆w₀ is typically much smaller than 1/4, ∆w₀ being ~0.1 for a pressure bump induced by a Jupitermass planet (Weber et al. 2018). Thus, the dust/gas ratio on the midplane at the pressure bump is only moderately increased, by a factor of (1 + 4η/∆w₀) < 2 with respect to the local, vertically integrated disk metallicity ƒ_dust. Once this moderate enrichment is achieved, a steady state dust flux is set and all the net dust flow is carried into the gap, preventing any further accumulation of dust at the pressure bump in the midplane. In particular, it is unlikely that planetesimals would start to form in the midplane near a giant planet gap if they could not form in absence of the planet.

4 The case of super-Earths, and of Jupiter and Saturn in the Solar System

If giant planets cannot block the flux of dust in the inner part of the disk, the case is considerably more hopeless for super-Earths, as they open much shallower gaps. This does not imply that close-in super-Earth grew efficiently in situ by pebble accretion. The small Stokes number of the dust and their uniform vertical distribution in the disk makes pebble accretion an inefficient three-dimensional process (Batygin & Morbidelli 2022). Indeed, as we have seen above, even for giant planets the accretion of dust has to occur together with the accretion of gas. Given that super-Earths, by definition, accreted only moderate quantities of gas, we do not expect that the accretion of dust delivered a substantial fraction of a super-Earth's solid mass. Planets with M_p ~ M_h ~ 70 M_⊕, as visible in Fig. 1, are not reproduced in our model starting from a 15 M_⊕ core, as the green curve in the figure shows. Indeed, within the framework of our picture, these objects require the accretion of a large amount of planetesimals or mutual merging of multiple super-Earths of smaller masses.

Jupiter and Saturn are also enriched in heavy elements relative to solar metallicity. As individual planets in the outer disk, condition (7) would not be satisfied because the gas would flow through their gaps in the inside-out direction. However, Jupiter and Saturn have the tendency to lock in a mean motion resonance within the nebula, which halts or reverses their migration direction (Masset & Snellgrove 2001; Morbidelli & Crida 2007; Griveaud et al. 2023). Once this happens, the gas flows in the outside-in direction through their common gap, carrying small-enough dust with it. The typical Stokes number of particles at ~5 au should be much larger than the threshold of 10⁻⁴ − 10⁻³ for transport into the gap (Fig. 3 and Weber et al. 2018), but particle fragmentation at the pressure bump can produce a small-end tail in the particle size distribution with Stokes numbers smaller than this threshold (Stammler et al. 2023). These particles would be accreted by the planets with the same efficiency as gas (i.e., up to 90%; Lubow et al. 1999)⁴. However, it is unlikely that the presence of these particles would have made the metallicity of the accreted gas super-solar given that, as fragments, they do not represent the bulk of the solid mass and the outer disk metallicity is not expected to be substantially enriched, unlike in the inner disk (Sect. 2.4). Thus, we conclude that the enrichment in heavy elements of Jupiter and Saturn is likely due to accretion of planetesimals (Venturini & Helled 2020) and volatile-element vapors (Guillot & Hueso 2006).

5 Conclusions

In this work, we have analyzed the dynamics of particles near the outer edge of a gap opened by a giant planet in the gas radial distribution. The edge of a gap is a pressure bump, but it also enhances the radial velocity of the gas. In the outer part, giant planet Type-II migration toward the star is typically faster than the radial velocity of the disk gas, and slower in the inner part. Thus, in the outer disk, giant planets are effective barriers against the flow of dust of any size because both the positive radial motion of the gas relative to the planet and the pressure bump prevent particle drift into the gap. Instead, in the inner part of the disk the radial flow of gas relative to the planet is in the outside-in direction and can entrain particles with low Stokes numbers into the gap. We find that for the radial entraining to ensue despite the existence of a pressure bump, the particles’ Stokes number has to be lower than 10⁻⁴−10⁻³, depending on planet mass, disk viscosity, and scale height. For particles whose size is limited by the velocity fragmentation threshold, the Stokes number scales inversely with the square of the gas sound speed, and therefore it decreases rapidly in the inner disk. Thus, in the inner disk all conditions are met for typical particles to flow into the gaps opened by giant planets.

As an aside, this implies that the so-called inside-out planet formation model (Chatterjee & Tan 2014) is unlikely to be operational. That model relies on the ability of the first planet (either formed at or migrated to the inner edge of the disk) to create a pressure bump where drifting dust particles accumulate until forming a second planet, and so on. If even a giant planet is not able to block the drift of particles in the innermost part of the disk, this process should not promote the formation of super-Earth systems.

Returning to giant planet metal enrichment, we have shown that particles, once in the gap, can only be accreted by a planet together with the gas because of the low Stokes number. This does not, however, mean that the accreted material has a stellar metallicity. In fact, if the dust is allowed to freely drift toward the star (i.e., it is not blocked farther out by a dynamical barrier), it naturally piles up in the inner disk, enhancing the local metallicity by an order of magnitude or more. For this reason, if hot Jupiters accreted a substantial fraction of their envelope in situ (see Bodenheimer et al. 2000; Batygin et al. 2016, but these works should be updated to account for the high opacity of the dust-rich gas demonstrated in this paper), the large enrichments in heavy elements deduced from measurements of the planet’s mass and radius could have been achieved by this process.

Although we based our analysis on an α-disk model, we have shown that it also holds if the flow of gas in the disk is dominated by angular momentum removal in magnetized winds. We have also shown that, because of the low Stokes number, the transport of gas and dust into the gap in a surface layer of the disk is sufficient to prevent dust pileup at the outer edge of the gap and ensures that, in a steady state, the full dust flux crosses the orbit of the planet. To sum up, this study has explored the multi-faceted relationship between particles, gas flow, and giant planet migration near gap edges, contributing to a better understanding of the conditions under which particles can enter these gaps. Our findings may help inform future research on the heavy element enrichment observed in hot Jupiters.

Acknowledgments

A.M. is grateful for support from the ERC advanced grant HolyEarth N. 101019380 and to Caltech for the visiting professor program that he could benefit from. K. B. is grateful to Caltech’s center for comparative planetology, the David and Lucile Packard Foundation, and the National Science Foundation (grant number: AST 2109276) for their generous support.

References

All Figures

Fig. 1 Reproduction of Fig. 7 from Thorngren et al. (2016), showing the mass in heavy elements of weakly irradiated transiting giant planets as a function of the total planet mass. Relationship (1) is also illustrated. The magenta curve shows the expected enrichment if one considers the accretion of gas with 1% metallicity and the accretion of planetesimals following Venturini & Helled (2020). The green curve represents the case with a gas metallicity of 30%, the most favorable case expected from Sect. 2.

In the text

Fig. 2 Value of the Stokes number of particles with zero radial velocity, as a function of the radial distance beyond the planet, normalized to the planet’s RH. The parameters used for the calculations are the same as those of the nominal simulation in Weber et al. (2018): q = 10−3, α = 3 × 10−3, and H/r = 0.05, so that the curve depicted here can be compared with the red curve in the bottom panel of Fig. 4 of Weber et al.. The critical Stokes number is the minimal value along the curve.

In the text

	Fig. 3 Critical Stokes number as a function of planetary mass, q (normalized to the stellar mass) and viscosity parameter, α, for H/r = 0.05 (top) and H/r = 0.03 (bottom). The dark red color in the top-left corner of the top panel denotes a region of parameter space where particles with any Stokes number penetrate into the gap.
In the text

	Fig. 4 Reedition of Fig. of Batygin & Morbidelli (2022), re-centered in the region of interest, from 0.01 to 10 au.
In the text

Fig. 5 Dust/gas radial velocity ratio, υr/ur, as given by Eq. (12) as a function of radial distance from the star. The validity of the curve extends to the location of the snow line (2–15 au, depending on the stellar accretion rate for the assumed α = 10−4; see Batygin & Morbidelli 2022), where St (and hence υr/ur) jumps by a factor equal to the square of the ratio of fragmentation velocity thresholds for icy and rocky particles.

In the text