© The Authors 2026

1 Introduction

The behaviour of solids within protoplanetary discs is key for understanding how planets form. Traditional models suggest that dust particles, initially sub-micron to micron in size, coagulate into aggregates that grow to ‘pebble’ sizes (e.g. Birnstiel 2024). Growth can be halted by such barriers as the fragmentation barrier (Dominik & Tielens 1997; Dullemond & Dominik 2005; Blum & Wurm 2008), the bouncing barrier (Zsom et al. 2010; Dominik & Dullemond 2024), and the radial-drift barrier (Whipple 1972; Weidenschilling 1977; Brauer et al. 2008). Depending on their material properties, some solids are more strongly coupled to sub-Keplerian gas discs than others. This coupling is the main driving force of planetesimal formation, producing the kilometre-sized precursors of protoplanets. These can form via resonant drag instabilities (Squire & Hopkins 2018; Magnan et al. 2024a), with the streaming instability being the most widely studied in the literature (Youdin & Goodman 2005; Johansen et al. 2007; Magnan et al. 2024b). From this point onwards, mutual collisions between planetesimals further grow the bodies until the efficient phase of pebble accretion sets in (Ormel & Klahr 2010; Johansen & Lambrechts 2017; Ormel 2017). Pebble accretion is extremely efficient mainly for two reasons. First, the marginal coupling of pebbles to the gas greatly enhances the effective accretion radius, as it causes them to spiral inwards towards the protoplanet rather than being scattered. Second, the previously mentioned drift barrier provides a nearly constant supply of pebble-sized particles drifting inwards from the outer disc.

The conventional way to calculate accretion rates is via a particle approach, in an unperturbed¹ gaseous disc (e.g. Visser & Ormel 2016). Pebbles are modelled as moving through the gas disc under the influence of gravity (from the star and planet) and drag. This method has been able to determine the effects of a planet’s eccentricity on growth (Liu & Ormel 2018), a planet’s angular momentum response (Visser et al. 2020), and even the growth timescale of a possible companion (Konijn et al. 2023).

The motion of solids in protoplanetary discs can also be modelled hydrodynamically by treating dust² as a pressureless fluid. Such an approach has been applied to both grid-based methods (e.g. Paardekooper & Mellema 2006; Mignone et al. 2007; Benítez-Llambay & Masset 2016; Stone et al. 2020; Huang & Bai 2022; Lesur et al. 2023) as well as smoothed-particle hydrodynamics (SPH) frameworks (e.g. Laibe & Price 2012; Price & Laibe 2015; Price et al. 2018). In contrast to Lagrangian particle schemes, the grid-based fluid approach avoids sampling limitations and particle shot noise (Genel et al. 2013; Commerçon et al. 2023) and enables momentum exchange between gas and solids to be handled straightforwardly in the solver.

Early investigations of pebble accretion in the context of hydrodynamical gas models perturbed by a planet were conducted by Morbidelli & Nesvorny (2012), using the FARGO code (Masset 2000). Building on this legacy, using the modern FARGO3D code (Benítez-Llambay & Masset 2016), Chrenko et al. (2024) investigated the pebble flows in a gas disc to study the torque on embedded planets in the 2D pebble accretion regime. In their study, 2D multi-fluid simulations were first performed with gas and pebbles. However, the smoothing of the planetary gravitational potential in these simulations prevented pebble accretion from being resolved. To overcome this limitation, the authors introduced a hybrid approach in which the pebbles are evolved as Lagrangian superparticles in a steady-state gaseous background obtained from hydrodynamical simulations. In this setup gas does not evolve simultaneously with the pebbles, and the back-reaction of the pebbles on the gas is therefore neglected. FARGO3D does, however, allow gas and multiple dust fluids to evolve self-consistently as a coupled system (Benítez-Llambay et al. 2019).

To date, most pebble accretion models have generally assumed a monodisperse particle population, corresponding to a single dust size. Recent work has demonstrated that adopting a polydisperse description can significantly affect the early stages of planet formation, for instance in the acoustic resonant drag instability (Paardekooper & Aly 2025a) or, more specifically, the streaming instability (Paardekooper et al. 2020; Matthijsse et al. 2025; Paardekooper & Aly 2025b). Since solids of different sizes drift at different velocities (Weidenschilling 1977), the accretion of material in later stages is also expected to differ substantially between the mono- and polydisperse cases.

Only a handful of studies have directly examined polydisperse pebble accretion (Lyra et al. 2023; Andama et al. 2022). Lyra et al. (2023) developed an analytical framework for polydisperse pebble accretion, finding that the onset of the Bondi (lower-mass) regime occurs at lower core masses than in the monodisperse case, reducing the need for planetesimal collisions to reach the pebble accretion phase. Their results indicate a modest decrease of 3/7 in accretion efficiency in the Hill (highermass) regime. An earlier investigation by Andama et al. (2022), based on a viscous 1D disc model, found that polydisperse Hill growth yields a higher final core mass (M_iso).

In this study, we developed a hydrodynamical, self-consistent, multi-fluid framework that models pebble accretion rather than Lagrangian particles in a static gas-disc. With this method, we investigated how disc evolution, solid-to-gas back-reaction, and a realistic polydisperse pebble population modify the accretion process.

This paper is structured as follows. In Sect. 2, we show how a fluid model for monodisperse pebble accretion works in an unperturbed disc. In Sect. 3, we discuss the perturbed and polydisperse model setup as well as the assumptions made. We then validate the framework with earlier studies on pebble accretion. In Sect. 4, we describe the results of our numerical simulations, first for an evolving the gas disc in a monodisperse setup, followed by an examination of the polydisperse case. In Sect. 5, we discuss the outcome and results, and in Sect. 6, we summarise the implications.

2 A fluid model for pebble accretion in an unperturbed gas disc

2.1 Governing equations: An accreting planet in the disc

Due to mass and momentum conservation, the governing equations of a 2D, monodisperse dusty disc consist of $${\mathrm{\partial }}_t{\mathrm{\Sigma }}_{\mathrm{g}}+\mathrm{\nabla }\cdot \left({\mathrm{\Sigma }}_{\mathrm{g}}{\mathbf{v}}_{\mathrm{g}}\right)=0,$$(1) $${\mathrm{\partial }}_t{\mathbf{v}}_{\mathrm{g}}+({\mathbf{v}}_{\mathrm{g}}\cdot \mathrm{\nabla }){\mathbf{v}}_{\mathrm{g}}=-\mathrm{\nabla }\mathrm{\Phi }-\frac{{\mathrm{\Sigma }}_{\mathrm{d}}}{{\mathrm{\Sigma }}_{\mathrm{g}}}\frac{{\mathbf{v}}_{\mathrm{g}}-{\mathbf{v}}_{\mathrm{d}}}{{\tau }_{\mathrm{s}}}-\frac{\mathrm{\nabla }P}{{\mathrm{\Sigma }}_{\mathrm{g}}}+\frac{\mathrm{\nabla }\cdot \mathsf{T}}{{\mathrm{\Sigma }}_{\mathrm{g}}},$$(2) $${\mathrm{\partial }}_t{\mathrm{\Sigma }}_{\mathrm{d}}+\mathrm{\nabla }\cdot \left({\mathrm{\Sigma }}_{\mathrm{d}}{\mathbf{v}}_{\mathrm{d}}\right)=\mathrm{\nabla }\cdot \left(D{\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\nabla }\left(\frac{{\mathrm{\Sigma }}_{\mathrm{d}}}{{\mathrm{\Sigma }}_{\mathrm{g}}}\right)\right),$$(3) $${\mathrm{\partial }}_t{\mathbf{v}}_{\mathrm{d}}+({\mathbf{v}}_{\mathrm{d}}\cdot \mathrm{\nabla }){\mathbf{v}}_{\mathrm{d}}=-\mathrm{\nabla }\mathrm{\Phi }-\frac{{\mathbf{v}}_{\mathrm{d}}-{\mathbf{v}}_{\mathrm{g}}}{{\tau }_{\mathrm{s}}},$$(4)

where Σ denotes the surface density, v the velocity, τ_s the stopping time, Φ = Φ_* + Φ_pl the gravitational potential (star and planet), ∇P the gas pressure gradient, and the subscripts g and d refer to gas and dust, respectively. We used a locally isothermal equation of state, $P=c_{\mathrm{s}}^2(r){\mathrm{\Sigma }}_{\mathrm{g}}$ and sound speed c_s ∝ r^−1/2 so that the disc has a constant aspect ratio H_g/r, where H is the scale height. The viscous stress tensor, T, is given by $$\mathsf{T}={\mathrm{\Sigma }}_{\mathrm{g}}\nu [\mathrm{\nabla }{\mathbf{v}}_{\mathrm{g}}+{\left(\mathrm{\nabla }{\mathbf{v}}_{\mathrm{g}}\right)}^T-\frac{2}{3}\mathsf{I}\mathrm{\nabla }\cdot {\mathbf{v}}_{\mathrm{g}}],$$(5)

where ν is the kinematic viscosity. The dust continuity equation includes dust diffusion with a diffusion coefficient D, which we set to D = ν throughout this work (i.e. we assume a Schmidt number of unity for simplicity³).

The planet potential is given by $${\mathrm{\Phi }}_{\mathrm{p}\mathrm{l}}=-\frac{G{M}_{\mathrm{p}\mathrm{l}}}{\sqrt{r^2+r_0^2-2r{r}_0\cos (\phi -{\phi }_0)+{\lambda }^2}},$$(6)

