© The Authors 2024

1 Introduction

Protoplanetary disks serve as the birthplaces of planets. Spatially resolved observations of these disks have the potential to unravel the physical processes that govern planet formation (see the comprehensive review by Bae et al. 2023). The Atacama Large Millimeter/submillimeter Array (ALMA) has made significant strides in this field by detecting numerous protoplanetary disks exhibiting substructures such as gaps, rings, spirals, and cavities (e.g., Andrews et al. 2018; Zhang et al. 2018; Long et al. 2018). While the origins of these substructures remain elusive, they could plausibly be the result of planet-disk interactions, thereby revealing the presence of hidden planets (e.g., Lin & Papaloizou 1986; Crida et al. 2006; Pinilla et al. 2012). These hypothetical planets would need to possess considerable masses (≳10 M_⊕) in order to perturb the surrounding gas within the disk, as well as orbiting at large distances (> 10 AU), in line with the observed substructure locations (we refer to Fig. 1 by Lodato et al. 2019, and references therein). Moreover, some observations reveal the presence of gaps in disks younger than 1 Myr (e.g., Sheehan & Eisner 2018). If these gaps are indeed a result of the gravitational influence of embedded protoplanets, then it suggests that the formation of massive protoplanets occurs during the early stages of disk evolution.

There is already strong evidence supporting the notion that some substructures have formed thanks to the presence of protoplanets. In PDS 70, the giant planets PDS 70 b and c were detected within a cavity, orbiting at ~22 and ~35 AU, respectively (Keppler et al. 2018; Haffert et al. 2019). Another more recent evidence comes from the spirals excited by the planet AB Aur b at ~93 AU from its central star (Currie et al. 2022). These detections do not imply that all detected substructures must be due to the presence of a planet. However, it is necessary to study how common wide-orbit planet formation is and whether it can also statistically explain the origin of the substructures for which no direct evidence of protoplanets has been found. To date, the early formation mechanisms of distant and massive planets remain an open question.

In the classical picture, the core of a gas giant forms by accreting planetesimals from its vicinity (Wetherill & Stewart 1989; Kokubo & Ida 1996; Ormel et al. 2010), which is followed by the accretion of a gaseous envelope before disk dispersal (Safronov 1972; Pollack et al. 1996), typically within a few million years (Haisch et al. 2001; Soderblom et al. 2014; Tychoniec et al. 2020). Nevertheless, the formation of gas giant cores in the outer regions via planetesimally driven scenarios is hindered by long formation times (e.g., Thommes et al. 2003; Ida & Lin 2004; Bitsch et al. 2015; Johansen & Bitsch 2019; Lorek & Johansen 2022). Given this limitation concerning planetesimal accretion, a new paradigm known as pebble accretion was proposed to explain the mechanism that increases the growth rates of forming planets (Ormel & Klahr 2010; Lambrechts & Johansen 2012; Johansen & Lambrechts 2017). Pebble accretion allows for a faster growth because pebbles (~mm and cm sized particles) drift radially from the outside of the disk toward the center, continuously replenishing the accreting zone. Moreover, the gas in the vicinity of the protoplanet exerts a drag force that drains kinetic energy from the pebbles. However, such drag force would not be sufficient to increase the accreting rate of 100 km-sized planetesimals.

Although pebble accretion is a prospective mechanism for rapid core formation, the formation of wide-orbit planets still faces challenges such as limited mass reservoirs (Ormel 2017) and planetary migration (Ward 1997; Johansen et al. 2019). The latter occurs due to the gravitational interaction between the gas and the protoplanet, causing an inward migration to the central star. Indeed, that migration may be rapid enough to prevent the retention of gas giants in wide orbits (Coleman & Nelson 2014). It is therefore necessary to identify the characteristics that can lead to the formation of distant planets via pebble accretion.

Recent studies of wide-orbit planet formation embraced the idea of pressure bumps or rings. These not only prevent the planet from migrating too quickly, but they also accumulate enough solid build-up for growth (e.g., Morbidelli 2020; Chambers 2021; Jiang & Ormel 2023). However, as previously discussed, the substructures observed in disks might be elicited by other planets. In this work, we therefore focus on the formation of the earliest planets and consider disks with monotonic pressure profiles. We analyze the evolution of individual protoplanets located in the outer regions, with initial masses of M₀ = 0.01 M_⊕. The feeding material of these bodies is limited by the depletion of the pebble reservoir. Even though detailed numerical calculations have been implemented to estimate the evolution of solids (e.g., Brauer et al. 2008; Birnstiel et al. 2010, 2012; Stammler & Birnstiel 2022), due to the computational cost of these simulations, the joint study of the evolution of protoplanets while considering pebble flux decay is a difficult task. Some authors have dealt with this issue by employing small pebbles and thereby assuming a tightly coupled evolution of solids and gas (e.g., Liu et al. 2019; Johansen et al. 2019), while others have treated the flux as a free parameter (e.g., Bitsch et al. 2019; Lambrechts et al. 2019; Ogihara & Hori 2020). In this work, we present a new analytical expression that is an exact solution to the mass conservation equation for pebbles undergoing radial drift in a viscous gas disk. We validate this expression against the more complex computational analysis as (similarly) conducted by Appelgren et al. (2023). The new analytical model allows us to explore a new pathway for gas accretion, where the protoplanet can start to accrete once the pebble accretion rate has dropped substantially, in addition to the conventional pathway that requires the core to reach the pebble isolation mass (Lambrechts et al. 2014). The aim of this study is to present the new pebble flux model and to implement it to study the formation of distant cores and gas giants.

This work is structured as follows. In Sect. 2, we introduce the model of pebble accretion for the outer regions, where we include the derivation of the new analytical model for the pebble flux (in Sect. 2.3). In Sect. 3, we show the furthest possible core formation in different disk models, which we link to the gaps observed in protoplanetary disks. In Sect. 4, we include a simple gas accretion prescription and a new pathway for gas accretion. Furthermore, we highlight the importance of considering pebble depletion to explain gas accretion in distant orbits. We discuss the implications of our findings and the limitations in Sect. 5. Finally, we summarize our work in Sect. 6.

2 Pebble accretion in the outer regions

In this section, we present the model used to elucidate the evolution of a protoplanet in the outer regions. We first describe the disk structure and then proceed to the derivation of a new analytical pebble flux. The outcome of the derivation is specifically given in Eqs. (32) and (36) under two different assumptions for the pebble Stokes number. Consequently, we describe the growth rates via pebble accretion, followed by the migration of the growing body.

2.1 Disk structure

The gas surface density profile of a disk without substructures (Lynden-Bell & Pringle 1974) is described as: $${\Sigma }_{\text{g}}(r,t)=\frac{{\stackrel{.}{ℳ}}_{\text{g},0}}{3\pi v_1{\tilde{r}}^{\gamma }}{T}^{-\frac{5/2-\gamma }{2-\gamma }}\exp \left(\frac{-{\tilde{r}}^{(2-\gamma )}}{T}\right),$$(1)

where ${\stackrel{.}{ℳ}}_{\text{g},0}$ is the initial gas accretion rate onto the star in the inner regions, v₁ is the viscosity at the characteristic disk size R₁, γ is the viscous power-law index, $\tilde{r}\equiv r/R_1$ is the dimensionless position, and T is the dimensionless time defined through $$T\equiv \frac{t-t_0}{t_{\text{s}}}+1.$$(2)

Here, t₀ is the initial time of the disk when ${\stackrel{.}{ℳ}}_{\text{g}}={\stackrel{.}{ℳ}}_{\text{g},0}$, and t_s is the viscous timescale of the gas at a radial distance of R₁, which characterizes the time span required for the gas to undergo substantial radial transport. The viscous timescale is defined as: $$t_{\text{s}}\equiv \frac{1}{3{(2-\gamma )}^2}\frac{R_1^2}{v_1}.$$(3)

We describe the viscosity, v, using the α-disk model (Shakura & Sunyaev 1973), $$v=\alpha c_{\text{S}}H,$$(4)

where c_s is the sound speed and H is the gas scale-height. These two quantities are defined as: $$c_{\text{S}}=c_{\text{s},1}{\left(\frac{r}{\text{AU}}\right)}^{-\frac{\zeta }{2}}\text{,}$$(5) $$H=\frac{c_{\text{s}}}{\Omega }\text{,}$$(6)

where c_s,1 is the sound speed at 1 AU, ζ is the negative power-law index of the temperature, and $\Omega =\sqrt{G{M}_{\star }/r^3}$ is the Keplerian frequency. Equation (4) relates γ and ζ such that γ = 3/2 − ζ.

To set the fiducial values, we assume that for a solar-mass star at an initial time of t₀ = 0.2 Myr, there is an accretion rate of ${\stackrel{.}{ℳ}}_{\text{g},0}\approx {10}^{-7}{M}_{\odot }{\text{yr}}^{-1}$ and an accretion coefficient of α ~ 0.01 (Hartmann et al. 1998). We choose c_s,1 = 650 m s⁻¹ at 1 AU (Johansen et al. 2019). In the outer regions where viscous heating can be ignored, γ ≈ 15/14 and ζ ≈ 3/7 (Ida et al. 2016). Even though most observed protoplanetary disks appear to be small (e.g., Barenfeld et al. 2017; Tobin et al. 2020), we are interested in understanding substructure formation far from the central star, which calls for large disks. Hence, we assumed a fiducial initial

disk size of R₁ = 100 AU. Employing these values, we compute the total mass of the disk as: $$\begin{array}{c}{M}_{\text{g}}(t)={\displaystyle {\int }_0^{\infty }2}\pi r{\Sigma }_{\text{g}}(r,t)\text{d}r=\frac{2}{3}\frac{{\stackrel{.}{ℳ}}_{g,0}}{v_1}\frac{R_1^2}{(2-\gamma )}{T}^{-\frac{1}{2(2-\gamma )}}\\ \approx 0.17M_{\odot }{\left(\frac{\alpha }{0.01}\right)}^{-1}\left(\frac{{\stackrel{.}{ℳ}}_{\text{g},0}}{{10}^{-7}{M}_{\odot }{\text{yr}}^{-1}}\right)\\ \times {\left(\frac{c_{\text{s},1}}{650 \text{m} {\text{s}}^{-1}}\right)}^{-2}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{\frac{1}{2}}{\left(\frac{R_1}{100\text{AU}}\right)}^{2-\gamma }{T}^{-\frac{1}{2(2-\gamma )}}.\end{array}$$(7)

Hence, the initial disk-mass for our fiducial values is M_g,0 ≈ 0.17 M_⊙. For the same values, but while changing the disk size to a larger disk of R₁ = 300 AU, we instead obtain M_g,0 ≈ 0.48 M_⊙. These disks are stable since the Toomre parameter, Q, is higher than 1 everywhere in the disks (Toomre 1964). We list the fiducial values of the parameters in Table 1.