where M_pl is the mass of the planet, the subscript 0 denotes the planet’s location, and λ is a gravitational potential smoothing parameter. The reason this smoothing factor is necessary in a fluid approximation is two-fold. First, it eliminates the problem of a (close-to) infinite potential close by the gravitational source. Second, it provides an analogy for a 3D representation in a 2D simulation, i.e. when a parcel of gas is high in the disc yet close by in the (r, φ)-plane. However, the pebble layer is much thinner (H_p ≪ H_g) and for most Stokes numbers considered in this work, remains well below the accretion radius for the planets in our parameter space (M_pl > 1.5 M_⊕).

For both problems, the smoothing factor creates a realistic solution; however, it misrepresents the mechanics close to the planet in a pebble accretion scenario (see Fig. 3 from Chrenko et al. (2024) for reference). Thus, for pebble accretion we want this factor λ to be as small as possible.

Pebble accretion is such a highly efficient phase of growth because the accretion radius R_acc greatly exceeds planetary radius. In the Hill regime, which is the focus of this work, R_acc (Ormel 2017) can be defined as $$R_{\mathrm{a}\mathrm{c}\mathrm{c}}={\left(\frac{G{M}_{\mathrm{p}\mathrm{l}}{\tau }_{\mathrm{s}}}{\mathrm{\Omega }}\right)}^{1/3}\approx R_{\mathrm{H}\mathrm{i}\mathrm{l}\mathrm{l}}{\mathrm{S}\mathrm{t}}^{1/3},$$(7)

where Ω is the angular velocity, R_Hill = a_pl (M_pl/3M_*)^1/3 is the Hill radius (Hill 1878), and a_pl is the orbital distance of the planet.

2.2 Dust streamlines

In this section, we assumed that the gas disc is static, inviscid, and does not feel the planet - a common assumption in pebble accretion studies (e.g. Ormel & Klahr 2010; Lambrechts & Johansen 2012; Visser & Ormel 2016). An axisymmetric steady solution can be constructed where the gas velocity v_g is purely azimuthal. The equilibrium gas velocity can be found by balancing gravity and the radial pressure gradient. This enables a description of the azimuthal gas velocity by relating it to the Keplerian velocity v_K via a non-dimensional value, η: $$v_{\mathrm{g},\phi }=v_{\mathrm{K}}(1-\eta ),$$(8)

with $$\eta =-\frac{1}{2}\frac{c_{\mathrm{s}}^2}{v_{\mathrm{K}}^2}\frac{\mathrm{d}\ln P}{\mathrm{d}\ln r}\approx {\left(\frac{H_{\mathrm{g}}}{r}\right)}^2.$$(9)

Note that η ≪ 1, so that v_g,φ ≈ v_K.

Unlike the gas, the dust does feel the planet’s influence, but we assume this effect to be small. The dust momentum equation (Eq. (4)) is given by $${\mathrm{\partial }}_t{\mathbf{v}}_{\mathrm{d}}+({\mathbf{v}}_{\mathrm{d}}\cdot \mathrm{\nabla }){\mathbf{v}}_{\mathrm{d}}=-\mathrm{\nabla }{\mathrm{\Phi }}_{\ast }-\frac{{\mathbf{v}}_{\mathrm{d}}-{\mathbf{v}}_{\mathrm{K}}}{{\tau }_{\mathrm{s}}}+\frac{{\mathbf{v}}_{\mathrm{g}}-{\mathbf{v}}_{\mathrm{K}}}{{\tau }_{\mathrm{s}}}-\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{p}\mathrm{l}},$$(10)

where the last two terms are treated as small perturbations. As with the gas, we seek a steady state. Without the perturbations, the dust equilibrium velocity is purely azimuthal, with v_d,φ = v_K.

For simplicity, we assumed St ≡ Ωτ_s, the Stokes number, to be constant. The perturbed dust velocity field can be found by linearising the dust momentum equation, Fourier-decomposing the planet potential, and adding all contributions (see Appendix A). The number of Fourier modes necessary for accurate results depends on the smoothing length T; smaller values of λ require more modes. From the resulting velocity field, dust streamlines are calculated.

In absence of the planet, we find that the perturbed velocity field is $$v_r^{\prime }=\frac{2r\mathrm{S}\mathrm{t}({\mathrm{\Omega }}_{\mathrm{g}}-{\mathrm{\Omega }}_{\mathrm{K}})}{1+{\mathrm{S}\mathrm{t}}^2},$$(11) $${\mathrm{\Omega }}^{\prime }=\frac{{\mathrm{\Omega }}_{\mathrm{g}}-{\mathrm{\Omega }}_{\mathrm{K}}}{1+{\mathrm{S}\mathrm{t}}^2},$$(12)

which is the usual radial drift solution in an unperturbed gas disc (Weidenschilling 1977). After including the planet’s perturbation (Appendix A), we integrated to find the velocity field of the pebbles. We show four different velocity fields in Fig. 1 for two different Stokes numbers (left to right) and two different softening parameters λ (top to bottom). The difference in velocity fields for the different Stokes numbers (St) is easily noticeable: a higher St drifts faster while a lower St is more coupled to the gas. The latter also reduces the planet’s influence. The difference for the smoothing factor λ is much less visible; however, there is certainly a difference. We see fewer streamlines accreting with higher λ, which is not unexpected since it essentially lowers the gravitational potential (Eq. (6)). This effect is more pronounced for the lower St since it is more coupled to the gas.

Fig. 1 Velocity fields for pebbles in an unperturbed gaseous disc, calculated by Fourier-decomposing the planet’s potential (as explained in Appendix A). Here, a 10 M⊕ planet is embedded at 1 AU around a 1 M⊙ star. Two Stokes numbers, St = 0.01 (left) and St = 0.1 (right), as well as two different softening parameters, λ = 0.005 (top) and λ = 0.0075 (bottom), are shown. The darker shade indicates the accreted pebbles, the dash-dotted green circle indicates the Hill radius, and the dashed red circle indicates the accretion radius Racc. Only one in five trajectories of the non-accreting pebbles is plotted to avoid overfilling the figure with streamlines.

2.3 Measuring accretion efficiency: The impact of softening parameter

One method for testing pebble accretion involves examining the accretion efficiency ε (Guillot et al. 2014; Lambrechts & Johansen 2014; Liu & Ormel 2018). This efficiency is the fraction of all radially drifting mass accreted by the planet: $$\epsilon =\frac{{\dot{M}}_{\mathrm{a}\mathrm{c}\mathrm{c}}}{{\dot{M}}_{\mathrm{r}\mathrm{a}\mathrm{d},\mathrm{p}\mathrm{e}\mathrm{b}}},$$(13)

where Ṁ_acc is the pebble accretion rate and Ṁ_rad,peb is the total radial mass flow of pebbles in the disc. We calculated this efficiency simply by examining the streamlines (Fig. 1), specifically those that fall into the accretion radius R_acc of the planet (Eq. (7)) and those that do not. We did this for different smoothing parameters λ and for multiple St. We compare our findings to those of Liu & Ormel (2018) in Fig. 2.

Our findings are mostly consistent with the expectations of Liu & Ormel (2018). However, once λ increases, it begins to deviate at the lower end of the Stokes numbers, lowering the efficiency ε. Intuitively, this makes sense since for lower St the decisive ‘accreting moment’ occurs closer to the planet; therefore, it is more heavily impacted by the smoothing factor. This is also visible in Fig. 1, where, for the lower St, fewer streamlines accrete onto the planet at λ = 0.0075 than at λ = 0.005.

Fig. 2 Accretion efficiency, ε, for a 10 M⊕ planet, as calculated by Liu & Ormel (2018) (dashed line), compared to the efficiency obtained with our method for different smoothing parameters, λ (solid-coloured lines).

3 Numerical hydrodynamical model for pebble accretion

We simulated the gas and dust in a global disc using the multi-fluid hydrodynamical FARGO3D code (Benítez-Llambay & Masset 2016; Benítez-Llambay et al. 2019). This is a finite-difference, grid-based hydrodynamical code capable of simulating multiple dust species with an evolving gas disc in multiple geometries. Because we are working in a disc, we used cylindrical coordinates (r,φ, z). Since we primarily examine a phase where 2D accretion takes place, i.e. the accretion radius surpasses the pebble scale height, we used a vertically integrated 2D setup in the (r, φ) grid. Opting for a 2D rather than 3D simulation directly helps keep the computation time manageable.

For the duration of this work, we focus on a 1 M_⊙ star, with a gas surface density at the location of the planet of Σ_g (a_pl = 1 AU) = 10⁻⁴ M_⊕ AU⁻². The disc viscously evolves with an α-description of α = 10⁻³. We obtain a constant aspect ratio of H_g/r = 0.05. The inner and outer radii of our simulation are situated at 0.6 AU and 1.6 AU, respectively. Our grid has dimensions of 400x2512 in the radial and azimuthal direction, respectively, such that the grid cells around the planet are square. Note that the code essentially uses dimensionless code units, but we adopt M_⊙, AU, and (GM_⊙/AU³)^−1/2 as the mass, distance, and time units in all our figures. All our results are therefore scalable with ease, we only keep the aforementioned units to make our figures more easily comprehensible for the reader.

We did not adopt a prescribed or constant pebble mass flux at the radial boundaries. Instead, at both the inner and outer radial boundaries, we applied standard Keplerian disc boundary conditions. In these, the density and azimuthal velocity of each fluid are relaxed to the analytic background disc profiles, while the radial velocity is antisymmetric, corresponding to a no-penetrating condition. These boundaries therefore act as a reservoir that maintains the background disc state near the domain edges without enforcing a fixed inflow rate. The pebble mass flux is thus not externally controlled but emerges self-consistently from the interaction between the disc and the planet. The reasoning behind this method is as follows: Imposing a constant pebble mass flux in a polydisperse disc is non-trivial, since different grain sizes experience distinct drift velocities and filtering efficiencies; therefore, we adopt these reservoir-type boundary conditions to allow the mass flux of each species to adjust self-consistently.

3.1 Accretion mechanism

Previous studies have raised concerns about whether pebble accretion can be accurately reproduced within a fluid approximation. Chrenko et al. (2024) for instance found multifluid simulations to be impractical, mainly due to the smoothing of the planetary gravitational potential. These difficulties are demonstrated in their Fig. 3.

In this work, we adopted an alternative strategy: the gas retains its usual gravitational smoothing (λ_g = 0.6H_g) to account for its vertical scale height, but each dust fluid experiences an unsmoothed potential (λ_d = 0). This is valid because the pebble scale height H_p is mostly smaller than the accretion radius in our parameter space (Birnstiel 2024; Fromang & Nelson 2009; Binkert 2023), implying 2D accretion. We can see this by assuming the scale height of the pebbles to be $$H_{\mathrm{p}}=H_{\mathrm{g}}\sqrt{\frac{\alpha }{\mathrm{S}\mathrm{t}+\alpha }},$$(14)

where the maximum scale height corresponds to our smallest Stokes number grain, i.e. St_min = 10⁻², in our case. This would correspond to a pebble scale height, which remains smaller than the Hill radius for our lowest simulated planet mass of 1.5 M_⊕.

This is consistent with a very thin pebble disc, which does not necessitate smoothing to account for vertical height. For every dust fluid, we created a separate potential field⁴ without λ. After which, we removed all the solid mass within R_acc.

Previous attempts to replicate pebble accretion in fluid simulations (Chrenko et al. 2017; Regály 2020; Pierens 2023) were based either on gas-accretion schemes analogous to those of Kley (1999) or on semi-analytical accretion prescriptions derived from particle-based studies (Liu & Ormel 2018; Ormel & Liu 2018). In contrast, our approach is more targeted towards pebble accretion by removing all material within R_acc, which we argue more accurately captures the essence of pebble accretion since all mass within the accretion radius is captured in that phase by the planet. Since we modelled the pebbles as a pressureless fluid, the streamlines outside R_acc remain unaffected by the accretion process.

3.2 Validating the approach: Comparison to earlier studies

To evaluate our pebble accretion mechanism, we benchmarked it against earlier studies. Liu & Ormel (2018) provides a useful point of comparison. They evaluated the accretion efficiency, ε (Eq. (13)), both locally and globally. Their work employed the conventional particle-based approach in a unperturbed gas disc. They defined the accretion efficiency as the ratio of the number of pebbles accreted over the total number of pebbles entering from the outer disc: $${\epsilon }_{\mathrm{L}\mathrm{O}18}=\frac{N_{\mathrm{h}\mathrm{i}\mathrm{t}}}{N_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}},$$(15)

where N_hit is the amount of pebbles accreted and N_total is the total amount of pebbles reaching the planet from the outer ring. A schematic of this method is shown on the left of Fig. 3.

Since we employed fluid approximation, we required a different method to calculate the accretion efficiency. Due to our chosen disc parameters, the radial mass flow M should be constant in r if no planet is embedded. When an accreting planet is added, the mass removal ensures |Ṁ_in| < |Ṁ_out|, where Ṁ_in and Ṁ_out are the total radial mass flows through rings situated in and outside the planet, respectively. A sketch of this method is shown on the right of Fig. 3. We easily calculate the efficiency from these values via $$\epsilon =\frac{{\dot{M}}_{\mathrm{o}\mathrm{u}\mathrm{t}}-{\dot{M}}_{\mathrm{i}\mathrm{n}}}{{\dot{M}}_{\mathrm{o}\mathrm{u}\mathrm{t}}},$$(16)

where the pebble radial mass flow in a disc is $${\dot{M}}_{\mathrm{r}\mathrm{a}\mathrm{d}}(r)={\int }_{\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{g}}d{\dot{M}}_{\mathrm{r}\mathrm{a}\mathrm{d}}={\int }_0^{2\pi }\mathrm{\Sigma }(r,\phi )v_r(r,\phi )rd\phi$$(17)

and v_r is the radial component of the dust velocity v_dust. We note that Eq. (17) describes the advective (bulk) radial mass flux for pebbles. In addition to this advective component, a diffusive flux arises due to the diffusion term in Eq. (3). In the present work, we do not explicitly analyse the diffusive contribution. In the absence of sharp gradients, the ratio of diffusive to advective flux is ∼α/St, meaning that for our parameters, the advective flux is always more than an order of magnitude larger. Once we reach the pebble isolation mass, small grains partially diffuse across the pressure bump generated by the planet (e.g. Bitsch et al. 2018; Ataiee et al. 2018), thereby modifying the net flux. In this work, however, we do not analyse this effect.

The orbital distances of Ṁ_in and Ṁ_out are located at a_in, and a_out, respectively. We chose a_in and a_out carefully to truly capture the difference between the incoming and passing mass flow. Choosing values too close to the planet for instance could artificially increase or decrease both the mass flow and the efficiency. Doing so via this method, we obtain the same results if we measure efficiency by comparing the removed disc mass with the incoming mass flow.

We compare our simulation to Liu & Ormel (2018)’s (dashed lines) in Fig. 4. For completeness we also show Eq. (33) of Lambrechts & Johansen (2014) (dotted lines). We use an unperturbed gas disc, consistent with the classical particle-based pebble accretion setup. We use a gas pressure gradient of η = 0.001875. The gas is kept unperturbed by allowing only the dust fluids feel the planetary potential (i.e. Φ_pl = 0 for the gas), so the gas remains static while the dust responds to the planet. Back-reaction is disabled to prevent the dust from influencing the gas.

We simulated four different planet masses - 1.5, 4.7, 10, and 30 M_⊕ - for completeness. Note that the latter planetary mass achieves an efficiency of ε > 1 in line with the expectation of Liu & Ormel (2018). We therefore included it mostly as a sanity check. The two studies included here for comparison use an inviscid disc, which is why we also account for inviscid simulations in our comparison. It is clear that the inviscid disc more closely resembles the expected values from Liu & Ormel (2018). Our viscous simulation shows a slight discrepancy for lower Stokes number pebble, which is not unexpected. For St ≲ H²_g/r², the viscous radial gas speed exceeds the pebble drift speed, likely affecting pebble accretion. We note lower efficiencies for the higher Stokes numbers compared to Liu & Ormel (2018). We attribute this to the fact that the power law relies on the low-Stokes approximation (their Appendix A.1), which our simulations do not use.

This comparison demonstrates the promising results of our method for intermediate planetary masses (M_pl ∈̃ [1,10] M_⊕), as well as its necessity for higher masses (M_pl ≳ 10 M_⊕). The latter is because the local particle approach can overestimate the efficiency (exceeding 1), and neither the local nor the global particle approach account for the perturbations of the planet on the gas disc (which we address in Sect. 4.1). We note, however, the shortcomings of our method for lower planetary masses; even for 1.5 M_⊕, the smallest two Stokes number grain sizes are not accreted since R_acc is smaller than the grid-cell size of our simulation. At these masses the gas is insignificantly perturbed (see Sect. 4.1), so the particle-based method is expected to be sufficient.

Fig. 3 Left : sketch of the accretion efficiency determined by Liu & Ormel (2018). Individual pebbles simulated either hit the planet and accrete or pass by the ring at r = apl. The efficiency is calculated using Eq. (15). Right : sketch of the accretion efficiency, as determined by our method. The dust radially drifts inwards with a mass flow Ṁ, calculated at two rings located at ain and aout using Eq. (17). The efficiency ε is subsequently calculated using Eq. (16).

Fig. 4 Accretion efficiency, ε, for multiple planetary masses in a disc where the gas is unperturbed by the planet and back-reaction is disabled. The results are taken from the repo package of Liu & Ormel (2018) (dashed lines); Eq. (33) of Lambrechts & Johansen (2014) (dotted lines); our simulation for a viscous disc with α = 10−3 (dots); and an inviscid disc (crosses). Note that the smallest two Stokes numbers are not plotted for the 1.51 M⊕ planet because the accretion radius for these Stokes numbers is smaller than the cell width at our resolution.

3.3 Approximating the size distribution for a polydisperse disc

In our polydisperse simulations, we considered a continuous size distribution. This is computationally challenging with a finite number of species. We tackled this problem by means of the Gauss-Legendre (GL) quadrature.

Following the methods of Paardekooper et al. (2020), we first defined a ‘size density’ σ, such that $${\mathrm{\Sigma }}_{\mathrm{d}}=\int \sigma (\mathrm{S}\mathrm{t})\mathrm{d}\mathrm{S}\mathrm{t},$$(18)

where Σ_d is the total surface density of the solids. The problem lies in the drag term in Eq. (2) (second term on the right-handside). For a polydisperse setup this becomes an integral, which we approximated as $$\frac{1}{{\mathrm{\Sigma }}_{\mathrm{g}}}\int \sigma \frac{{\mathbf{v}}_{\mathrm{d}}-{\mathbf{v}}_{\mathrm{g}}}{\mathrm{S}\mathrm{t}}\mathrm{d}\mathrm{S}\mathrm{t}\approx \frac{1}{{\mathrm{\Sigma }}_{\mathrm{g}}}\sum _n{w}_n\sigma ({\mathrm{S}\mathrm{t}}_n)\frac{{\mathbf{v}}_{\mathrm{d}}({\mathrm{S}\mathrm{t}}_n)-{\mathbf{v}}_{\mathrm{g}}}{{\mathrm{S}\mathrm{t}}_n},$$(19)