The radial and temporal dependence of the gas flux is: $${\stackrel{.}{ℳ}}_{\text{g}}(r,t)={\stackrel{.}{ℳ}}_{\text{g},0}{T}^{-\frac{5/2-\gamma }{2-\gamma }}\exp \left(-\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right)\times \left[1-2(2-\gamma )\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right].$$(8)

This flux changes direction at $R_{\text{t}}=R_1{\left[\frac{T}{2(2-\gamma )}\right]}^{1/(2-\gamma )}$. definition, the gas flux is also $${\stackrel{.}{ℳ}}_{\text{g}}\equiv -2\pi r{v}_{\text{r},\text{g}}{\Sigma }_{\text{g}},$$(9)

where v_r,g is the radial velocity of the gas. The negative sign is due to the definition of a positive flux as being directed toward the star. Knowing Σ_g and ${\stackrel{.}{ℳ}}_{\text{g}}$, the radial velocity of the gas is: $$v_{\text{r},\text{g}}=-\frac{{\stackrel{.}{ℳ}}_{\text{g}}}{2\pi r{\Sigma }_{\text{g}}}=-\frac{3}{2}\frac{v}{r}\times \left[1-2(2-\gamma )\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right].$$(10)

In Fig. 1, we show the gas surface density and the gas flux described by Eqs. (1) and (8), respectively. When the initial disk size is set to R₁ = 100 AU, the outward flux initially emerges at a radius greater than 50 AU and it then gradually moves outward, reaching 100 AU in 1 Myr. Considering the momentum redistribution of the gas is therefore clearly relevant for protoplanets forming at early stages and large distances.

Regarding the orbital motion, the gas is rotating at sub-Keplerian velocities, such that υ_ϕ,g ≈ υ_K − ∆υ, where ∆υ is: $$\Delta \upsilon =-\frac{1}{2}\left(\frac{H}{r}\right)\frac{\partial \ln P}{\partial \ln r}{c}_{\text{s}}=\frac{1}{2}\left(\frac{H}{r}\right)\chi c_{\text{s}}.$$(11)

We denote the negative logarithmic pressure gradient in the midplane as χ, $$\chi \equiv -\frac{\partial \ln P}{\partial \ln r}=\gamma +\frac{\zeta }{2}+\frac{3}{2}+(2-\gamma )\frac{{\tilde{r}}^{(2-\gamma )}}{T}={\chi }_0+(2-\gamma )\frac{{\tilde{r}}^{(2-\gamma )}}{T}.$$(12)

Here, χ₀ ≈ 2.79 is the pressure gradient in the inner regions of the protoplanetary disk. This pressure gradient is increased in the outer disk by the last term of Eq. (12).

Fig. 1 Gas structure in the outer regions of the protoplanetary disk at different times, for our fiducial model with an initial disk size of R1 = 100 AU and an initial mass of 0.17 M⊙. Top: gas surface density profile from Eq. (1). Bottom: radial gas flux from Eq. (8). The dots illustrate the location where the gas moves outwards with time.

Fig. 2 Limiting solid growth barriers across the protoplanetary disk for our fiducial model. Stokes number limited by fragmentation from Eq. (15), indicated by the green line, and limited by the radial drift from Eq. (16), indicated by the blue line (top). The red dashed-dotted line illustrates initial St when we set constant Stχ ≡ St ⋅ χ (see Sect. 2.3). Initial particle size for each growth barrier from Eq. (14), indicated by the dashed lines, and the growth timescale for reaching the corresponding size from Eq. (17), indicated by the solid lines (bottom). The dotted gray lines denote St ≈ 0.03 and τg ≈ 0.2 Myr at R1 = 100 AU set by radial drift.

2.2 Evolution of pebbles

Since the solid particles initially orbit at Keplerian speed and the gas at sub-Keplerian one, the gas exerts a drag force on the solids, which makes them drift radially in an efficient way when they are pebble-size bodies. The radial velocity of pebbles is described as (Weidenschilling 1977): $${\upsilon }_{\text{r},\text{p}}=\frac{{\upsilon }_{\text{r},\text{g}}}{1+{\text{St}}^2}-\frac{2\Delta \upsilon \text{St}}{1+{\text{St}}^2}.$$(13)

Here, St is the Stokes number, υ_r,g is the radial velocity of the gas from Eq. (10) and ∆υ is the Keplerian velocity reduction of

the gas from Eq. (11). The first term in Eq. (13) describes the advection mode of transport that occurs when the gas flux drags along solid particles, while the second term corresponds to the radial drift toward higher pressure. In the Epstein regime, St is defined as (Weidenschilling 1977; Drążkowska et al. 2023): $$\text{St}=\frac{\pi }{2}\frac{{\text{a}}_{\text{s}}{\rho }_{\text{s}}}{{\Sigma }_{\text{g}}},$$(14)

where a_s is the particle size and ρ_s is its internal density. We chose ρ_s = 1 g cm⁻³ as a nominal value.

Regarding the growth of solids, the µm-sized primordial dust particles will easily stick together to form larger particles. In order to determine the particle size, we consider for simplicity that solid growth can be limited either by fragmentation or radial drift (for other growth barriers, see review by Testi et al. 2014). To describe how pairwise collisions can result in fragmentation, the fragmentation velocity υ_f is employed. We adhere to υ_f ~ 1 m s⁻¹, as reported by laboratory experiments done by Güttler et al. (2010) for silicate grains. The maximum Stokes number is then described by (Ormel & Cuzzi 2007; Birnstiel et al. 2009): $${\text{St}}_{\text{frag}}\approx \frac{1}{3{\alpha }_{\text{t}}}{\left(\frac{{\upsilon }_{\text{f}}}{c_{\text{s}}}\right)}^2,$$(15)

where α_t is the midplane turbulence. Since weak gas turbulence has been inferred from dust observations in the outer regions, we chose α_t = 10⁻⁴ as a fiducial value¹ (see review by Pinte et al. 2023, and the references therein). The Stokes number of particles limited by the radial drift is (Ida et al. 2016): $${\text{St}}_{\text{drift}}\approx \frac{\sqrt{3\pi }}{80}\frac{v_{\text{K}}}{\Delta v}{Z}_0.$$(16)

The maximum pebble-size at a certain location is the minimum between the fragmentation and drift limit. In the top panel of Fig. 2, we show the fragmentation and drift limits for fiducial values across the protoplanetary disk and we note that at R₁ = 100 AU, growth is limited by radial drift with St = 0.03.

The timescale required to grow from 1 µm to the particle size a_s is (Takeuchi & Lin 2005; Brauer et al. 2008; Sato et al. 2016): $${\tau }_{\text{g}}\approx \frac{4}{\sqrt{3\pi }}\frac{1}{Z_0\Omega }\ln \left(\frac{a_{\text{s}}}{\mu \text{m}}\right).$$(17)

Here Z₀ is the initial metallicity of the disk and ω is the Keplerian frequency. In the bottom panel of Fig. 2, we show that the dust growth timescale is approximately 0.2 Myr at R₁ = 100 AU. Therefore, already at t ≈ 0.2 Myr pebbles have formed and they drift from the outer to the inner regions of the disk, thus a pebble flux, ${\stackrel{.}{ℳ}}_{\text{p}}$, can be defined to describe the carried mass (Lambrechts & Johansen 2014). The expression of the pebble flux is given by: $${\stackrel{.}{ℳ}}_{\text{p}}\equiv -2\pi r{\upsilon }_{\text{r},\text{p}}{\Sigma }_{\text{p}},$$(18)

where υ_r,p is the radial velocity of pebbles from Eq. (13) and Σ_P is the pebble surface density, which we calculate below.

2.3 Derivation of the analytical pebble flux

Since the dust mass distribution is dominated by the mass of the largest particles, especially in the drift-limited regime (Birnstiel et al. 2012), we assume that the pebble-to-gas ratio is equal to the total metallicity. The pebble-to-gas ratio is defined as: $$Z=\frac{{\Sigma }_{\text{p}}}{{\Sigma }_{\text{g}}}.$$(19)

The pebble flux can be written in terms of the metallicity by dividing Eqs. (9) and (18), $${\stackrel{.}{ℳ}}_{\text{p}}=Z\frac{{\upsilon }_{\text{r},\text{p}}}{{\upsilon }_{\text{r},\text{g}}}{\stackrel{.}{ℳ}}_{\text{g}}.$$(20)

For simplicity, we define now an auxiliary quantity $h=\frac{{\upsilon }_{\text{r},\text{p}}}{{\upsilon }_{\text{r},\text{g}}}{\stackrel{.}{ℳ}}_{\text{g}}$ so that ${\stackrel{.}{ℳ}}_{\text{p}}=Zh$. From Eqs. (8), (10), and (13), we get the ratio of particle speed to gas speed and, thus, h(r, t) as: $$\begin{array}{l}\frac{{\upsilon }_{\text{r},\text{p}}}{{\upsilon }_{\text{r},\text{g}}}=\frac{1}{1+{\text{St}}^2}\left(1-\frac{2\Delta \upsilon \text{St}}{{\upsilon }_{\text{r},\text{g}}}\right)\hfill \\ =\frac{1}{1+{\text{St}}^2}\left[1+\frac{2}{3}\frac{\chi \text{St}}{\alpha }\frac{1}{1-2(2-\gamma )\frac{{\tilde{r}}^{(2-\gamma )}}{T}}\right],\hfill \end{array}$$(21) $$\begin{array}{c}h(r,t)=\frac{v_{\text{r},\text{p}}}{v_{\text{r},\text{g}}}{\stackrel{.}{ℳ}}_{\text{g}}\\ =\frac{1}{1+{\text{St}}^2}\left[1-2(2-\gamma )\frac{{\tilde{r}}^{(2-\gamma )}}{T}+\frac{2}{3}\frac{\chi \text{St}}{\alpha }\right]\\ \times {\stackrel{.}{ℳ}}_{\text{g},0}{T}^{-\frac{5/2-\gamma }{2-\gamma }}\exp \left(-\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right).\end{array}$$(22)

Here, $\tilde{r}=r/R_1$ and T (see Eq. (2)) are the dimensionless spatial and time variables. As shown in Eq. (12), χ depends on both r and t, and generally St can also depend on r and t. Defining then: $$b(r,t)=\frac{2}{3}\frac{\chi (r,t)\mathrm{St}(r,t)}{\alpha },$$(23)

we can rewrite h as: $$\begin{array}{c}h(r,t)=\frac{1}{1+{\text{St}}^2}\left[1+b(r,t)-2(2-\gamma )\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right]\\ \times {\stackrel{.}{ℳ}}_{\text{g},0}{T}^{-\frac{5/2-\gamma }{2-\gamma }}\exp \left(-\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right).\end{array}$$(24)