with weights w_n and nodes St_n. These integration nodes and weights (St_n, w_n) correspond to the roots and coefficients of the Legendre polynomial of order n. This approach minimises integration errors compared to simple midpoint or logarithmic binning schemes. Matthijsse et al. (2025) showed a lognormal distribution approximated with an error of ∼10⁻⁴ using only five bins with the GL-method (their Fig. 1). By contrast, this same amount of bins gives an error of ∼10⁻² using a log-uniform method. Furthermore, the growth rate they measured shows higher accuracy for a GL method with five bins, compared to the log-uniform distribution with 40 bins (their Table 2).

Following Matthijsse et al. (2025), we generated the nodes and weights using the scipy.special.roots.legendre routine (Virtanen et al. 2020) and applied them to an MRN-like dust size distribution (Mathis et al. 1977; Draine & Lee 1984) with F_MRN(a) ∝ a^−3.5, corresponding to a mass distribution of ∝a^−0.5 ∝ St^−0.5. We sampled Stokes numbers in the range St ∈ [10⁻²,10⁰]. This distribution is shown in Fig. 5. Note that, due to the logarithmic x-axis and the chosen weights, it might not appear to be an MRN-distribution, which is why we show the dashed line as a visual aid.

If we then define the densities of the individual pebble species as $${\mathrm{\Sigma }}_{\mathrm{d},n}\equiv w_n\sigma ({\mathrm{S}\mathrm{t}}_n),$$(20)

we can use the back-reaction sum in FARGO3D and thereby simulate an analogue to a continuous distribution, albeit technically discrete.

Fig. 5 Distribution of density weights over different St numbers for an MRN-distribution spanning St ∈ [10−2, 100]. The dots represent the calculated St values, and the dotted lines showcase the weights that every dot represents. The dashed line shows the density fraction divided by the weight of the Stokes number, closely following the MRN-distribution (∝a−0.5 ∝ St−0.5).

Fig. 6 Top : accretion efficiency, ε, for three different planet masses for a static (dots) and evolving disc (crosses). Overplotted (dashed lines) are the results of Liu & Ormel (2018), which represent static discs. Bottom : relative difference between εsta, and εevo, the accretion efficiency for a static and evolving disc, respectively. Note that the evolving disc is more efficient for the lower Stokes numbers.

4 Results

We ran our simulation for different scenarios to compare the full effect of an evolving gaseous disc, dust-to-gas feedback, and a polydisperse disc. First, we ran our original monodisperse simulations with and without an evolving gaseous disc (Sect. 4.1) and dust-to-gas feedback (Sect. 4.2). We then introduced a fiducial model for polydisperse accretion with an evolving gas disc in the Hill regime (Sect. 4.3) and compared our findings with previous analytical (polydisperse) results (Sect. 4.4).

4.1 Monodisperse: Perturbing the gas disc

For the classic case of monodisperse pebble accretion, we turned on the evolving gas disc for the same (viscous) simulations as in Fig. 4 to examine its effect on the accretion efficiency. Our result are shown in Fig. 6.

We find that gas evolution actually improves the efficiency for the lower Stokes numbers. This effect is more pronounced for higher planetary mass. Contrary to the lower Stokes numbers, the efficiency decreases at the highest Stokes numbers for an evolving disc. This is illustrated in Fig. 7, where we compare static and evolving discs. The figure shows results for two Stokes numbers, St = 0.1 (left) and St = 0.889 (right), each with six simulations: three planet masses (1.5, 4.7, and 10 M_⊕), for both perturbed and unperturbed gas discs. In the first five rows, we show from top to bottom the azimuthally averaged gas surface density Σ_gas, pebble surface density Σ_pebbles, radial mass flow Ṁ_rad for both pebbles and gas, as well as the dimensionless η-parameter (Eqs. (8)-(9)). The latter was easily calculated using an isothermal setup, and the density and local sound speed are direct outputs of FARGO3D. We subsequently calculated η via Eq. (9), where we used a finite-difference derivative for the pressure P. The colours denote planet mass, while the dotted lines show the static disc simulations.

The bottom row with four streamline plots only show the 10 M_⊕ simulations, where the difference between the static and evolving discs are illustrated by the trajectories of the pebbles. The black arrows in the evolving gas plots denote the direction of the horseshoe orbits. Since the static gas disc does not feel the planet, there are no horseshoe orbits there.

For the lowest simulated planet mass (1.51 M_⊕), the inclusion of an evolving gas disc has a negligible effect on all plotted quantities. This is expected, as such a low-mass planet does not significantly perturb the surrounding gas (i.e. for our simulations of gas-only discs, the isolation mass⁵ lies between 15 and 17.5 M_⊕). However, at 4.70 M_⊕, still well below the pebble isolation mass, the influence of the planet on the disc becomes apparent. A slight deviation is visible in the gas density (top row), accompanied by a modest change in the dust surface density (second row), a measurable change in the η-parameter (fifth row), and a slight change in mass flow at r < a_pl (corresponding to a differing accretion efficiency of the planet).

The latter might be the most interesting effect. For St = 0.1, we see that the change in mass flow, which is equal to the accretion on the planet (Eq. (16)), is higher for the evolving disc compared to the static disc. The planet accretes more due to the gas disc. For higher St, however, we see the opposite effect: the evolving disc decreases the pebble accretion rate. This effect is small (of the order of ∼10%) but measurable.

We believe this effect to be - at least partly - qualitatively due to the horseshoe orbits created by the planets gravity in the gas (Murray & Dermott 1999). Overplotted in the bottom streamline plots of Fig. 7, we illustrate those horseshoe orbits schematically. For the high Stokes numbers, we see that the pebbles on the edge of accretion are pushed outwards, deflecting away from the planet and thus decreasing efficiency. Whereas for the lower Stokes, we see the opposite effect. Pebbles on their respective way back for a second turn (bottom left of the plot) are pushed inwards towards the accreting planet, increasing the efficiency.

Our findings - decreasing efficiency for high St but increasing efficiency for low St - contradicts the results of Kuwahara & Kurokawa (2020a,b), who find low St to be less efficient for perturbed discs. But we note these are different findings. First, our lowest simulated St is 0.011, whereas their low St are <10⁻³. Second, they simulate 3D, local, inviscid discs, whereas our simulations are global, 2D, and for a viscous disc (α = 10⁻³). Therefore, the two studies are hard to compare.

Fig. 7 Two monodisperse setups for St = 0.1 (left) and St = 0.889 (right). Top five rows: Three different planet masses (1.5, 4.7, and 10 M⊕), denoted by colour. The dotted lines signify a static (as opposed to evolving) gas disc. The planet is situated at apl = 1 AU. Back-reaction is disabled for these results. First row : azimuthally averaged gas surface density. Second row : azimuthally averaged pebble surface density. Third row : total radial pebble mass flow, as described in Eq. (17). Fourth row : total radial gas flow. Fifth row : dimensionless η-parameter relating vg,φ and vκ, as in Eqs. (8) and (9). Bottom plots : pebble streamlines for a static (left), and an evolving disc (right) for a 10 M⊕ planet. These pebble lines are indicated by the grey lines in the azimuthally averaged and/or summed graphs above. The dash-dotted circle indicates the Hill sphere. The arrows on the right plots illustrating the evolving gas disc signify the direction of the horseshoe orbits. Pebble streamlines arise in the outer disc, with the darker ones accreting onto the planet. The accreting fraction for an evolving gas disc is higher by approximately 16% for St = 0.1 and lower by approximately 12% for St = 0.889.

4.2 Negligible effect of back-reaction

We assessed the impact of including dust back-reaction onto the gas and found its effect on the pebble accretion efficiency to be negligible compared to the evolving disc with no back-reaction. The primary impact of back-reaction is rather seen in the gas dynamics, where we find that the flow reverses direction for discs with larger grains (St ≳ 0.3), consistent with the results of Dipierro et al. (2018). Within the explored planet-mass range for this work, however, the accretion rate and efficiency remain approximately unchanged, and we therefore omit back-reaction from the previous section on monodisperse evolving discs.

At larger planet masses (M_pl ≳ 17.5 M_⊕), back-reaction begins to impact dust dynamics as well, as the solid-to-gas ratio f_s/g reaches near-unity values. We did incorporate back-reaction, as discussed below, for a polydisperse disc, as it is essential in this case. Without it we would practically simulate multiple monodisperse simulations rather than the coupled polydisperse system we seek to investigate.

4.3 Polydisperse: A fiducial model

For our polydisperse simulations, we used a fiducial model similar to the monodisperse model previously used. We included the back-reaction of the pebbles on the gas, as well as seven pebble ‘fluids’ representing a continuous distribution. The size bins and corresponding weights were carefully chosen with the method described in Sect. 3.3. The values for the chosen disc parameters are shown in Table 1. This fiducial simulation and the simulations hereafter were run for many orbits, until the disc reached a stationary structure.

This fiducial model replicates 2D Hill accretion, which interests us for multiple reasons. First, it enables 2D simulations, which greatly reduce computation time. Second, this regime marks the onset of the planet perturbing the gas disc, which is not captured in the classical particle-approach framework (e.g. Liu & Ormel 2018) and assumes an unperturbed gas disc. Studying this regime allows us to investigate the gradual decline of pebble accretion and the associated disc response in a polydisperse setting.

The result of our fiducial model is plotted in Fig. 8. On top, we plot the surface density of the gas as well as the pebbles. The first row shows the azimuthally averaged surface density of the pebbles in red. Also plotted are the individual pebble species with nodes St_n where we assume $\sigma (\mathrm{S}\mathrm{t})=\sum _n{\sigma }_n\delta (\mathrm{S}\mathrm{t}-{\mathrm{S}\mathrm{t}}_n)$. Already we see a pile-up of pebbles outside the planet similar to the 10 M_s planets in Fig. 7. Pebble density increases at this location, such that the solid-to-gas ratio in the pile-up increases to f_s/g ≈ 0.016.