Since h(r, t) is a known function of the underlying evolution of the α-disk and the speed of the pebbles, our goal now is to derive Z(r, t) so that we can get the analytical form of ${\stackrel{.}{ℳ}}_{\text{p}}$ in Eq. (20). Neglecting the diffusivity of solid particles within the gas, the pebble flux fulfills its own continuity equation as: $$r\frac{\partial {\Sigma }_{\text{p}}}{\partial t}-\frac{1}{2\pi }\frac{\partial {\stackrel{.}{ℳ}}_{\text{p}}}{\partial r}=0.$$(25)

We can rewrite this equation in terms of Z by replacing Σ_P> = ZΣ_g and ${\stackrel{.}{ℳ}}_{\text{p}}=Zh$ from Eqs. (19) and (20). We make the change of variables from r and t to the dimensionless parameters $\tilde{r}$ and T, respectively, to simplify the equation. Applying the chain rule, we have: $$\frac{\partial }{\partial t}\equiv \frac{1}{t_{\text{s}}}\frac{\partial }{\partial T}\equiv 3{(2-\gamma )}^2\frac{v_1}{R_1^2}\frac{\partial }{\partial T},\frac{\partial }{\partial r}\equiv \frac{1}{R_1}\frac{\partial }{\partial \tilde{r}}.$$(26)

The continuity equation for the pebbles, therefore, is $$3{(2-\gamma )}^2\frac{v_1}{R_1}\tilde{r}\frac{\partial \left(Z\cdot {\Sigma }_{\text{g}}\right)}{\partial T}-\frac{1}{2\pi }\frac{1}{R_1}\frac{\partial (Z\cdot h)}{\partial \tilde{r}}=0.$$(27)

We simplify the equation by multiplying away R₁ and expanding the partial derivatives to obtain: $$3{(2-\gamma )}^2{v}_1\tilde{r}\left(Z\frac{\partial {\Sigma }_{\text{g}}}{\partial T}+{\Sigma }_{\text{g}}\frac{\partial Z}{\partial T}\right)-\frac{1}{2\pi }\left(Z\frac{\partial h}{\partial \tilde{r}}+h\frac{\partial Z}{\partial \tilde{r}}\right)=0.$$(28)

To solve the partial differential equation, we substitute Σ_g and h from Eqs. (1) and (24) and compute their (dimensionless) temporal and spatial derivatives respectively. The time derivative of Σ_g $$\frac{\partial {\Sigma }_{\text{g}}}{\partial T}=\frac{{\stackrel{.}{ℳ}}_{\text{g},0}}{3\pi v_1{\tilde{r}}^{\gamma }}{T}^{-\frac{5/2-\gamma }{2-\gamma }}\exp \left(-\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right)\frac{1}{T}\left(-\frac{5/2-\gamma }{2-\gamma }+\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right).$$(29)

Computing $\frac{\partial h}{\partial \tilde{r}}$ is not trivial and, therefore, we must make an approximation. The largest pebbles in the outer regions do not usually exceed St ~ 0.1. Therefore, for simplicity, $\frac{1}{{\text{St}}^2+1}\approx 1$. Consequently, the spatial derivative of h is simplified as: $$\begin{array}{c}\frac{\partial h}{\partial \tilde{r}}=\left\{\frac{\partial b}{\partial \tilde{r}}-(2-\gamma )\frac{{\tilde{r}}^{(1-\gamma )}}{T}\left[1+b+2(2-\gamma )\left(1-\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right)\right]\right\}\\ \times {\stackrel{.}{ℳ}}_{\text{g},0}{T}^{-\frac{5/2-\gamma }{2-\gamma }}\exp \left(-\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right).\end{array}$$(30)

By replacing the calculated partial derivatives in Eq. (28) and simplifying terms, the general evolution equation for $Z(\tilde{r},T)$ becomes: $$\begin{array}{c}\frac{1}{2}\left[1+b-2(2-\gamma )\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right]\frac{\partial Z}{\partial \tilde{r}}-{(2-\gamma )}^2{\tilde{r}}^{(1-\gamma )}\frac{\partial Z}{\partial T}\\ =Z\frac{{\tilde{r}}^{(1-\gamma )}}{T}\left[\frac{2-\gamma }{2}b-\frac{T}{2{\tilde{r}}^{(1-\gamma )}}\frac{\partial b}{\partial \tilde{r}}\right].\end{array}$$(31)

To compute the solution of this partial differential equation (PDE), first we need to obtain the analytical expression of b(r, t) from Eq. (23). Since the value of St likely changes only weakly with distance in the outer regions (see top panel of Fig. 2), we first assume that St is constant. Assuming that the initial pebble-to-gas ratio remains equal to the initial dust-to-gas ratio Z₀ across the disk, that is, $Z(\tilde{r},1)=Z_0$, the expression of Z as solution of Eq. (31) is: $$\begin{array}{c}Z(\tilde{r},T)=Z_0{T}^{\frac{1}{2(2-\gamma )}+\frac{\text{St}}{3\alpha }}\\ \times \exp \left\{-\left[\frac{1}{2(2-\gamma )}\left(2{\chi }_0+\frac{3\alpha }{\text{St}}\right)+\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right]\left[T^{\frac{\text{St}}{3\alpha }}-1\right]\right\}.\end{array}$$(32)

For the derivation of this result, we refer to Appendix A. Another possible analytical solution for Eq. (31) comes from assuming that St ⋅ χ is constant, as consequently, b = b₀ is constant (see Eq. (23)). The governing PDE can be simplified to: $$\begin{array}{c}\frac{1}{2}\left[1+b_0-2(2-\gamma )\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right]\frac{\partial Z}{\partial \tilde{r}}-{(2-\gamma )}^2{\tilde{r}}^{(1-\gamma )}\frac{\partial Z}{\partial T}\\ =Z\frac{{\tilde{r}}^{(1-\gamma )}}{T}\frac{2-\gamma }{2}{b}_0.\end{array}$$(33)

We introduce an ad hoc assumption to simplify the equation by setting $\frac{\partial Z}{\partial \tilde{r}}=0$. Consequently, the PDE reduces to an ordinary differential equation (ODE) since the remaining dependence on $\tilde{r}$ cancels out, as follows: $$\begin{array}{c}-{(2-\gamma )}^2{\tilde{r}}^{(1-\gamma )}\frac{\partial Z}{\partial T}=Z\frac{{\tilde{r}}^{(1-\gamma )}}{T}\frac{2-\gamma }{2}{b}_0\Rightarrow \\ -(2-\gamma )\frac{\partial Z}{\partial T}=\frac{Z}{T}\frac{b_0}{2}.\end{array}$$(34)

Then, we integrate the ODE for the initial conditions $Z(\tilde{r},1)=Z_0$, $$Z(T)=Z_0{T}^{-\frac{b_0}{2(2-\gamma )}}.$$(35)

The solution is consistent with the assumption $\frac{\partial Z}{\partial \tilde{r}}=0$ and, since the solution of the PDE, which is subject to the specified initial conditions, is unique, Eq. (35) is, in fact, the solution of the PDE from Eq. (33). We also include the derivation of Eq. (33) with the method of characteristics from Appendix B to show that we get the same result. Since $b_0=\frac{2}{3}\frac{\chi \cdot \text{St}}{\alpha }$, the solution becomes: $$Z(T)=Z_0{T}^{-\frac{1}{2-\gamma }\frac{\chi \cdot \text{ St }}{3\alpha }},$$(36)

where χ ⋅ St is constant and is abbreviated as "constant St_χ" below. Later in this paper, we explore how this assumption compares to a model that includes fragmentation limited particle sizes. In order to implement constant St_χ we first set constant $b_0=\frac{2}{3}\frac{\chi \left(R_1,t_0\right)\cdot \mathrm{St}\left(R_1,t_0\right)}{\alpha }$. Then, the variation of St in both space and time is described by $\mathrm{St}(r,t)=\frac{3}{2}\frac{b_0\alpha }{\chi (r,t)}$ with χ defined in Eq. (12) (see top panel Fig. 2). In Appendix C, we demonstrate an alternative derivation of the metallicity, not involving the continuity equation, that yields the same result as in Eq. (36).

To summarize, for constant values of St and St_χ, we can derive the evolution of the metallicity via Eqs.(32) and (36), respectively². In Fig. 3, we show that the two expressions are similar until the metallicity drops to ~10% of its original value. After that, the constant St case displays a faster decrease than constant St_χ.

The analytical pebble flux is computed by replacing the metallicity in Eq. (20). Even though the full analytical expression of ${\dot{ℳ}}_{\text{p}}$ is nontrivial, we note from Eq. (21) that at r = 0, the pebble-to-gas flux ratio is described by the simple relation ${\stackrel{.}{ℳ}}_{\text{p}}/{\stackrel{.}{ℳ}}_{\text{g}}=(1+b)Z$. In Fig. 4, we compare the full analytical pebble flux expression with a numerical simulation as conducted by Appelgren et al. (2023), excluding disk formation and photo-evaporation, and employing the same disk temperature profile as in this work. We also include the case where Z= Z₀ does not evolve (for a detailed comparison with alternative analytical approaches from the literature, see Appendix E). Overall, we see that the model with constant Z = Z₀ strongly overestimates the pebble flux after a few hundred thousand years. Otherwise, both of the new analytical expressions imitate the behavior of a more complex computer simulation. After 1 Myr, the constant St_χ case nevertheless describes the pebble flux decay more accurately. In Fig. 5, we compute the cumulative mass of pebbles at two different locations and, overall, the derived expressions properly estimate the crossing mass.

Given that both constant St and constant St_χ assumptions provide an approximated crossing mass value of the flux in the outer regions (and given that a constant St_χ properly describes the decay of the flux), we adhere to the latter to model the planetary growth.

Fig. 3 Evolution of the analytical expression for metallicity Z according to constant St (Eq. (32)), indicated by dashed lines, and constant Stχ (Eq. (36)), indicated by solid lines. For both models, we apply the initial condition Z0 = 0.01 at 0.2 Myr and therefore the initial lines overlap. Then, until the metallicity drops to ~10% at approximately 1 Myr, both expressions yield a similar evolution. However, the metallicity decreases faster for the constant St case in the following Myrs. The employed fiducial values are listed in Table 1.

Fig. 4 Comparison between the analytical model (solid lines) with the numerical simulations from Appelgren et al. (2023; dashed lines) at different times. We assume Z0 = 0.008 and ${\dot{ℳ}}_{\text{g},0}=6\times {10}^{-8}{M}_{\odot }{\text{yr}}^{-1}$ to match the simulation, and St= 0.03 in the three analytical expressions. Left: assuming constant metallicity Z(r, t) = Z0 significantly overestimates the pebble flux over the lifespan of the disk. Center: assuming constant Stχ (see Eq. (32)) gives a relatively good match to the overall behavior of the pebble flux. Right: assuming constant St gives an accurate description over the first ~1 Myr, but it severely underestimates the flux at late times, since the value of St is limited by the growth timescale (see Fig. 2) in the outer disk in the simulations.

Fig. 5 Cumulative mass of drifting pebbles crossing 50 AU and 100 AU according to different analytical models and according to the numerical simulation from Appelgren et al. (2023). In the analytical models shown here, we assume that pebbles have grown to completion and started drifting at t0 = 0.2 Myr and that the initial metallicity is Z0 = 0.008. Assuming constant Z over time (blue line) significantly overestimates the crossing mass compared to the simulation (black line). In contrast, the cumulative masses from the new analytical models (red and yellow lines) approach the simulated case. We note that at 50 AU, the cumulative masses from these analytical models are slightly underestimated when comparing to the simulated case since Z0(50 AU) > 0.008 in the simulation.

2.4 Growth via pebble accretion

The planetesimals in the outer regions of the protoplanetary disk most likely form by the streaming instability (SI; Johansen et al. 2014). Lyra et al. (2023) found that growth via pebble accretion is possible directly after the planetesimal formation by SI. Thus, we examine here the evolution of a typical initial protoplanet mass of ~ 0.01 M_⊕ that could form directly by SI at large distances (Liu et al. 2020).

To formulate the growth rate via pebble accretion, we first consider whether the planetary mass is sufficient to accrete from the complete vertical extent of the pebble layer (2D accretion) or not (3D accretion). Secondly, we consider if the relative velocity between the protoplanet and pebbles is dominated by the sub-Keplerian gas flow (Bondi Regime) or by the Keplerian shear (Hill regime).

Regarding the first aspect, the 2D regime is relevant when the pebble accretion radius, R_acc, is larger than the pebble scale-height (Dubrulle et al. 1995; Johansen et al. 2014): $$H_{\text{p}}=H\sqrt{\frac{{\alpha }_{\text{t}}}{{\alpha }_{\text{t}}+\text{St}}},$$(37)

where α_t is the aforementioned midplane turbulence. The characteristics favoring the 2D accretion scenario are a large planetary mass and the settling of pebbles. In contrast, in the 3D regime, the protoplanet only has access to a fraction of the pebble layer. Therefore, the efficiency in the 3D regime is lower than in the 2D regime. From Johansen & Lambrechts (2017), the expressions for the growth rate of the protoplanet in each regime are: $${\dot{M}}_{2\text{D}}=2R_{\text{acc}}{\Sigma }_{\text{p}}\delta \upsilon ,$$(38) $${\dot{M}}_{3\text{D}}=\pi R_{\text{acc}}^2{\rho }_{\text{p}}\delta \upsilon ,$$(39)

where Σ_P is the pebble surface density that can be analytically derived from Eqs. (19) and (36), $${\Sigma }_{\text{p}}(r,t)=\frac{{\stackrel{.}{ℳ}}_{\text{g},0}}{3\pi v_1{\tilde{r}}^{\gamma }}{Z}_0{T}^{-\frac{1}{2-\gamma }\left(\frac{5}{2}-\gamma +\frac{\chi \text{st}}{3\alpha }\right)}\exp \left(\frac{-{\tilde{r}}^{(2-\gamma )}}{T}\right),$$(40)

assuming that all dust grows into pebbles. In Eq. (39), ${\rho }_{\text{p}}=\frac{1}{\sqrt{2\pi }}\frac{{\Sigma }_{\text{p}}}{H_{\text{p}}}$ is the pebble density in the midplane. Both rates depend on δυ, the approach velocity between the pebbles and the protoplanet, defined as: $$\delta \upsilon \equiv \Omega R_{\text{acc}}+\Delta \upsilon ,$$(41)

where ∆υ is the sub-Keplerian velocity reduction of the gas from Eq. (11). For the transition between 3D and 2D to be continuous, ${\dot{M}}_{2\text{D}}={\dot{M}}_{3\text{D}}$ must hold at a certain accretion stage. From this equality, we get that the transition occurs when: $$\frac{R_{\text{acc}}}{H_{\text{p}}}=\sqrt{\frac{8}{\pi }}\approx \mathrm{1.6.}$$(42)

The analytical form of the accretion radius R_acc will depend on whether the protoplanet is accreting in the Bondi or Hill regime. The Hill and Bondi radius of a protoplanet are defined as: $$R_{\text{H}}\equiv r{\left(\frac{M}{3M_{\star }}\right)}^{1/3},$$(43) $$R_{\text{B}}\equiv \frac{GM}{\Delta {\upsilon }^2},$$(44)

where r and M are the position and mass of the protoplanet. Johansen & Lambrechts (2017) derived an expression for the effective accretion radius in each regime, as follows: $$R_{\text{acc}}={\left(\frac{\text{St}}{0.1}\right)}^{1/3}{R}_{\text{H}}\text{\hspace{1em} (Hill regime),}$$(45) $$R_{\text{acc}}={\left(\frac{4\text{St}\Delta \text{v}}{\Omega R_{\text{B}}}\right)}^{1/2}{R}_{\text{B}}\text{\hspace{1em} (Bondi regime), }$$(46)

and by equating the two accretion radii, we can determine the transitional mass between the two regimes as: $$M_{\text{t}}=\frac{25}{144}\frac{\Delta {\upsilon }^3}{G\Omega }\frac{1}{\text{St}}.$$(47)

The accretion will occur in the Hill regime if M ≥ M_t and in the Bondi regime if M < M_t. In Fig. 6, we show the accretion rate calculated via Eqs. (38) or (39), depending on Eq. (42), and by replacing Eqs. (45) or (46), depending on Eq. (47). This approach ensures a smooth and continuous growth rate during the transition from one regime to another.

The maximum core mass that the protoplanet can attain is known as the pebble isolation mass, M_iso (Lambrechts et al. 2014). According to the 3D simulations in Bitsch et al. (2018), if the protoplanet grows massive enough to carve out a gap of a depth of ~10–20%, the pressure gradient in the outer edge is reversed and pebbles are pushed outward, causing pebble accretion to cease and enabling gas accretion to occur. The authors derived a scaling law of M_iso for protoplanets orbiting solar-mass stars, such that: $$\begin{array}{c}{M}_{\text{iso}}(r)=25M_{\oplus }{\left(\frac{H/r}{0.05}\right)}^3\times \left[0.34{\left(\frac{\log {10}^{-3}}{\log {\alpha }_{\text{t}}}\right)}^4+0.66\right]\\ \times \left[1-\frac{-{\chi }_0+2.5}{6}\right],\end{array}$$(48)

where α_t is the midplane turbulence and χ₀ s the negative logarithmic pressure gradient in Eq. (12). The pebble isolation mass increases for small pebbles (Bitsch et al. 2018), but we neglected this effect due to our choice of relatively large pebbles.

Fig. 6 Initial pebble accretion rate onto a protoplanet as a function of its position and mass. The separation between 3D and 2D accretion regimes and between Bondi and Hill regimes are indicated. In the outer regions, we need to consider 3D and 2D as well as Bondi and Hill regimes. Employed fiducial values are listed in Table 1.

2.5 Migration

Protoplanets undergo inward radial migration while they grow. Planets that are not massive enough to open a gap in the gas fall into the type-I migration regime. We describe the migration speed of these protoplanets by the standard scaling law derived in Tanaka et al. (2002), $$\dot{r}=-k_{\text{mig}}\frac{M}{M_{\star }}\frac{{\Sigma }_{\text{g}}{r}^2}{M_{\star }}{\left(\frac{H}{r}\right)}^{-2}{\upsilon }_{\text{K}},$$(49)

where υ_K is the Keplerian velocity and k_mig the constant prefactor that was fitted using 3D numerical simulations in D'Angelo & Lubow (2010) $$k_{\text{mig }}=-2(1.36+0.62\gamma +0.43\zeta ).$$(50)

Here, γ and ζ are the already discussed power-law indexes of the viscosity and the midplane temperature or sound speed, respectively.

When a planet reaches a certain mass threshold, M_gap, it creates a density gap along its orbit that leads to reduced migration rates. If the gas does not flow through this gap, the planet is forced to migrate at the same speed as the viscous accretion of the gas (Lin & Papaloizou 1986). This type of migration is referred to as type II migration. However, hydrodynamical simulations have demonstrated that gas will cross the gap (e.g., Dürmann & Kley 2015). In Kanagawa et al. (2018), a new physical model was proposed in which the torque exerted by the gas that crosses the gap depends on the surface density at the bottom of the gap, Σ_gap rather than the unperturbed density, Σ_g as in the type-I regime. Then, Σ_gap decreases as the protoplanet’s mass increases, thereby slowing down migration. The migration effectively slows down when the gap depth (Σ_gap/Σ_g) is reduced to approximately 50%. Consequently, the protoplanet reaches M_iso slightly before M_gap. In Johansen et al. (2019), it was suggested that a relative gap height of around 85% is sufficient to reach M_iso. They found that M_gap ≈ 2.3 M_iso and provided the modified migration equation as follows: $$\dot{r}=\frac{{\Sigma }_{\text{gap }}}{{\Sigma }_{\text{g}}}\cdot {\dot{r}}_{\text{I}}=\frac{{\dot{r}}_{\text{I}}}{1+{\left[M/\left(2.3M_{\text{iso}}\right)\right]}^2},$$(51)

where ${\dot{r}}_{\text{I}}$ represents the classical type-I migration rate from Eq. (49) and M_iso is described in Eq. (48).

3 Formation of wide-orbit cores

In this section, we analyze the evolution of protoplanets that grow via pebble accretion while migrating until they reach the pebble isolation mass, M_iso. We also included gas accretion (see Sect. 4).

3.1 Growth tracks with the new pebble flux model

We applied the new analytical models derived in Sect. 2.3 for calculating the growth tracks of protoplanets and compared them with the results obtained by the assumption of constant metallicity (or coupling between solids and gas) over time.

For that purpose, we initialized the protoplanets of M₀ = 0.01 M_⊕ formed at t_0,p = 0.2, 0.5 or 0.8 Myr and located at r₀ = 20, 50 or 80 AU. We placed them within a disk for the three different descriptions of the pebble flux (or metallicity) and we varied the pebble Stokes number as St = 0.01, 0.03, or 0.06. For the rest of the disk parameters, we adhered to the values motivated in Sect. 2 and listed in Table 1. We calculated the planetary growth and migration as described in Sects. 2.4 and 2.5, respectively, and stopped the calculations either at the end of the disk lifetime (we set t_f = 5 Myr) or earlier if they have reached the value of M_iso from Eq. (48).

We show the results of these calculations in Fig. 7. For St = 0.01, the growth tracks of protoplanets started at t_0,p = 0.2 Myr are similar in the three models. However, for protoplanets formed at t_0,p = 0.5 and 0.8 Myr, the growth is overestimated when using the model with a constant Z. For St = 0.03, the growth in the outer regions is overestimated even for protoplanets formed at t_0,p = 0.2 Myr. In contrast, models with constant St_χ and St have similar outcomes. Finally, for St = 0.06, the Z constant model only gives proper results when t_0,p = 0.2 Myr and r₀ = 20 AU.