The second panel shows the radial mass flow of all pebble species (Eq. (17)). Also included here is the gas mass flow. Note that the gas moves outwards in contrast to the gas in Fig. 7. The reason for this discrepancy is that back-reaction is enabled; the gas reacts to the inward moving mass and, therefore, moves outwards. We note this replicates results of Dipierro et al. (2018), who found a reversal of gas radial flow when the mass is dominated by larger, less coupled grains, similar to the distribution in our fiducial model. We highlight the drop in mass flow at a_pl. Similar to the efficiency, this difference is the exact mass accreted by the planet (Ṁ_acc = ∆Ṁ_rad).

The third panel shows the η-parameter, relating v_g,φ and v_K (Eqs. (8)-(9)). We see a dip outside the planet’s orbit at the same location where the pebble pile-up is visible (top panel) as well as a small decrease in radial flow (second panel).

Figure 9 shows the total accretion on the planet decomposed over the different pebble species. This is the same Ṁ_acc denoted in the second panel, as the difference in Ṁ_rad outside and inside the planet signifies the mass accreted onto the planet. The dashed blue lines show the accretion rate for the same simulation for the fiducial model but with a static disc. We see that the larger grains are suppressed by the gas, whereas for the lower St, we observe a small increase due to gas evolution. Overlaid in red is the original MRN-distribution, as seen in Fig. 5. Comparing the two we see the accretion is lob-sided towards higher Stokes numbers. The largest grains are the highest contributors to mass accretion at this stage.

4.4 Comparison to analytical polydisperse pebble accretion

In Fig. 10, we compare the accretion rates as a function of planet mass obtained from our simulations across three panels. The top panel shows the absolute mass accretion rate, Ṁ_acc, for different scenarios.

The solid grey and dashed lines represent the analytical predictions for monodisperse (St = 0.889) and polydisperse (St ∈ [0.01,1]) discs, respectively. These analytical estimations were calculated using the efficiency, ε, taken from Liu & Ormel (2018), assuming the same incoming pebble mass flux as in our fiducial disc model (i.e. without the low Stokes approximation). For the polydisperse analytical estimation, we integrated the efficiency over Stokes for an MRN-distribution, where St ∈ [0.01,1].

Our numerical results for static discs (blue symbols) follow the same scaling as the analytical estimations but lie systematically lower. This offset is consistent with the slightly lower efficiencies obtained in our framework for most Stokes numbers compared to the analytical (see Fig. 4).

The numerical results for evolving discs (red and orange symbols) show the same scaling as the static discs, until isolation sets in at M_pl ≳ 17.5 M_⊕. Before isolation, the results for the evolving disc lie slightly lower than the static values. We attribute this to the negative impact of the evolving gas on the higher Stokes numbers (see Fig. 7).

To further investigate the impact of the size distribution, the middle panel shows the ratio between the polydisperse and monodisperse accretion rates (Ṁ_poly/Ṁ_mono). For a static disc, we find ratios of $${\left(\frac{{\dot{M}}_{\mathrm{p}\mathrm{o}\mathrm{l}\mathrm{y}}}{{\dot{M}}_{\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{o}}}\right)}_{\mathrm{n}\mathrm{u}\mathrm{m}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{l}}\approx 0.78\mathrm{\&}{\left(\frac{{\dot{M}}_{\mathrm{p}\mathrm{o}\mathrm{l}\mathrm{y}}}{{\dot{M}}_{\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{o}}}\right)}_{\mathrm{a}\mathrm{n}\mathrm{a}\mathrm{l}\mathrm{y}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{l}}\approx \mathrm{0.72.}$$(21)

Notably, both these ratios are significantly higher than 3/7, as found by Lyra et al. (2023). As detailed in Appendix B, this discrepancy arises because previous analytical results rely on the low Stokes approximation - which underestimates the pebble flux once St ≳ 0.1 - since it neglects the (1 + St²) term from the drift velocity (Eq. (11)).

The bottom panel highlights the relative impact of gas evolution by plotting the ratio between the evolving and static disc accretion rates (Ṁ_evo/Ṁ_sta). The monodisperse (solid) line shows a clear decrease, when the ratio is below 1. This reflects our earlier finding that perturbed gas flow decreases the efficiency for the larger grains with St ≳ 0.3 (see Figs. 6-7). The polydisperse (dashed) ratio shows values closer to unity. This similarly reflects our findings, as it also contains the smaller grains for which accretion is enhanced by gas perturbations.

Finally, all situations including evolving gas discs illustrate a sharp drop at M_pl ≳ 17.5 M_⊕, marking the onset of pebble isolation (e.g. Lambrechts et al. 2014; Bitsch et al. 2018). In this post-isolation regime, the accretion rate for our polydisperse disc is nearly an order of magnitude higher than the monodisperse case. This is because the polydisperse population includes smaller grains that are capable of diffusing across the planet-generated pressure bump, whereas the inward flux of larger monodisperse pebbles is effectively terminated. The exact mechanics of how the pebble isolation mass is influenced by disc parameters is beyond the scope of this study and will addressed in future work.

Fig. 8 Snapshot after 6000 orbits of the fiducial simulation (Table 1) implementing our multi-fluid pebble accretion mechanism. Top left : gas density of the disc. Top right : total pebble density. First row : azimuthally averaged surface density Σ for the different dust species. Second row : radial mass flow Ṁrad of the gas and the different dust species, as explained in Eq. (17). Note that the gas flows outwards due to the back-reaction of the inward moving pebbles. Third row: dimensionless η-parameter, as explained in Eqs. (8)-(9), relating vg,φ with vK.

Fig. 9 Pebble accretion rate for the 10 M⊕ planet from Fig. 8 (blue solid line), and in the static gas disc scenario (blue dashed line), decomposed over the different pebble species. Overlaid in red is the original MRN-distribution, as depicted earlier in Fig. 5.

Fig. 10 Pebble accretion rates and relative accretion ratios as a function of planet mass. Top panel : mass accretion rate Macc for different scenarios. The blue triangles and lines represent a static gas disc, while the red and orange crosses and lines indicate the evolving gas disc. The grey lines denote our analytical results, and the solid lines denote the monodisperse populations (St=0.889 for the numerical data; St=1 for the analytic baselines). The enlarged orange cross at Mpl = 10 M⊕ denotes the fiducial simulation depicted in Fig. 8. Middle panel : ratio between poly- and monodisperse accretion rates, showing the numerical results for a static and evolving disc, as well as the analytic reference for a static disc. Bottom panel : ratio between the evolving and static accretion rates, for the monodisperse (solid) and polydisperse (dashed) populations, highlighting the impact of gas evolution on accretion efficiency. The sharp decline and increase at Mpl ≳ 17.5 M⊕ in all panels corresponds to the onset of pebble isolation.

5 Discussion

5.1 (Multi-)Fluid prescription versus particle approach

Our results demonstrate that pebble accretion can be modelled accurately using a multi-fluid hydrodynamic framework. The accretion efficiencies we obtain are consistent with those from previous particle-based studies, validating the multi-fluid method as a reliable alternative. This agreement confirms that, despite the conceptual differences between Lagrangian particle tracking and Eulerian fluid descriptions, the essential physics of pebble accretion are well captured in a multi-fluid treatment.

An important advantage of the multi-fluid framework emerges at higher planet masses, close to the isolation mass. As the planetary mass increases, the gas disc becomes increasingly perturbed, altering the local pressure gradient and gas velocity field. In particle-based approaches, such perturbations can complicate the interpretation of pebble trajectories and introduce several numerical and physical challenges. Specifically, implementing gas-particle back-reaction in a Lagrangian framework requires mapping the momentum exchange of discrete particles onto the Eulerian gas grid. This often leads to particle shot noise unless an extremely large number of particles is used (e.g. Peirano et al. 2006). Furthermore, as the planet mass increases and gaps or pile-ups form, particles tend to cluster in high-density regions (e.g. the pressure bump), leaving other parts of the disc undersampled and making it numerically difficult to maintain a converged feedback calculation.

In a polydisperse scenario, these challenges are amplified because each size bin requires a sufficient number of particles to avoid sampling errors, which becomes computationally prohibitive when approximating a continuous distribution. In contrast, the multi-fluid method treats each pebble species as a continuous density field, allowing the momentum exchange and back-reaction to be handled self-consistently within the fluid solver without shot noise, naturally incorporating the feedback between the gas and multiple pebble species.

This makes the multi-fluid approach particularly well suited for studying later stages of growth, where gap formation and strong pressure perturbations modify pebble fluxes. In this regime, the assumption of an unperturbed background disc becomes invalid, and a hydrodynamic treatment becomes essential. Therefore, the framework developed here provides a natural pathway towards investigating pebble isolation mass, where the planet-induced pressure bump halts inward drift. While a detailed study of isolation mass lies beyond the scope of this paper, our results demonstrate that the present method is well positioned to address this problem in future work.

5.2 Change in pebble accretion efficiency for perturbed gas discs

Allowing the gas disc to evolve modifies the accretion efficiency in a systematic, near-linear manner (Fig. 6). For St ≲ 0.3, the efficiency in the evolving disc is higher than in the static case, whereas for St ≳ 0.3 it is reduced. The transition occurs around St ∼ 0.3 and becomes more pronounced with increasing planetary mass. This trend is robust across the explored mass range; therefore, we believe it to reflect a structural change in the flow rather than stochastic variability.