Overall, Fig. 7 demonstrates that it is necessary to consider the decay of the pebble flux to prevent overestimating the growth. For the subsequent calculations, we therefore adhere to the Z(t) model with the constant St_χ from Eq. (36).

Fig. 7 Numerically integrated growth tracks for three different pebble fluxes or metallicity models (row-to-row) and different Stokes numbers (column-to-column). In the constant Stχ model, St = 0.01, 0.03 and 0.06 are employed to set the initial St at radial distance R1 (see top panel Fig. 2). Protoplanets with initial mass of M0 = 0.01 M⊕ are placed at r0 = 20, 50 or 80 AU at different formation times (values indicated in the legend). We stop the integration either at the end of the disk lifetime of tf = 5 Myr or earlier if they reach Miso, indicated by the gray dashed-dotted line. When Z is constant, which is equivalent to assuming that the pebble-to-gas flux ratio in the inner regions is constant, the growth tracks are independent of t0,p. Consequently, the growth is overestimated compared to constant St and Stχ, especially when St = 0.03 or 0.06 and when t0,p = 0.5 Myr or 0.8 Myr. Models with constant Stχ and St from Eqs. (36) and (32), respectively, yield similar results.

3.2 Location of the furthest cores in different scenarios

Given that the protoplanet reaches M_iso by reversing the pressure gradient when carving out a gap in its vicinity, the formation of the core is related to the formation of gaps. Hence, if we can demonstrate that the formation of distant cores is possible, we can elucidate a plausible origin for the observed gaps far away from their central star. For this purpose, we investigated how far out cores of giant planets can form by reaching M_iso in different disks. We then compared the location of the furthest cores with the observed gap locations in protoplanetary disks.

Here, we additionally varied the Stokes number (again 0.01, 0.03, and 0.06), the initial metallicity of Z₀ = 0.01 and 0.02, and the midplane turbulence of α_t = 10⁻⁴ and 10⁻⁵ (see Table 2). Placing protoplanets of 0.01 M_⊕ at any location in the disk (r₀ ≤ R₁), we find the location of the furthest cores for different formation times: t_0,p, = 0.2, 0.3, 0.5 and 0.7 Myr. We performed the calculations for the disk sizes R₁ = 100 and 300 AU. We note that when R₁ = 300 AU, we obtain St_dritft ≈ 0.02 and t₀ ≈ 0.7 Myr for our fiducial values (see Fig. 2). For simplicity, we only varied a single parameter at a time and for the remaining disk parameters, we adhered to the values motivated in Sect. 2. However, we address this assumption in the discussion.

In Fig. 8, we show the location of furthest cores for different scenarios (different bars) for different formation times of the protoplanet (different subfigures) and for different disk sizes (different color bars). We also include the location of the observed gaps in protoplanetary disks (horizontal green lines). First, we see that the location of the furthest core depends on the disk properties, including the disk size. Closer to the initial time, t₀, the pebble flux is higher and therefore the growth is faster when the protoplanet forms early on. However, the level of impact of varying the formation time will depend on the disk parameters. On the one hand, when t_0,p = 0.2 Myr, the furthest cores form when St = 0.06 (abbreviation "hst"). In these scenarios where there is a strong but short-lasting pebble flux, the furthest distance drops quickly when varying t_0,p (especially for the smaller disk); for example, for the case of St = 0.06, Z₀ = 0.02 (hst.hz) and disk size R₁ = 100 AU, initially the furthest core can form beyond 50 AU, but if the protoplanets form 0.5 Myr later, the furthest core do not reach 10 AU. On the contrary, in scenarios where St = 0.01, the furthest cores are initially less distant but do not depend so strongly on the formation time due to the weak but longer-lasting pebble flux. Therefore, when varying St, planet formation faces the dilemma that increasing the pebble flux will accelerate the decay and consequently requires protoplanets to start their growth earlier.

We also see that a higher Z₀ leads to the formation of more distant cores. Doubling Z₀ doubles the growth rate of protoplanets (see Eqs. (38)–(40)) but does not affect the decay time of the flux, nor the migration rate. Therefore, increasing metallicity has a positive outcome on the formation of distant cores; for instance, when t_0,p = 0.2 Myr, in the fiducial simulation ("fid" and "hz") the most distant core moves from ~30 to ~50 AU when R₁ = 100 AU, and from ~70 to ~120 AU when R₁ = 300 AU.

Decreasing the turbulence, α_t, increases the location of the most distant core. As the pebble scale height, H_p, is smaller for lower turbulence, the 3D accretion rate is higher and the transition from 3D to 2D accretion occurs earlier (see Eq. (42)).

Decreasing α_t has a higher impact for the case where St = 0.01; since the pebble scale-height is larger and the pebble flux is weaker, the 3D accretion stage is longer.

Fig. 8 Location of the furthest cores formed by reaching Miso, marked by the height of the colored bar, for different scenarios (specified in Table 2) and for protoplanets initialized at different formation times, t0,p, in each panel. Horizontal green lines in the background illustrate the location of the observed gaps in protoplanetary disks that could be caused by planetary cores placed in the same location (Data from http://ppvii.org/chapter/12/figure7.txt/). Whent0,p = 0.2 Myr (top-left panel), favorable disk parameters for forming distant cores up to ~ R1/2 are high St, high Z0 and low αt. “hst” runs with St = 0.06 also show that when initial protoplanets are injected late (t0,p > 0.5 Myr, bottom panels) in disks with R1 = 100 AU, the pebble flux by that point is too weak to grow wide-orbit cores. We note that for R1 = 300 AU, the t0,p = 0.7 Myr case is more consistent with the pebble growth timescale in outer disk (see Fig. 2). Overall, many observed gaps in protoplanetary disks could be caused by planetary cores, but the most distant gaps could only form in disks with high initial metallicity and a large disk size (see Sect. 5.2 for discussion).

Fig. 9 Location and formation time of the furthest cores simulated in different scenarios from Table 2 when protoplanets emerge at t0,p = 0.2 Myr (left). Scenarios with the same St are highlighted by their colored groups. Final location of the furthest cores after gas accretion (Sect. 4.1) and type II migration (Sect. 2.5) at the end of the disk lifetime of 5 Myr (right). An example growth track is shown at the bottom-right to facilitate the readability of the axes. In all scenarios, protoplanets migrate several AU during gas accretion. This phenomenon is visually depicted as a significant horizontal displacement from the diagonal line rcore = rf. The displacement is more pronounced in disks with high St (purple and blue groups) due to an earlier formation of the cores.

4 Formation of wide-orbit gas giants

In this section, we include a simple model for gas accretion and study the planetary evolution after reaching M_iso. In addition, we implement an alternative pathway for gas accretion that a planet can enter before reaching the pebble isolation mass but after the pebble accretion rate has decayed substantially.

4.1 Simple prescription for gas accretion

We assume that gas accretion starts with the contraction of the gaseous envelope at a rate suggested by Ikoma et al. (2000), $${\dot{M}}_{\text{KH}}={10}^{-5}{M}_{\oplus }{\text{yr}}^{-1}{\left(\frac{M}{10M_{\oplus }}\right)}^4{\left(\frac{\kappa }{0.1 {\text{m}}^2 {\text{kg}}^{-1}}\right)}^{-1},$$(52)

where κ is the opacity of the envelope. We take κ = 0.005 m² kg⁻¹ as in Johansen et al. (2019). Since the contraction accelerates at higher mass, the planet might eventually become so massive that growth will be restricted by the supply of gas flowing into its Hill sphere. Once this occurs, the growth rate of the protoplanet will be equal to the rate of gas supply, as described by Tanigawa & Tanaka (2016) and Ida et al. (2018): $${\dot{M}}_{\text{disk }}=0.29{\left(\frac{H}{r}\right)}^{-2}{\left(\frac{M}{M_{\star }}\right)}^{4/3}{\Sigma }_{\text{g}}{r}^2\Omega \frac{{\Sigma }_{\text{gap }}}{{\Sigma }_{\text{g}}}\text{, }$$(53)

where Σ_gap/Σ_g is the same as in Sect. 2.5.

The growth rate cannot be larger than the global accretion rate of the gas flux within the disk ${\stackrel{.}{ℳ}}_{\text{g}}$ from Eq. (8). Indeed, Lubow & D’Angelo (2006) estimated that the maximum accretion rate onto the protoplanet is approximately 80% of the gas flux and therefore: $$\dot{M}=\min \left[{\dot{M}}_{\text{KH}},{\dot{M}}_{\text{disk}},0.8\left|{\stackrel{.}{ℳ}}_{\text{g}}\right|\right].$$(54)

4.2 Evolution of distant cores

We first analyze the gas accretion applied to the furthest cores found in Fig. 8 of the set of simulations in Table 2 when the protoplanet forms at t_0,p = 0.2 Myr. In Fig. 9, we plot the core position vs the core formation time. We also display the final location at the end of the disk lifetime of 5 Myr of these cores after they have accreted gas. The figure shows clearly that when St = 0.03 or 0.06, the cores form at early epochs, causing the gas giants to migrate tens of AU before the disk lifetime. In the case of a large disk, the cores that form at around ~ 100 AU end up becoming gas giants close to 10 AU. As these cores form when Σg is still high, which the migration rate scales linearly with (see Eq. (49)), after they reach M_iso they still undergo very significant migration despite the gap-opening. When St = 0.01, the core forms at later epochs, and consequently the protoplanets undergo less (but still significant) migration when accreting gas. Thus, since the most significant migration occurs within the first Myr, these results remain consistent, even for a shorter disk lifetime of 3 Myr.

The formation of wide-orbit gas giants further out than 10 AU is clearly very challenging: the high M_iso in the outer regions requires a high surface density of pebbles to grow and this can only occur in disks with strong but short-lasting pebble fluxes. Hence, the core forms at early epochs and undergoes a long migration path due to the high surface density of the gas.

According to direct imaging surveys, the occurrence of stars hosting at least one giant planet with a mass between 0.5 and 14 M_Jup orbiting at 20–300 AU is only around 1% (Vigan et al. 2017). However, the current model we have developed does not adequately account for the formation of these distant gas giants. We go on to describe an alternative pathway for gas accretion below, which entails pebble flux decay and we investigate the feasibility of explaining wide-orbit gas giants via this pathway.

Fig. 10 Two pathways for gas accretion based on the requirement that pebble heating must cease for gas accretion to commence. (1) The pebble isolation pathway: the protoplanet reaches Miso and the pebbles are trapped at the outer edge of the gap. (2) The pebble decay pathway: the pebble flux decays due to radial drift. If the protoplanet previously attained sufficient mass, its envelope contracts by radiative heat loss, resulting in gas accretion. The core mass distinguishes giant planets formed by the two pathways.

4.3 The pebble decay pathway for gas accretion

The pebble isolation mass is very high in the outer regions; for example, to reach M_iso at 50 AU the protoplanet must form a core of ~ 50 M_⊕. During pebble accretion, the impacting pebbles heat the atmosphere, preventing contraction and efficient gas accretion (Lambrechts et al. 2014). Hence, the conventional pathway for gas accretion is that it is necessary to reach M_iso to cool down the envelope and accrete gas. Nevertheless, when the pebble flux reduces with time due to the radial drift of pebbles, a protoplanet that did not reach M_iso can still accrete gas efficiently. We shall denote the two distinct formation channels as the “pebble isolation pathway” and the “pebble decay pathway”, as illustrated in Fig. 10.

We use a simple description to implement the pebble decay pathway: if the time required for doubling the mass of a protoplanet via pebble accretion is longer than a certain threshold time, τ_th, pebble accretion stops and the protoplanet starts accreting gas. The mass-doubling timescale is $$\tau \approx \frac{M}{\dot{M}},$$(55)

where M is the mass of the protoplanet and $\dot{M}$ its growth rate via pebble accretion. When the pebble flux decays, $\dot{M}$ decreases and, consequently, τ will approach τ_th. Even small protoplanets that have experienced only limited growth can fulfill τ > τ_th. However, when M is very small, gas accretion is inefficient and the protoplanet will not grow within the lifetime of the protoplanetary disk.

The threshold time, τ_th, is a priori unknown and it will depend on the mass loading of heavy elements in the envelope, as it might take longer to cool down if the envelope is highly polluted (Ikoma et al. 2000). In Lambrechts et al. (2014), these authors calculated the minimal accretion rates required to sustain a stable gas envelope and from their calculations, we infer that τ_th is likely in the range between 10 and 100 Myr. We take τ_th = 10 Myr as the fiducial value.

We show in Fig. 11 the pathway that protoplanets would take depending on their initial position, r₀, and formation time, t_0,p, for different Stokes numbers when the disk lifetime is extended up to 5 Myr. For St = 0.03 and 0.06, a protoplanet initially placed at 30 AU needs to form earlier than 0.75 and 0.5 Myr, respectively, to reach M_iso. For St = 0.01, the required formation time extends beyond 1.2 Myr. However, a protoplanet placed at 100 AU will never reach M_iso when St = 0.01; furthermore, when St = 0.03 and 0.06 the required formation times drop down to 0.3 Myr and 0.25 Myr. Above the line that separates the protoplanets that reach M_iso and the ones that do not, in the St = 0.03 and 0.06 scenarios there is a relatively narrow region where the body can efficiently accrete gas³ due to pebble flux decay.

Figure 12 displays the evolution of 1000 protoplanets with initial conditions randomly chosen from the ranges shown in Fig. 11. We compare the core masses with the total metal amount of the giant planets in the Solar System: Jupiter from 25 to 45 M_⊕ (Wahl et al. 2017), Saturn from 18 to 20 M_⊕ (Mankovich & Fuller 2021), and Uranus and Neptune from 11 to 13 M_⊕ and from 13 to 15.5 M_⊕ respectively (Helled et al. 2011). A discussion on Solar System formation is given in Sect. 5.4. The main takeaway from Fig. 12 is that considering two paths for gas accretion results in two types of gas giants. Gas giants formed via the pebble isolation pathway have a metal-rich core. As their core forms at early epochs, they undergo a strong migration and end up orbiting ~ 5 AU at the furthest. In contrast, the final position of gas giants formed via the pebble decay pathway can extend beyond 5 AU up to 40 AU. They have a smaller core between 1.5 and 8 M_⊕. In addition, their formation is rather unusual, as the protoplanet must attain sufficient growth through pebble accretion to be able to accrete gas, while avoiding excessive migration and aligning with the decrease in pebble flux.

Fig. 11 Gas accretion pathway depending on the initial position and formation time of a protoplanet of 0.01 M⊕ in disks with different pebble Stokes numbers. We employ a disk lifetime of tf = 5 Myr, an initial metallicity of Z0 = 0.01 and a midplane turbulence of at = 10−4. The rest of the parameter values are listed in Table 1. Left: protoplanets cannot accrete gas via the pebble decay pathway when the pebbles are small and do not drift significantly. Center: under limited initial conditions, some gas giants form via the pebble decay pathway due to the short-lasting flux for a slightly higher value of St. Right: for even higher values of St the decay of the flux occurs earlier, and therefore, protoplanets must emerge at even earlier stages to reach Miso. The pebble decay pathway is less pronounced in this case.

Fig. 12 Population plot of core formation for 1000 protoplanets randomly initialized from Fig. 11 (top). We indicate Miso with gray dashed-dotted lines. We zoom in in the outer regions where the protoplanets reach Miso, and we indicate the metal content of the giant planets in the Solar System (see Sect. 5.4 for discussion). A medium or high pebble Stokes number is required to form a planet with JupiterΫs metallicity. Protoplanets that accrete gas via the pebble decay pathway (purple dots) acquire a significantly less massive core. Population plot including gas accretion for tf = 5 Myr (bottom). The pebble isolation mass Miso, 1 MJup and 10 MJup are indicated as reference lines, as well as the giant planets in the Solar System and the protoplanets PDS 70 b and c are indicated with red error bars (data from Wang et al. 2021). With a medium Stokes number and a short-lasting pebble flux, the formation of wide-orbit gas giants located up to 40 AU is possible due to the pebble decay pathway.

5 Implications and limitations

In this section, we establish connections between our primary findings and the observations of protoplanetary disks and planetary systems. Additionally, we highlight the limitations that require further exploration in future studies. First, we delve into the evolutionary patterns concerning the overall solid mass of protoplanetary disks. Subsequently, we investigate the plausibility of planetary cores being responsible for the observed gaps within these disks. Moreover, we investigate the potential origin of the PDS 70 system, which encompasses two wide-orbit Jovian protoplanets, as well as the formation of wide-orbit gas giants observed in mature planetary systems. Lastly, we conclude by exploring the prospective scenarios of the formation of the Solar System.

5.1 Solid mass evolution in protoplanetary disks

Our analytical solution to the pebble flux problem allows us to also express the temporal evolution of the total solid mass as: $$M_{\text{S}}(t)={\displaystyle {\int }_0^{\infty }2}\pi r{\Sigma }_{\text{p}}(r,t)\text{d}r=M_{\text{s},0}{T}^{-\frac{1}{2(2-\gamma )}\left(1+\frac{2}{3}\frac{\chi \text{St}}{\alpha }\right)},$$(56)

where the initial solid mass is $M_{\text{s},0}=M_{\text{g},0}\cdot Z_0=\frac{2}{3}\frac{{\stackrel{.}{ℳ}}_{g,0}}{v_1}\frac{R_1^2}{(2-\gamma )}{Z}_0$, the pebble surface density, Σ_p, is given by Eq. (40), and the rest of the parameters are specified in Sect. 2. This expression is valid for t ≥ t₀, t₀ being the time at which the solids grow up to the fragmentation or drift limit at R1. We ignore the dust mass loss in the earliest phases (t < t₀) and, therefore, a decrease in the initial metallicity Z₀ needs to be considered to apply properly to large disks of, for instance, 300 AU.

Figure 13 compares the solid mass evolution with the masses of Class 0 and I disks in the star-forming region Perseus derived by Tychoniec et al. (2020), as well as for Class II sources in Lupus derived by Ansdell et al. (2016). Their estimated median values for the dust mass in Class 0 and I disks are below the ones we calculated for a disk with R₁ = 100 AU. However, if we assume a more typical disk size of R₁ = 30 AU, the analytical expression lines up with the observed masses. This preliminary comparison is done only for Sun-like stars and the further exploration of parameter variations is left to future studies (see also Appelgren et al. 2023).

Fig. 13 Comparison between the evolution of the solid mass reservoir according to Eq. (56) for different initial disk sizes R1. We indicate the median dust masses derived by Tychoniec et al. (2020) for Class 0 and Class I sources in Perseus, and for Class II sources in Lupus by Ansdell et al. (2016). The dashed lines indicate the median value for each evolutionary class. We assume that the evolutionary phases Class 0 and I last approximately 0.1 and 0.5 Myr, respectively (Dunham et al. 2014). We assume that Z0 = 0.01 remains constant until solids grow to the drift or fragmentation limit at R1 as in Fig. 2. For R1 = 10, 30,100 AU, we estimate t0 = 0.007, 0.03, 0.2 Myr and St = 0.02, 0.03, 0.03, respectively. Gray dotted lines indicate solid mass decrease solely due to gas accretion onto the star, i.e., Ms(t) = Z0Mg(t).

5.2 Considering whether distant core formation explains the observed gaps in disks

The emergence of the observed substructures in the outer regions of protoplanetary disks has raised questions about their origin, since gas giants are rarely encountered in wide orbits in mature planetary systems (see discussion in Drążkowska et al. 2023). Nevertheless, in all our scenarios, once the cores form in the outer regions via the pebble isolation pathway (see Fig. 10), they migrate tens of AUs until the disk dissipates. That could explain why we observe multiple substructures beyond 10 AU that differ from an orbital location of < 10 AU, where most giant exoplanets are situated (see Fig. 1 in Lodato et al. 2019).

Bae et al. (2023) analyzed a sample of 62 protoplanetary disks observed at NIR and/or mm wavelengths containing rings and gaps, finding a maximum frequency of these substructures to occur at 20–50 AU. However, some rings were also found further than ~ 100 AU from their central star (see their Fig. 3d). When it comes to determining whether we can explain these substructures, we assume that protoplanets carve out gaps once they reach M_iso and neglect the plausible impact of fast migration regarding the gap-opening (Kanagawa et al. 2020). We hereby find that the furthest gap-forming cores form at distances from 20 AU to 80 AU in the most realistic scenarios (see Fig. 8). We note that we only considered the t_0,p = 0.7 Myr case for R₁ = 300 AU, since it takes approximately ~ 0.7 Myr for the solids to grow to the drift-limit at 300 AU. Overall, we find that the formation of distant cores is determined by the disk characteristics and formation of the initial protoplanet.

Disk parameters. Our fiducial disk model of Z₀ = 0.01 and R₁ = 100 AU contains an initial pebble mass reservoir of approximately 570 M_⊕. In this case, we see gap formation as far out as 30 AU. In order to form cores up to 50 AU or 80 AU, we find necessary either higher metallicity of Z₀ = 0.02 or larger disk size of R₁ = 300 AU, respectively, which both result in a higher solid mass reservoir (see also previous study by Ndugu et al. 2019). In addition, we find that weak turbulence enhances distant core formation and that the outcome depends strongly on the pebble Stokes number. When there is a strong short-lasting flux, distant cores form only at the early epoch (t < 1 Myr). On the contrary, less distant cores can form at a later epoch with a weak but long-lasting pebble flux instead (see Fig. 9). Protoplanets formed via the pebble decay pathway by a short-lasting pebble flux are also expected to carve out gaps at later stages after undergoing significant gas accretion (beyond 1 Myr or 2.5 Myr for St = 0.06 or St = 0.03, respectively). Such gaps may exhibit less observational prominence, since they would form only after significant depletion of the solid mass.

Initial protoplanets. we find that distant and rather early-formed initial protoplanets are required to form wide-orbit cores (see Fig. 11). The formation time is further constrained if the initial protoplanet’s mass is lower. These initial protoplanets could directly form by streaming instability if a local enrichment in metallicity of Z ~ 0.015 and a high Stokes number up to at least St = 0.01 − 0.1 are satisfied (Johansen et al. 2009; Bai & Stone 2010). A more recent study carried by Li & Youdin (2021) found that the streaming instability can trigger planetesimal formation across a broader range of parameters, even extending to subsolar metallicities. Therefore, planetesimal formation might also take place in the outer regions. The water ice line, where the inward drifting pebbles sublimate, is a favorable location for fulfilling these criteria (e.g., Ros & Johansen 2013; Ida & Guillot 2016; Drążkowska & Alibert 2017; Schoonenberg et al. 2018). If the ice lines of more volatile species (e.g., N₂, CO) are also favorable locations for forming planetesimals, this could explain the origin of distant protoplanets (Qi et al. 2013).

Finally, it has been suggested that envelope pollution due to the sublimation of accreted pebbles could decrease the critical metal mass for runaway gas accretion (e.g., Lambrechts et al. 2014; Brouwers & Ormel 2020). While the contraction of a polluted envelope could explain why ice giants like Uranus and Neptune were able to accrete a minor fraction of H₂/He below the pebble isolation mass, such a contraction among the polluted envelopes is not generally relevant for understanding the formation of gas giants in our study.

5.3 Explaining PDS 70 b and c and other wide-orbit gas giants

According to our model, wide-orbit gas giants should be rare, since the protoplanet must attain sufficient growth through pebble accretion to be able to accrete gas while avoiding excessive migration and aligning with the decrease in pebble flux. With the inclusion of the pebble decay pathway for gas accretion, Fig. 12 demonstrates that it is possible to explain the existence of gas giants such as PDS 70 b, located at ${20.8}_{-1.1}^{+1.3}\text{AU}$, and PDS 70 c, located at ${34.3}_{-3.0}^{+4.6}\text{AU}$ (Wang et al. 2021). While the scenario of both protoplanets forming via the pebble decay pathway appears to be less common, it suggests the possibility that this pathway may have played a role (see also Jiang & Ormel 2023, for an alternative formation channel in rings).

Our current study’s main limitation when reproducing the PDS 70 system is that we consider individual protoplanets while, in reality, simultaneously growing bodies could interact with each other. Such interactions can lead to scattering processes that hinder their growth (Levison et al. 2015; Bitsch et al. 2019), or reduce material delivery to inner bodies due to the gap-opening of an outer protoplanet (Weber et al. 2018, but see also Stammler et al. 2023). In addition, having a pair of giant planets could slow down their joint migration (Griveaud et al. 2023). Further studies are required to understand how gravitational interaction between simultaneously evolving bodies impacts the viability of the pebble decay pathway.

Overall, we reproduced gas giants up to 40 AU in an initial disk size of R₁ = 100 AU. Subsequent studies may explore whether a larger disk and the inclusion of photoevaporation (which could slow down the late migration) can reproduce giant planets even further out or whether other mechanisms, such as formation in rings or gravitational instability, would have to be invoked to explain such systems as HR 8799 (Marois et al. 2010).

5.4 Explaining the Solar System

In our simulations, most of the cores accrete solid material from the outer regions and end up as cold gas giants. This implies that the giant planets of the Solar System may have started their formation at distant locations, as suggested first by Bitsch et al. (2015) and explored further by Pirani et al. (2019). This distant formation is actually in line with the super-solar nitrogen abundance of Jupiter, since most of the nitrogen freezes outside of the N₂ snowline beyond 30 AU (Bosman et al. 2019; Öberg & Wordsworth 2019). In addition, the metal content of Jupiter is between 25 M_⊕ and 45 M_⊕ (Wahl et al. 2017). If most of the metal content is accreted during core formation, then according to Fig. 12, the protoplanet must reach M_iso between ~24 AU and 46 AU. Since modeling the formation of such massive cores growing solely via pebble accretion is challenging, late planetesimal accretion is usually deemed necessary to complete the total metal amount of Jupiter (proposed by Shiraishi & Ida 2008). However, Eriksson et al. (2022) demonstrated that during the gas accretion phase of the protoplanet, the accretion of planetesimals formed at planetary gap edges is a rather inefficient process. Therefore, an alternative mechanism for the late solid enrichment or the formation of a more massive core is necessary. In Fig. 8, we could reproduce within specific scenarios cores as massive as Jupiter’s via pebble accretion alone. According to Fig. 12, however, most of these bodies would become more massive than Jupiter once the process of gas accretion is terminated. This might be due to our simplified gas accretion prescription. Effects that have not been considered in our model such as photoevaporation (e.g., Ercolano & Rosotti 2015) could also hamper the gas accretion process.

Regarding the formation of Saturn, the estimated metal content falls within the range of 18 M_⊕ −20 M_⊕ (Mankovich & Fuller 2021), suggesting that it may have also reached M_iso. It is plausible that a weak long-lasting pebble flux could account for the formation of gas giants even at distances of approximately 10 AU (see Fig. 12). However, a stronger pebble flux is required to reach Jupiter’s metal content. Consequently, a formation scenario of both Jupiter and Saturn can be described in a model with a long-lasting pebble flux with St = 0.01 and a higher solid mass reservoir with Z₀ = 0.02 (see Fig. 8). Jupiter, being the first to form when the pebble flux was stronger, acquired a higher metal content and experienced more inward migration. On the other hand, Saturn could have gradually accreted from a less intense flux, ultimately settling further away. An alternative possibility is that Saturn formed through the pebble decay pathway, yet this seems less probable given its considerable metal content. Nevertheless, it is possible for protoplanets to accumulate additional metals through alternative mechanisms, such as the accretion of gas enriched in volatiles deposited by drifting pebbles (Schneider & Bitsch 2021a,b). Finally, for Uranus and Neptune, the most likely scenario is that they did not undergo efficient gas accretion, as they should have never reached M_iso (Lambrechts et al. 2014). This notion aligns well with the proposed formation of Jupiter and Saturn via a long-lasting pebble flux, as the lack of depletion would prevent Uranus and Neptune from accreting gas via the pebble decay pathway.

6 Summary

In this paper, we report the discovery of a new analytical expression to describe the temporal decay of the pebble flux. We first derived two analytical forms of the evolution of the metallic-ity, provided in Eqs. (32) or (36), under the assumptions that, respectively, the pebble Stokes number St or St_χ ≡ St ⋅ χ (where χ is the logarithmic gas pressure gradient) are constant both in time and space. We advocate to use constant St_χ as this agrees best with numerical simulations that include both radial drift and fragmentation in limiting the dust size. We demonstrated that for particles with St ≳ 0.01, core growth rates are significantly overestimated when using a simplified model that assumes a constant pebble-to-gas ratio at all radii and times.

We then used our derived pebble flux model to study the formation of distant planetary cores that reach the pebble isolation mass and therefore are able to carve out gaps in the outer regions of protoplanetary disks. In a large disk of 100 AU in size, we found that a Moon-sized protoplanet, that emerges within the first ~0.5 Myr and is located beyond 50 AU, can grow to become the core of a gas giant at 20–50 AU (see Fig. 8). The most distant cores form when there is a strong, but short-lasting, pebble flux with a Stokes number of ≳0.03 (see Sect. 3.2). An initial metallicity of ≳0.01 and low turbulence of ≲ 10⁻⁴ also promote the formation of distant cores. As these cores form early, they undergo a fast inward migration while they accrete gas, and by the end of the disk lifetime, they become giant planets orbiting at radii within 10 AU. In larger disks (e.g., R₁ = 300 AU), cores more massive than 50 M_⊕ could form beyond 50 AU, but it is still a challenge to explain the formation of planetary cores beyond 80 AU (see Sect. 5.2).

We also explored an alternative pathway for triggering gas accretion that we named the “pebble decay pathway” (see Fig. 10). We have demonstrated that this pathway could explain the formation of wide-orbit gas giants such as PDS 70 b and c. This pathway is only possible when the pebble flux decays significantly before the disk lifetime, implying that pebbles must grow up to large sizes (St ≳ 0.03) to undergo fast radial drift and deplete. Given the appropriate disk parameters and initial conditions of the protoplanets, the location of these gas giants ranges from 1 AU up to 40 AU (see Fig. 12). According to our model, to form a wide-orbit gas giant the protoplanet must have attained sufficient mass by the time the pebble flux has decayed in order to accrete gas efficiently, but not be too massive to avoid excessive migration. Hence, these planets are rare, in agreement with the low occurrence of ~1% of distant gas giants from direct imaging surveys. Since they never reach M_iso, we predict a low bulk metallicity content for these wide-orbit gas giants.

Acknowledgements

N.G. thanks Federico Finkel for his comments on the mathematical derivations. The authors also thank the anonymous referee for the comments that helped to improve the manuscript. A.J. is supported by the Swedish Research Council (Project grant 2018-04867), the Danish National Research Foundation (DNRF Chair grant DNRF159), and the Knut and Alice Wallenberg Foundation (Wallenberg Academy Fellow grant 2017.0287). A.J. further thanks the European Research Council (ERC Consolidator grant 724 687-PLANETESYS), the Göran Gustafsson Foundation for Research in Natural Sciences and Medicine, and the Wallenberg Foundation (Wallenberg Scholar KAW 2019.0442) for research support. M.L. acknowledges the ERC starting grant 101041466-EXODOSS.

Appendix A Derivation of the analytical metallicity for a constant St

We assume that St is constant and that $Z(\tilde{r},1)=Z_0$. Substituting $\chi \left(\tilde{r},T\right)$ from Eq. (12), we get: $$\begin{array}{l}(\tilde{r},T)=\frac{2}{3}\frac{\text{St}}{\alpha }\left({\chi }_0+(2-\gamma )\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right)\\ =b_0\left(1+\frac{2-\gamma }{{\chi }_0}\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right),\end{array}$$(A.1)

where $b_0=\frac{2}{3}\frac{{\chi }_0\text{St}}{\alpha }$. The spatial derivative of $b\left(\tilde{r},T\right)$ is: $$\frac{\partial b}{\partial \tilde{r}}=\frac{2}{3}\frac{\text{St}}{\alpha }{(2-\gamma )}^2\frac{{\tilde{r}}^{(1-\gamma )}}{T}=b_0\frac{{(2-\gamma )}^2}{{\chi }_0}\frac{{\tilde{r}}^{(1-\gamma )}}{T}.$$(A.2)

Replacing $\frac{\partial b}{\partial \tilde{r}}$ in the general continuity Eq. (31) and multiplying the equation by the term ${\left[(2-\gamma ){\tilde{r}}^{(1-\gamma )}\right]}^{-1}$, we get the governing $$\begin{array}{l}\frac{1}{2}\left[\frac{1+b_0}{2-\gamma }\frac{1}{{\tilde{r}}^{(1-\gamma )}}+\left(\frac{b_0}{{\chi }_0}-2\right)\frac{\tilde{r}}{T}\right]\frac{\partial Z}{\partial \tilde{r}}-(2-\gamma )\frac{\partial Z}{\partial T}\\ =\frac{Z}{T}\frac{b_0}{2}\left[1-\frac{2-\gamma }{{\chi }_0}\left(1-\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right)\right].\end{array}$$(A.3)

This equation is a first-order linear PDE, which, in turn, can simply be expressed as: $$A(\tilde{r},T)\frac{\partial Z}{\partial \tilde{r}}+B(\tilde{r},T)\frac{\partial Z}{\partial T}=C(\tilde{r},T,Z).$$(A.4)

The corresponding Lagrange-Charpit system is $$\frac{d\tilde{r}}{A(\tilde{r},T)}=\frac{dT}{B(\tilde{r},T)}=\frac{dZ}{C(\tilde{r},T,Z)}.$$(A.5)

From the first equality, we get: $$\frac{d\tilde{r}}{\frac{1}{2}\left[\frac{1+b_0}{2-\gamma }\frac{1}{{\tilde{r}}^{(1-\gamma )}}+\left(\frac{b_0}{{\chi }_0}-2\right)\frac{\tilde{r}}{T}\right]}=-\frac{dT}{2-\gamma }.$$(A.6)

By making the change of variable $x={\tilde{r}}^{(2-\gamma )}$, $$\frac{dx}{dT}=(2-\gamma ){\tilde{r}}^{1-\gamma }\frac{d\tilde{r}}{dT}=-\frac{1}{2}\left[\frac{1+b_0}{2-\gamma }+\left(\frac{b_0}{{\chi }_0}-2\right)\frac{x}{T}\right].$$(A.7)

This equation is a homogeneous first-order ODE, meaning that it takes the form $\frac{dx}{dT}=f\left(\frac{x}{T}\right)$ and therefore it can be solved by: $$\frac{dT}{T}=\frac{du}{F(u)-u},$$(A.8)

where $u=\frac{x}{T}$ and $F(u)=\frac{dx}{dT}$. Integrating the equation and substituting $x={\tilde{r}}^{(2-\gamma )}$, we get the solution $$\ln T=-\frac{2{\chi }_0}{b_0}\ln \left(\frac{1+b_0}{2(2-\gamma )}+\frac{b_0}{2{\chi }_0}\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right)+p,$$(A.9)

where p is an invariant. Rearranging the equation, we have: $$\tilde{r}(T)={\left[-\frac{{\chi }_0}{2-\gamma }\left(1+\frac{1}{b_0}\right)T+p{T}^{-\frac{b_0}{2{\chi }_0}+1}\right]}^{\frac{1}{2-\gamma }},$$(A.10)

or $$p=\left[\frac{{\chi }_0}{2-\gamma }\left(1+\frac{1}{b_0}\right)+\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right]T^{\frac{b_0}{2{\chi }_0}}.$$(A.11)

From the second equality in Eq. (A.5), we get: $$-\frac{dT}{2-\gamma }=\frac{dZ}{\frac{Z}{T}\frac{b_0}{2}\left[1-\frac{2-\gamma }{{\chi }_0}\left(1-\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right)\right]}.$$(A.12)

By replacing $\tilde{r}$ from Eq. (A. 10) and rearranging the equation, $$\left(\frac{1}{2(2-\gamma )}+\frac{b_0}{2{\chi }_0}\right)\frac{dT}{T}-\frac{b_0}{2{\chi }_0}p{T}^{-\frac{b_0}{2{\chi }_0}-1}dT=\frac{dZ}{Z}.$$(A.13)

Integrating the ODE, we get: $$\ln T^{\frac{1}{2(2-\gamma )}+\frac{b_0}{2{\chi }_0}}+p{T}^{-\frac{b_0}{2{\chi }_0}}+f(p)=\ln Z.$$(A.14)

By rearranging and replacing the invariant p from Eq. (A.11), the general solution of the PDE is: $$\begin{array}{l}Z(\tilde{r},T)=T^{\frac{1}{2(2-\gamma )}+\frac{b_0}{2{\chi }_0}}\exp \left[\frac{{\chi }_0}{2-\gamma }\left(1+\frac{1}{b_0}\right)+\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right]\\ \times f\left(\left[\frac{{\chi }_0}{2-\gamma }\left(1+\frac{1}{b_0}\right)+\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right]T^{\frac{b_0}{2{\chi }_0}}\right).\end{array}$$(A.15)

To get the form of $Z(\tilde{r},T)$ that fulfills the initial condition $Z(\tilde{r},1)=Z_0$, since we have: $$\begin{array}{l}Z(\tilde{r},1)=\exp \left[\frac{{\chi }_0}{2-\gamma }\left(1+\frac{1}{b_0}\right)+{\tilde{r}}^{(2-\gamma )}\right]\\ \times f\left(\left[\frac{{\chi }_0}{2-\gamma }\left(1+\frac{1}{b_0}\right)+{\tilde{r}}^{(2-\gamma )}\right]\right)=Z_0,\end{array}$$(A.16)

f (p) must take the form f (p) = Z₀e^−p. Setting everything together, we obtain the solution: $$\begin{array}{l}Z(\tilde{r},T)=Z_0{T}^{\frac{1}{2(2-\gamma )}+\frac{b_0}{2{\chi }_0}}\\ \times \exp \left\{-\left[\frac{{\chi }_0}{2-\gamma }\left(1+\frac{1}{b_0}\right)+\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right]\cdot \left[T^{\frac{b_0}{2{\chi }_0}}-1\right]\right\}.\end{array}$$(A.17)

By replacing b₀, we have: $$\begin{array}{l}Z(\tilde{r},T)=Z_0{T}^{\frac{1}{2(2-\gamma )}+\frac{\text{St}}{3\alpha }}\\ \times \exp \left\{-\left[\frac{1}{2(2-\gamma )}\left(2{\chi }_0+\frac{3\alpha }{\text{St}}\right)+\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right]\cdot \left[T^{\frac{\text{St}}{3\alpha }}-1\right]\right\}.\end{array}$$(A.18)

Appendix B Derivation of the analytical metallicity for a constant St_χ

Here, we follow the same procedure as in previous appendix.

Multiplying the general equation by the term ${\left[(2-\gamma ){\tilde{r}}^{(1-\gamma )}\right]}^{-1}$ we get: $$\frac{1}{2}\left[\frac{1+b_0}{2-\gamma }\frac{1}{{\tilde{r}}^{(1-\gamma )}}-2\frac{\tilde{r}}{T}\right]\frac{\partial Z}{\partial \tilde{r}}-(2-\gamma )\frac{\partial Z}{\partial T}=\frac{Z}{T}\frac{b_0}{2}.$$(B.1)

The corresponding Lagrange-Charpit system is: $$\frac{d\tilde{r}}{\frac{1}{2}\left[\frac{1+b_0}{2-\gamma }\frac{1}{{\tilde{r}}^{(1-\gamma )}}-2\frac{\tilde{r}}{T}\right]}=-\frac{dT}{2-\gamma }=\frac{dZ}{\frac{Z}{T}\frac{b_0}{2}}.$$(B.2)

The second equality can be directly integrated to get: $$Z(\tilde{r},T)=T^{-\frac{b_0}{2(2-\gamma )}}f(p).$$(B.3)

We can rewrite the first equality from Eq. (B.2) as: $$\frac{d\tilde{r}}{dT}=-\frac{1}{2(2-\gamma )}\left[\frac{1+b_0}{2-\gamma }\frac{1}{{\tilde{r}}^{(1-\gamma )}}-2\frac{\tilde{r}}{T}\right].$$(B.4)

Making again the change in the variable $x={\tilde{r}}^{(2-\gamma )}$, $$\frac{dx}{dT}=(2-\gamma ){\tilde{r}}^{(1-\gamma )}\frac{d\tilde{r}}{dT}=-\frac{1}{2}\left[\frac{1+b_0}{2-\gamma }-2\frac{x}{T}\right],$$(B.5)

which is also an homogeneous first-order ODE that can be solved as described in Eq. (A.8). The solution is then: $$p=\frac{1}{2-\gamma }\ln T+\frac{2}{1+b_0}\frac{{\tilde{r}}^{(2-\gamma )}}{T}.$$(B.6)

The general solution for the PDE is therefore: $$Z(\tilde{r},T)=T^{-\frac{b_0}{2(2-\gamma )}}f\left(\frac{1}{2-\gamma }\ln T+\frac{2}{1+b_0}\frac{{\tilde{r}}^{(2-\gamma )}}{T}\right).$$(B.7)

To fulfill the initial condition $Z(\tilde{r},T)=Z_0$, since $$Z(\tilde{r},1)=f\left(\frac{2}{1+b_0}{\tilde{r}}^{(2-\gamma )}\right)=Z_0\text{, }$$(B.8)

f (p) must take the form of f (p) = Z₀. Hence, the solution is $$Z(\tilde{r},T)=Z_0{T}^{-\frac{b_0}{2(2-\gamma )}},$$(B.9)

and by replacing $b_0=\frac{2}{3}\frac{\chi \cdot \text{St}}{\alpha }$, we get $$Z(\tilde{r},T)=Z_0{T}^{-\frac{1}{2-\gamma }\frac{\chi \cdot \text{St}}{3\alpha }}.$$(B.10)

Appendix C An alternative method for determining the evolving metallicity

We describe an alternative way of deriving the evolution of the metallicity. First, we use simplified expressions for the gas surface density and gas flux, commonly used to describe the inner regions, $${\Sigma }_{\text{g}}(r,t)=\frac{{\stackrel{.}{ℳ}}_{g,0}}{3\pi v_1{\tilde{r}}^{\gamma }}{T}^{-\frac{5/2-\gamma }{2-\gamma }}\text{,}$$(C.1) $${\stackrel{.}{ℳ}}_g(t)={\stackrel{.}{ℳ}}_{g,0}{T}^{-\frac{5/2-\gamma }{2-\gamma }}.$$(C.2)

Recalculating $\frac{\partial \ln P}{\partial \ln r}$ with the simplified Σ_g(r, t) (see Eq. (12)), we get that χ = χ₀. The ratio between the radial velocities from Eqs. (10) and (13) is: $$\frac{{\text{υ}}_{\text{r},\text{p}}}{{\text{υ}}_{\text{r},\text{g}}}=\frac{1+\frac{2}{3}\frac{\chi \cdot \text{St}}{\alpha }}{1+{\text{St}}^2},$$(C.3)