We believe this can be explained, at least in part, qualitatively by examining the streamlines of the pebbles in the evolving models compared to the static discs (bottom plots of Fig. 7). When the gas is allowed to respond to the planet, a horseshoe region develops that is absent in the static prescription. Particles with low Stokes numbers, which remain tightly coupled to the gas, follow these horseshoe trajectories. During the inward leg of the horseshoe turn, material that has already passed (a_pl) from the outer disc is redirected towards the planet’s vicinity. This essentially gives the pebbles a second chance to be accreted. This ‘second-chance’ mechanism mirrors the behaviour of scattered aerodynamically large pebbles, as found by Huang & Ormel (2023). While they focussed on the 3D accretion of larger pebbles (St ≫ 1), our 2D multi-fluid results suggest a qualitatively similar promotive effect for more tightly coupled grains (St ≲ 0.3).

In contrast, particles with higher Stokes numbers are only partially coupled to the gas and do not complete the same turnaround motion. Instead, they are more readily displaced away from the planet’s feeding region during the horseshoe exchange, reducing the effective accretion cross section. While a fully quantitative description is lacking here (it would require a dedicated analysis of particle trajectories and torque balance), the correlation between the emergence of the horseshoe flow and the sign change in (ε_sta – ε_evo)/ε_sta strongly suggests that the modified co-orbital dynamics are responsible for the St-dependent efficiency shift.

5.3 Polydisperse accretion and its effects

The polydisperse treatment reveals that pebble accretion is intrinsically biased towards higher Stokes numbers, even more so than the underlying disc distribution (Fig. 9). While the MRN-distribution already favours larger particles in terms of mass content, the accretion process itself amplifies this asymmetry. This amplification is a direct consequence of the scaling of the accretion radius with the Stokes number; in the Hill regime this scaling follows R_acc ∝ St^1/3. Consequently, the effective distribution of accreted material becomes more top-heavy than the background population in the disc.

A key result of our polydisperse framework is that the ratio between polydisperse and monodisperse accretion (Eq. (21)) is noticeably higher than the 3/7 value previously estimated by Lyra et al. (2023). As detailed in Appendix B, this discrepancy arises because earlier analytical results rely on a low Stokes approximation for the incoming pebble flux. This approximation remains valid for very small grains but fails to capture the correct flux for pebbles with St ≳ 0.1. Since our distribution includes pebbles up to St = 1, the removal of this approximation yields a higher ratio that aligns with our numerical findings.

Finally, the net impact of gas evolution on the total accretion rate is intrinsically linked to the assumed size distribution. In our fiducial model, which adopts an MRN-distribution with St ∈ [10⁻², 1], we find that perturbing the gas disc overall reduces the total pebble accretion rate compared to a static disc. As shown in Sect. 4.1, the perturbed gas flow promotes accretion for small pebbles (St ≲ 0.3) but hinders it for larger pebbles (St ≳ 0.3). Because the MRN-distribution is mass-dominated by these larger Stokes numbers, their reduced efficiency outweighs the gains made by the smaller, more tightly coupled grains. This suggests that the influence of a planet on its own growth rate is highly sensitive to the pebble population. A disc dominated by significantly smaller grains could potentially show a reverse trend, where gas perturbations cause a net increase in the total mass accretion rate.

5.4 Limitations and future work

Several limitations of the present study should be acknowledged. First, our simulations are strictly 2D. While our 2D simulations capture the essential radial and azimuthal dynamics of pebble accretion, they neglect the complex vertical structure of the gas flow around the planet. For low-mass planets, such as those investigated in this study, 3D gas dynamics become essential. Specifically, realistic 3D gas flows are characterised by a strong midplane outflow and complex circulation patterns that are not captured in a vertically integrated 2D framework.

Since small pebbles are more tightly coupled to the gas, their trajectories are heavily influenced by these 3D flow patterns. The characteristic midplane outflow could significantly modify the accretion efficiency for these small grains, potentially altering our conclusions regarding the ’promotive’ effect of the horseshoes (see Fig. 7). Extending this framework to three dimensions is a next step to fully quantify the accretion rates for low-mass protoplanets. This would enable the simultaneous treatment of both the 3D gas flow structures, the vertical settling, and the stratification of the various pebble species.

Second, turbulent diffusion was implemented in a simplified manner. Although this is sufficient for the parameter space explored here, a more detailed turbulence model may become important when studying pebble isolation or gap formation, where subtle changes in pressure gradients can determine whether particles are trapped or continue drifting inwards.

A particularly promising avenue for future work is the study of pebble isolation mass within the multi-fluid framework. Because the method tracks multiple pebble species simultaneously, all coupled to the gas and therefore each other, it provides direct insight into how different Stokes numbers respond to the formation of a pressure bump as well as how the pressure bump reacts to the subsequent dust traffic jam. This allows one to determine not only when accretion halts, but also which size ranges are most efficiently filtered. Such a self-consistent treatment is difficult to achieve in simplified or monodisperse models.

6 Conclusions

• We demonstrate that polydisperse pebble accretion can be accurately modelled using a multi-fluid hydrodynamic framework. This approach yields results consistent with classical particle-based studies, provided that the planetary mass is sufficient for the accretion radius R_acc to be well-resolved by the numerical grid cells of the dust fluid;

• Perturbing the gas disc modifies the accretion efficiency. For St ≲ 0.3, the efficiency increases relative to the staticdisc case, whereas for St ≳ 0.3 it decreases. This trend becomes more pronounced with increasing planetary mass. We emphasize that this conclusion is based on 2D simulations. A realistic 3D gas flow would significantly influence the pebble accretion efficiency for small pebbles;

• In our assumption of an MRN-distribution with St ∈ [10⁻², 1], evolving the gas disc lowers the total accretion rate, since the mass distribution is dominated by pebbles with St ≳ 0.3;

• We find that back-reaction remains negligible for the planet masses studied here (≤ 10 M_⊕), as the solid-to-gas mass ratio f_s/g never exceeds ∼0.016. Consequently, the accretion rate differs only by a few percent;

• Pebble accretion in a polydisperse disc with an MRN-distribution is intrinsically biased towards higher Stokes numbers. This bias arises from two factors: first, the MRN-distribution is inherently mass-dominated by larger grains, and second, the accretion process itself amplifies this asymmetry because the accretion radius scales with the Stokes number, providing larger pebbles with a significantly greater accretion cross section;

• We find the ratio between the poly- and monodisperse accretion rates to increase noticeably above previous estimations of (Ṁ_poly/Ṁ_mono) = 3/7. This increase occurs once St_max ≳ 0.1. For this ratio we find an analytical value of ∼0.72 and a value of ∼0.78 for our simulations, both for an MRN-distribution where St ∈ [10⁻²,1].

Data availability

The data that support the findings of this study are openly available in 4TU.ResearchData under the name: ‘Repository supporting the publication: A multi-fluid approach for pebble accretion’ at https://doi.org/10.1051/0004-6361/202558434.

Acknowledgements

The authors thank the anonymous referee for helpful suggestions that greatly improved this manuscript. We also thank Jip Matthijsse, Benjamin Silk, and Hossam Aly for their input and helpful discussions. TJK further thanks Michiel Lambrechts, Anders Johansen, Ayumu Kuwahara, Wladimir Lyra & Satoshi Okuzumi for their fruitful discussions while on visit in Copenhagen. The authors acknowledge the use of computational resources of the DelftBlue supercomputer, provided by Delft High Performance Computing Centre (DHPC). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon Europe research and innovation programme (Grant Agreement No. 101054502). This work made use of several open-source software packages. We acknowledge FARGO3D (Benítez-Llambay & Masset 2016), numpy (Harris et al. 2020), matplotlib (Hunter 2007), and scipy (Virtanen et al. 2020).

References

Appendix A Dust streamline calculations

The velocity fields in Figs. 1 & 7 use the linear perturbation theory. The derivation of how this is done is shown below.

We start from the dust momentum conservation Eq. 4, in the monodisperse limit: $$\begin{array}{rl}{\mathrm{\partial }}_t\mathbf{v}+(\mathbf{v}\cdot \mathrm{\nabla })\mathbf{v}=& -\mathrm{\nabla }{\mathrm{\Phi }}_{\ast }-\frac{\mathbf{v}-{\mathbf{v}}_{\mathrm{K}}}{{\tau }_{\mathrm{s}}}+\frac{{\mathbf{v}}_{\mathrm{g}}-{\mathbf{v}}_{\mathrm{K}}}{{\tau }_{\mathrm{s}}}-\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{p}}\\ \equiv & -\mathrm{\nabla }{\mathrm{\Phi }}_{\ast }-\frac{\mathbf{v}-{\mathbf{v}}_{\mathrm{K}}}{{\tau }_{\mathrm{s}}}+{\mathrm{\Pi }}^{\prime }-\mathrm{\nabla }{\mathrm{\Phi }}^{\prime },\end{array}$$(A.1)

where we omitted the subscript ’d’ for clarity. The gas disc is taken to be static, with a purely azimuthal velocity v_gφ = (1 – η)v_K. For η ≪ St, and a sufficiently small planet mass, we can treat Π′ and Φ′ as small perturbations. It is straightforward to see that when Π′ = Φ′ = 0, the equilibrium velocity field is v = v_K.

We can describe linearised equations in terms of the angular momentum perturbation $L^{\prime }=r{v}_{\phi }^{\prime }$: $$({\mathrm{\partial }}_t+\frac{v_{\phi }}{r}{\mathrm{\partial }}_{\phi })v_r^{\prime }-\frac{2v_{\phi }}{r^2}{L}^{\prime }=-{\mathrm{\partial }}_r{\mathrm{\Phi }}^{\prime }-\frac{v_r^{\prime }}{{\tau }_{\mathrm{s}}},$$(A.2) $$({\mathrm{\partial }}_t+\frac{v_{\phi }}{r}{\mathrm{\partial }}_{\phi })L^{\prime }+\frac{v_{\phi }}{2}{v}_r^{\prime }=r{\mathrm{\Pi }}^{\prime }-{\mathrm{\partial }}_{\phi }{\mathrm{\Phi }}^{\prime }-\frac{L^{\prime }}{{\tau }_{\mathrm{s}}}.$$(A.3)

A.1 Fourier decomposition

As is standard in disc-planet interactions (Goldreich & Tremaine 1980), we assume the perturbations to be of the form: $$X^{\prime }(r,\phi ,t)={\hat{X}}_m(r)\exp (im\phi -im{\mathrm{\Omega }}_{0}t),$$(A.4)

where m is the (integer) azimuthal wave number and Ω₀ is the angular velocity of the planet. By choosing this form of a perturbation, they are steady in the frame co-rotating with the planet. Implementing this in Eqs. A.2-A.3: $$\sigma {\hat{v}}_m-\frac{2v_{\phi }}{r^2}{\hat{L}}_m=-{\mathrm{\partial }}_r{\hat{\mathrm{\Phi }}}_m,$$(A.5) $$\sigma {\hat{L}}_m+\frac{v_{\phi }}{2}{\hat{v}}_m={\delta }_{m0}r\hat{\mathrm{\Pi }}-im{\hat{\mathrm{\Phi }}}_m,$$(A.6)

where σ ≡ im(Ω - Ω₀) + 1/τ_s and Ω = v_φ/r denotes the angular velocity. This system can be straightforwardly solved: $${\hat{v}}_m=-\frac{1}{r}\frac{2\mathrm{\Omega }im{\hat{\mathrm{\Phi }}}_m-2\mathrm{\Omega }{\delta }_{m0}r\hat{\mathrm{\Pi }}+\sigma r{\mathrm{\partial }}_r{\hat{\mathrm{\Phi }}}_m}{{\sigma }^2+{\mathrm{\Omega }}^2},$$(A.7) $${\hat{L}}_m=-\frac{\sigma im{\hat{\mathrm{\Phi }}}_m-\sigma {\delta }_{m0}r\hat{\mathrm{\Pi }}+\frac{1}{2}\mathrm{\Omega }r{\mathrm{\partial }}_r{\hat{\mathrm{\Phi }}}_m}{{\sigma }^2+{\mathrm{\Omega }}^2}.$$(A.8)

From here we can describe the total velocity and angular momentum perturbations as: $$v_r^{\prime }=\sum _{m=0}^{\mathrm{\infty }}{\hat{v}}_m(r)\exp (im\phi -im{\mathrm{\Omega }}_{0}t),$$(A.9) $$L^{\prime }=\sum _{m=0}^{\mathrm{\infty }}{\hat{L}}_m(r)\exp (im\phi -im{\mathrm{\Omega }}_{0}t).$$(A.10)

In the absence of a planet, we only retain m = 0 and find: $$v_r^{\prime }=\frac{2r\mathrm{\Omega }{\tau }_{\mathrm{s}}({\mathrm{\Omega }}_{\mathrm{g}}-{\mathrm{\Omega }}_{\mathrm{K}})}{1+{\mathrm{\Omega }}^2{\tau }_{\mathrm{s}}^2},$$(A.11) $${\mathrm{\Omega }}^{\prime }=\frac{{\mathrm{\Omega }}_{\mathrm{g}}-{\mathrm{\Omega }}_{\mathrm{K}}}{1+{\mathrm{\Omega }}^2{\tau }_{\mathrm{s}}^2}.$$(A.12)

These describe the usual radial and azimuthal dust drift. For Φ′ ≠ 0, the number of terms required in the Fourier series depends on the smoothing length; smaller values of λ require additional terms to achieve a converged solution.

Fig. A.1 Example of streamlines created using this method. Illustrated is a local co-rating shearing box of a 10 M⊕ planet.

A.2 Series in Stokes number

For dust particles with Stokes number St = Ωτ_s ≪ 1, we can use an alternative to the Fourier decomposition, which can be an advantage for very small smoothing lengths. Combine the steady state (in the frame corotating with the planet) versions of Eqs. A.2 and A.3: $$\begin{array}{r}(\mathrm{\Omega }-{\mathrm{\Omega }}_0{)}^2{\mathrm{\partial }}_{\phi }^2{u}^{\prime }+\frac{2(\mathrm{\Omega }-{\mathrm{\Omega }}_0)}{\tau }{\mathrm{\partial }}_{\phi }{u}^{\prime }+u^{\prime }({\mathrm{\Omega }}^2+\frac{1}{{\tau }^2})=\\ 2\mathrm{\Omega }{\mathrm{\Pi }}^{\prime }-\frac{2\mathrm{\Omega }}{r}{\mathrm{\partial }}_{\phi }{\mathrm{\Phi }}^{\prime }-\frac{{\mathrm{\partial }}_r{\mathrm{\Phi }}^{\prime }}{\tau }-(\mathrm{\Omega }-{\mathrm{\Omega }}_0){\mathrm{\partial }}_{r,\phi }^2{\mathrm{\Phi }}^{\prime }.\end{array}$$(A.13)

Since the only derivatives of the unknown function u′ are with respect to φ, we can treat this as an ODE. While it is possible to write down the general solution, it is of not much use in practice. Assuming the Stokes number to be constant, we can however develop a series in $\mathrm{S}\mathrm{t}\ll 1,u^{\prime }=\sum _n{u}_n^{\prime }{\mathrm{S}\mathrm{t}}^n$, with $$u_0^{\prime }=0,$$(A.14) $$u_1^{\prime }=-2\eta r\mathrm{\Omega }-\frac{{\mathrm{\partial }}_r{\mathrm{\Phi }}^{\prime }}{\mathrm{\Omega }},$$(A.15) $$u_2^{\prime }=-\frac{2}{r\mathrm{\Omega }}{\mathrm{\partial }}_{\phi }{\mathrm{\Phi }}^{\prime }+(\mathrm{\Omega }-{\mathrm{\Omega }}_0)\frac{{\mathrm{\partial }}_{r,\varphi }^2{\mathrm{\Phi }}^{\prime }}{{\mathrm{\Omega }}^2},$$(A.16) $$u_n^{\prime }=-u_{n-2}^{\prime }-\frac{2(\mathrm{\Omega }-{\mathrm{\Omega }}_0)}{\mathrm{\Omega }}{\mathrm{\partial }}_{\phi }{u}_{n-1}^{\prime }-\frac{(\mathrm{\Omega }-{\mathrm{\Omega }}_0{)}^2}{{\mathrm{\Omega }}^2}{\mathrm{\partial }}_{\phi }^2{u}_{n-2}^{\prime },$$(A.17) $$v_0^{\prime }=-\eta r\mathrm{\Omega },$$(A.18) $$v_1^{\prime }=-\frac{1}{r\mathrm{\Omega }}{\mathrm{\partial }}_{\phi }{\mathrm{\Phi }}^{\prime },$$(A.19) $$v_n^{\prime }=-\frac{u_{n-1}^{\prime }}{2}-\frac{(\mathrm{\Omega }-{\mathrm{\Omega }}_0{)}^2}{2{\mathrm{\Omega }}^2}{\mathrm{\partial }}_{\phi }^2{u}_{n-1}^{\prime }.$$(A.20)

See Fig. A.1 as an example for these calculated streamlines.

Appendix B Analytical ratio between polydisperse and monodisperse accretion rate

We find the (analytical) ratio between poly- and monodisperse accretion to be ∼0.72 (Eq. 21), while previous analytical results found an exact factor of 3/7 for this ratio (Lyra et al. 2023). This appendix gives a short explanation of why we expect a higher ratio, due only to the low Stokes number approximation used previously.

The monodisperse accretion rate can be written as the total radial mass flow of pebbles times the accretion efficiency: $${\dot{M}}_{\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{o}}=2\pi r{v}_r{\mathrm{\Sigma }}_{\mathrm{p}}\epsilon ({\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{o}}),$$(B.1)

where v_r is the radial drifting velocity of the pebbles: $$v_r=\frac{2\mathrm{S}\mathrm{t}}{1+{\mathrm{S}\mathrm{t}}^2}\eta v_{\mathrm{K}},$$(B.2)

which is the usual drift solution (Weidenschilling 1977), similar to Eq. 11. We can calculate the polydisperse accretion rate in a similar manner. We just take the size density σ(St) defined in Eq. 18 to find: $${\dot{M}}_{\mathrm{p}\mathrm{o}\mathrm{l}\mathrm{y}}=4\pi r\eta v_{\mathrm{K}}{\int }_{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}\frac{\sigma (\mathrm{S}\mathrm{t})\mathrm{S}\mathrm{t}\epsilon (\mathrm{S}\mathrm{t})}{1+{\mathrm{S}\mathrm{t}}^2}\mathrm{d}\mathrm{S}\mathrm{t}.$$(B.3)

The ratio between these two accretion rates is then: $$\frac{{\dot{M}}_{\mathrm{p}\mathrm{o}\mathrm{l}\mathrm{y}}}{{\dot{M}}_{\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{o}}}=\frac{1}{{\mathrm{\Sigma }}_{\mathrm{p}}}\frac{1+{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{o}}^2}{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{o}}\epsilon ({\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{o}})}{\int }_{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}\frac{\sigma (\mathrm{S}\mathrm{t})\mathrm{S}\mathrm{t}\epsilon (\mathrm{S}\mathrm{t})}{1+{\mathrm{S}\mathrm{t}}^2}\mathrm{d}\mathrm{S}\mathrm{t}.$$(B.4)

We assume ε and σ to be power laws, so we can take ε ∝ St^a and σ ∝ St^b to find: $$\frac{{\dot{M}}_{\mathrm{p}\mathrm{o}\mathrm{l}\mathrm{y}}}{{\dot{M}}_{\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{o}}}=\frac{1+b}{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{b+1}-{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{i}\mathrm{n}}^{b+1}}\frac{1+{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{o}}^2}{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{o}}^{1+a}}{\int }_{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}\frac{{\mathrm{S}\mathrm{t}}^{1+a+b}}{1+{\mathrm{S}\mathrm{t}}^2}\mathrm{d}\mathrm{S}\mathrm{t},$$(B.5)

since ${\mathrm{\Sigma }}_{\mathrm{P}}=\int \sigma (\mathrm{S}\mathrm{t})\mathrm{d}\mathrm{S}\mathrm{t}$.

In our assumption for an MRN-distribution (b = −1/2) and taking the St-dependency for the efficiency from Liu & Ormel (2018) and Lambrechts & Johansen (2014) (a = −1/3), this means: $$\frac{{\dot{M}}_{\mathrm{p}\mathrm{o}\mathrm{l}\mathrm{y}}}{{\dot{M}}_{\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{o}}}=\frac{1/2}{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{1/2}-{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{i}\mathrm{n}}^{1/2}}\frac{1+{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{a}\mathrm{x}}^2}{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{2/3}}{\int }_{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}\frac{{\mathrm{S}\mathrm{t}}^{1/6}}{1+{\mathrm{S}\mathrm{t}}^2}\mathrm{d}\mathrm{S}\mathrm{t},$$(B.6)

assuming St_mono = St_max.

For [St_min, St_max] = [0,1], this ratio would be approximately ∼0.66. Our result of Eq. 21 is higher, because we have a finite St_min. If we put in our distribution of St ∈ [10⁻²,1], the ratio becomes ∼0.72, as is also portrayed by the grey line in the middle panel of Fig. 10.

If we use the small Stokes approximation thereby neglecting the (1 + St²) term from Eq. B.6 we find: $$\begin{gathered} \frac{{\dot{M}}_{\mathrm{p}\mathrm{o}\mathrm{l}\mathrm{y}}}{{\dot{M}}_{\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{o}}}=\frac{1/2}{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{1/2}-{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{i}\mathrm{n}}^{1/2}}\frac{1}{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{2/3}}{\int }_{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\mathrm{S}\mathrm{t}}^{1/6}\mathrm{d}\mathrm{S}\mathrm{t}, \\ =\frac{3}{7}\frac{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{7/6}-{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{i}\mathrm{n}}^{7/6}}{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{1/2}-{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{i}\mathrm{n}}^{1/2}}\frac{1}{{\mathrm{S}\mathrm{t}}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{2/3}}. \end{gathered}$$(B.7)

Which for St_min = 0 is exactly 3/7, the original solution of Lyra et al. (2023).

We plotted the effect of St_max on this ratio in Fig. B.1. The blue line keeps the theoretical value of St_min = 0, while we plot the orange line with St_min = 0.01. The vertical dash-dotted line is a visual aid of St_max = 1, our maximum value. We also plot two horizontal dashed lines as a visual aid. The top line denotes our analytical value from Eq. 21, while the bottom line denotes the analytical value of 3/7 found by Lyra et al. (2023).

Fig. B.1 Ratio between the poly- and monodisperse accretion rates (Eq. B.6) for two different Stmin. The horizontal dashed lines signify (Ṁpoly/Ṁmono) = 0.72, our analytical result from Eq. 21, and (Ṁpoly/Ṁmono) = 3/7, the result of Lyra et al. (2023). The vertical dash-dotted line is a visual aid for Stmax = 1, our considered maximum value.

All Tables

All Figures

Fig. 1 Velocity fields for pebbles in an unperturbed gaseous disc, calculated by Fourier-decomposing the planet’s potential (as explained in Appendix A). Here, a 10 M⊕ planet is embedded at 1 AU around a 1 M⊙ star. Two Stokes numbers, St = 0.01 (left) and St = 0.1 (right), as well as two different softening parameters, λ = 0.005 (top) and λ = 0.0075 (bottom), are shown. The darker shade indicates the accreted pebbles, the dash-dotted green circle indicates the Hill radius, and the dashed red circle indicates the accretion radius Racc. Only one in five trajectories of the non-accreting pebbles is plotted to avoid overfilling the figure with streamlines.

In the text

	Fig. 2 Accretion efficiency, ε, for a 10 M⊕ planet, as calculated by Liu & Ormel (2018) (dashed line), compared to the efficiency obtained with our method for different smoothing parameters, λ (solid-coloured lines).
In the text

Fig. 3 Left : sketch of the accretion efficiency determined by Liu & Ormel (2018). Individual pebbles simulated either hit the planet and accrete or pass by the ring at r = apl. The efficiency is calculated using Eq. (15). Right : sketch of the accretion efficiency, as determined by our method. The dust radially drifts inwards with a mass flow Ṁ, calculated at two rings located at ain and aout using Eq. (17). The efficiency ε is subsequently calculated using Eq. (16).

In the text

Fig. 4 Accretion efficiency, ε, for multiple planetary masses in a disc where the gas is unperturbed by the planet and back-reaction is disabled. The results are taken from the repo package of Liu & Ormel (2018) (dashed lines); Eq. (33) of Lambrechts & Johansen (2014) (dotted lines); our simulation for a viscous disc with α = 10−3 (dots); and an inviscid disc (crosses). Note that the smallest two Stokes numbers are not plotted for the 1.51 M⊕ planet because the accretion radius for these Stokes numbers is smaller than the cell width at our resolution.

In the text

	Fig. 5 Distribution of density weights over different St numbers for an MRN-distribution spanning St ∈ [10−2, 100]. The dots represent the calculated St values, and the dotted lines showcase the weights that every dot represents. The dashed line shows the density fraction divided by the weight of the Stokes number, closely following the MRN-distribution (∝a−0.5 ∝ St−0.5).
In the text

Fig. 6 Top : accretion efficiency, ε, for three different planet masses for a static (dots) and evolving disc (crosses). Overplotted (dashed lines) are the results of Liu & Ormel (2018), which represent static discs. Bottom : relative difference between εsta, and εevo, the accretion efficiency for a static and evolving disc, respectively. Note that the evolving disc is more efficient for the lower Stokes numbers.

In the text

Fig. 7 Two monodisperse setups for St = 0.1 (left) and St = 0.889 (right). Top five rows: Three different planet masses (1.5, 4.7, and 10 M⊕), denoted by colour. The dotted lines signify a static (as opposed to evolving) gas disc. The planet is situated at apl = 1 AU. Back-reaction is disabled for these results. First row : azimuthally averaged gas surface density. Second row : azimuthally averaged pebble surface density. Third row : total radial pebble mass flow, as described in Eq. (17). Fourth row : total radial gas flow. Fifth row : dimensionless η-parameter relating vg,φ and vκ, as in Eqs. (8) and (9). Bottom plots : pebble streamlines for a static (left), and an evolving disc (right) for a 10 M⊕ planet. These pebble lines are indicated by the grey lines in the azimuthally averaged and/or summed graphs above. The dash-dotted circle indicates the Hill sphere. The arrows on the right plots illustrating the evolving gas disc signify the direction of the horseshoe orbits. Pebble streamlines arise in the outer disc, with the darker ones accreting onto the planet. The accreting fraction for an evolving gas disc is higher by approximately 16% for St = 0.1 and lower by approximately 12% for St = 0.889.

In the text

Fig. 8 Snapshot after 6000 orbits of the fiducial simulation (Table 1) implementing our multi-fluid pebble accretion mechanism. Top left : gas density of the disc. Top right : total pebble density. First row : azimuthally averaged surface density Σ for the different dust species. Second row : radial mass flow Ṁrad of the gas and the different dust species, as explained in Eq. (17). Note that the gas flows outwards due to the back-reaction of the inward moving pebbles. Third row: dimensionless η-parameter, as explained in Eqs. (8)-(9), relating vg,φ with vK.

In the text

	Fig. 9 Pebble accretion rate for the 10 M⊕ planet from Fig. 8 (blue solid line), and in the static gas disc scenario (blue dashed line), decomposed over the different pebble species. Overlaid in red is the original MRN-distribution, as depicted earlier in Fig. 5.
In the text

Fig. 10 Pebble accretion rates and relative accretion ratios as a function of planet mass. Top panel : mass accretion rate Macc for different scenarios. The blue triangles and lines represent a static gas disc, while the red and orange crosses and lines indicate the evolving gas disc. The grey lines denote our analytical results, and the solid lines denote the monodisperse populations (St=0.889 for the numerical data; St=1 for the analytic baselines). The enlarged orange cross at Mpl = 10 M⊕ denotes the fiducial simulation depicted in Fig. 8. Middle panel : ratio between poly- and monodisperse accretion rates, showing the numerical results for a static and evolving disc, as well as the analytic reference for a static disc. Bottom panel : ratio between the evolving and static accretion rates, for the monodisperse (solid) and polydisperse (dashed) populations, highlighting the impact of gas evolution on accretion efficiency. The sharp decline and increase at Mpl ≳ 17.5 M⊕ in all panels corresponds to the onset of pebble isolation.

In the text

	Fig. A.1 Example of streamlines created using this method. Illustrated is a local co-rating shearing box of a 10 M⊕ planet.
In the text

	Fig. B.1 Ratio between the poly- and monodisperse accretion rates (Eq. B.6) for two different Stmin. The horizontal dashed lines signify (Ṁpoly/Ṁmono) = 0.72, our analytical result from Eq. 21, and (Ṁpoly/Ṁmono) = 3/7, the result of Lyra et al. (2023). The vertical dash-dotted line is a visual aid for Stmax = 1, our considered maximum value.
In the text