Planet formation via pebble accretion in externally photoevaporating discs

ABSTRACT

We demonstrate that planet formation via pebble accretion is sensitive to external photoevaporation of the outer disc. In pebble accretion, planets grow by accreting from a flux of solids (pebbles) that radially drift inwards from the pebble production front. If external photoevaporation truncates the outer disc fast enough, it can shorten the time before the pebble production front reaches the disc outer edge, cutting off the supply of pebble flux for accretion, hence limiting the pebble mass reservoir for planet growth. Conversely, cloud shielding can protect the disc from strong external photoevaporation and preserve the pebble reservoir. Because grain growth and drift can occur quickly, shielding even on a short time-scale (<1 Myr) can have a non-linear impact on the properties of planets growing by pebble accretion. For example, a |$10^{-3}\, \mathrm{ M}_{\oplus }$| planetary seed at 25 au stays at 25 au with a lunar mass if the disc is immediately irradiated by a 10³ G₀ field, but grows and migrates to be approximately Earth-like in both mass and orbital radius if the disc is shielded for just 1 Myr. In NGC 2024, external photoevaporation is thought to happen to discs that are <0.5 Myr old, which coupled with the results here suggests that the exact planetary parameters can be very sensitive to the star-forming environment. Universal shielding for time-scales of at least |${\sim} 1.5\,$| Myr would be required to completely nullify the environmental impact on planetary architectures.

planets and satellites: formation, protoplanetary discs, circumstellar matter

1 INTRODUCTION

There are now more than 5000 confirmed exoplanets exhibiting very diverse properties and planetary system architectures (e.g. Gaudi, Meyer & Christiansen 2021). Understanding this diversity of planets requires understanding planet formation and evolution. However, this is not straightforward as planet formation involves the complex interplay of many processes (see reviews by e.g. Drazkowska et al. 2022; Raymond & Morbidelli 2022; Emsenhuber, Mordasini & Burn 2023). One promising model for the formation of terrestrial planets and cores of giant planets (in the core accretion scenario) is pebble accretion, in which mm–cm-sized pebbles are mainly responsible for planetary growth by their accretion on to larger planetesimals and cores (Lambrechts & Johansen 2012, 2014). The pebble accretion mechanism is an efficient way of resolving issues associated with the traditional core accretion model: the growth barrier when particles reach mm–cm due to radial drift (Weidenschilling 1977) and slow planetesimal accretion with multikilometre-sized bodies (Pollack et al. 1996; Coleman & Nelson 2014, 2016a, b). Pebbles are accreted on to growing planetary embryos at a much faster rate than planetesimals due to the greater kinetic energy dissipation by gas drag when pebbles enter the gravitational reach of a core (Ormel & Klahr 2010; Lambrechts & Johansen 2012). There has also been some success in being able to form compact terrestrial planetary systems such as Trappist-1 via pebble accretion (Ormel, Liu & Schoonenberg 2017; Coleman et al. 2019).

One potential factor contributing to planet diversity is the birth environment. Most stars form in large stellar clusters with high stellar density (Lada & Lada 2003; Fatuzzo & Adams 2008; Lee & Hopkins 2020; Winter et al. 2020a), and there has been growing attention on how the dense cluster environment can influence the planet-forming disc evolution and possibly the final planetary systems. One way the stellar cluster environment can influence a protoplanetary disc is via star–disc gravitational encounters (for a recent review, see Cuello, Ménard & Price 2023). The other environmental impact on discs is due to ‘external’ ultraviolet (UV) irradiation by massive stars in the cluster (for a recent review, see Winter & Haworth 2022). This external irradiation drives a thermal wind from the outer disc in a process called external photoevaporation. In addition to depleting the disc mass, if material from the outer disc is removed by external photoevaporation more quickly than its resupply via viscous spreading, then the disc can be truncated (Clarke 2007; Rosotti et al. 2017; Eisner et al. 2018; Coleman & Haworth 2022).

There have been numerous direct observations of discs being affected by their environments. When the radiation field is sufficiently strong, the circumstellar disc is enshrouded in a cometary wind, with the cusp directed towards the exciting UV source. These ‘proplyds’ have been observed across a range of UV radiation field environments in the Orion Nebular Cluster and NGC 1977 (e.g. O’dell & Wen 1994; Henney & O’Dell 1999; Bally, Youngblood & Ginsburg 2012; Kim et al. 2016). There is also recent evidence of external photoevaporation happening at very early times in NGC 2024 (<0.5 Myr; Haworth et al. 2021), competing with even the earliest evidence for planet formation (Sheehan & Eisner 2018; Segura-Cox et al. 2020). One disc in σ Ori is observed with unusually high [O i] 6300 Å emission, which has also been interpreted as due to an external photoevaporative wind (Störzer & Hollenbach 1998; Rigliaco et al. 2009; Ballabio, Haworth & Henney 2023). Additionally, statistical comparisons of disc properties in clusters such as masses (Ansdell et al. 2017; Eisner et al. 2018), radii (Eisner et al. 2018; Boyden & Eisner 2020; Otter et al. 2021), and disc fractions (Reiter & Parker 2019; van Terwisga et al. 2020) also show evidence of external photoevaporation.

Analytical and semi-analytical models of the external photoevaporation of discs were developed after the discovery of proplyds (Johnstone, Hollenbach & Bally 1998; Störzer & Hollenbach 1999; Adams et al. 2004). Clarke (2007) coupled viscous evolution of discs with mass-loss rate calculations from semi-analytical models of photoevaporation. However, it is specifically far-ultraviolet (FUV) light that typically drives external photoevaporative winds, so to accurately calculate the mass-loss rate it is necessary to determine the temperature in the photodissociation region (PDR) from microphysics that cannot be solved analytically. Haworth et al. (2018a) developed a publicly available grid of mass-loss rates (called FRIED) computed with PDR hydrodynamics simulations, which has enabled the community to run on-the-fly modelling of disc evolution with external photoevaporation (e.g. Haworth et al. 2018b; Concha-Ramírez et al. 2019, 2021; Winter, Clarke & Rosotti 2019a; Winter et al. 2019b, 2020a, b; Sellek, Booth & Clarke 2020; Parker, Nicholson & Alcock 2021a; Concha-Ramírez, Wilhelm & Portegies Zwart 2023; Qiao et al. 2022; Wilhelm et al. 2023). Studies that included both effects of dynamical encounters and external photoevaporation of discs in dynamically evolving stellar clusters (Scally & Clarke 2001; Winter et al. 2018; Concha-Ramírez et al. 2019, 2021, 2023; Wilhelm et al. 2023) generally conclude that external photoevaporation is overall more significant in terms of causing mass-loss, but both mechanisms are likely to impact disc evolution, and the interplay between these processes is still not well understood. In this paper, we focus on planet formation in the context of external photoevaporation, but the additional role of encounters should be included in future.

In recent years, there have been efforts to more directly couple disc evolution calculations to star formation simulations. For example, Winter et al. (2019a), Concha-Ramírez et al. (2019, 2021), and Parker et al. (2021b) coupled n-body simulations of stellar cluster dynamics with external photoevaporation and disc evolution models. In those calculations, all stars/discs are assumed to have formed at the start of the simulation. They calculate the time-varying FUV field incident upon each disc (and hence the time-varying external photoevaporative mass-loss rate from FRIED) due to the aggregate of stars in the cluster, assuming geometric dilution of each star’s intrinsic FUV luminosity. However, this does not account for ongoing star formation, shielding by the star-forming cloud, or feedback processes such as ionizing radiation, winds, and supernovae that can act to disperse the cloud. Concha-Ramírez et al. (2023) coupled external photoevaporation with a simulation following the collapse of a giant molecular cloud in an SPH (smoothed particle hydrodynamics) hydrodynamical model that self-consistently formed stars, and found that ongoing star formation is necessary for massive discs to exist beyond the early cluster age. Recently, Qiao et al. (2022) and Wilhelm et al. (2023) coupled star cluster formation and feedback simulations with disc evolution models. These directly account for the impact of cloud shielding on the UV field incident upon discs and the dispersal of the shielding matter over time by feedback. They found that the effect of shielding in dense stellar clusters can prevent high external photoevaporative mass-loss at least in the early stage of disc evolution. However, typical shielding time-scales will depend on the nature of the star formation event, such as the cloud virial ratio, which was already known to determine how effectively and rapidly feedback can disperse the cloud (Dale, Ercolano & Bonnell 2012a, b, 2013a; Parker et al. 2021a; Qiao et al. 2022; Wilhelm et al. 2023).

Overall, it is now well recognized that in an evolving star-forming region, external photoevaporation and cloud shielding times can have significant impacts on disc evolution. However, what is much less studied is the impact of external photoevaporation on planet formation. Recently, Winter et al. (2022) found that the gas accretion and migration of wide-orbit giant planets in a disc can be suppressed by FUV-induced photoevaporation. There are also population synthesis models incorporating mass-loss in gas (interpolated from FRIED grid) and dust entrained through external photoevaporative winds (e.g. Burn et al. 2022; Emsenhuber et al. 2023). Burn et al. (2022) found that dust entrainment in photoevaporative winds can be an important mechanism for removing solids. However, one unclear but important aspect of all this is the competition between early planet formation (substructures in 0.5 Myr discs suggest possible early planet formation; Sheehan & Eisner 2018; Segura-Cox et al. 2020) and external photoevaporation. If a planet grows via pebble accretion, external photoevaporation could possibly limit the inwardly drifting pebble mass budget by depleting the outer disc. Conversely, cloud shielding could protect the disc for some time before being exposed to a high radiation field (Qiao et al. 2022; Wilhelm et al. 2023) that could give planet formation a head start. Therefore, the impact of external photoevaporation and shielding might affect not just the formation of wide-orbit giant planets, but also compact planetary systems. We therefore need to understand how important the shielding time-scale is for planet formation theoretically, and to observationally determine how long it is in practice. The objective of this paper is to provide the theoretical component. We couple a model of planet formation via pebble accretion in a viscously evolving disc with external photoevaporation induced by a parametrized time-varying FUV radiation field. We aim to isolate the impacts of external photoevaporation and shielding on planet formation by pebble accretion.

2 NUMERICAL METHOD

We investigate the formation of planets via pebble accretion in externally photoevaporating protoplanetary discs. To isolate the impact of external photoevaporation and cloud shielding on the planet formation process in a controlled way, we apply a state-of-the-art planet formation code with parametrized time-varying external irradiation field. The disc evolves viscously (Section 2.1) and is also subject to mass-loss via external photoevaporation due to the parametrized time-varying external FUV radiation field strength (Section 2.1.2). For simplicity, in each calculation we only consider a single planetary embryo in the disc. This can then evolve by accreting pebbles (Section 2.2) and a gaseous envelope (Section 2.3) while migrating in the type I migration regime (Section 2.4).

In total, we vary three parameters in this study: initial shielding time t_sh (0–3 Myr), the maximum FUV field strength that discs are exposed to after the initial shielding time F_FUV,max (10–10 000 G₀), and the initial semimajor axis of planetary embryos a_init (10–100 au). Each simulation has just a single planetary embryo in the disc, and is run for a total of 10 Myr, unless the planet migrates interior to 1|$\, {\rm au}$| or the disc mass drops below 0.1 M_⊕. For the purpose of isolating the impact of external photoevaporation, we do not consider internal photoevaporation at this stage (see Section 2.1.1 for detailed discussions).

2.1 Disc evolution

We adopt a one-dimensional (1D) viscous disc model where we evolve the gas surface density by solving the standard diffusion equation with an additional term to describe mass-loss by external photoevaporation:

$$\begin{eqnarray*} \dot{\Sigma }(r) = \frac{1}{r}\frac{\mathrm{ d}}{\mathrm{ d}r}\left[3 r^{1/2}\frac{\mathrm{ d}}{\mathrm{ d}r}(\nu \Sigma r^{1/2}) \right] - \dot{\Sigma }_{\mathrm{ PE}}(r), \end{eqnarray*}$$

(1)

where |$\dot{\Sigma }_{\mathrm{ PE}}(r)$| is the rate of change of surface density due to external photoevaporative wind induced by FUV radiation (see Section 2.1.1). We use the standard α model for the disc viscosity (Shakura & Sunyaev 1973)

$$\begin{eqnarray*} \nu = \alpha {c_{\mathrm{ s}}}^{2}/\Omega , \end{eqnarray*}$$

(2)

where c_s is the local sound speed, Ω is the angular velocity, and α is the viscosity parameter, which is set at α = 10⁻³ for all simulations.

For the initial gas surface density profile, we use the similarity solution of Lynden-Bell & Pringle (1974)

$$\begin{eqnarray*} \Sigma = \Sigma _{0}\left(\frac{R}{R_\mathrm{ C}}\right)^{-1}\exp {\left(-\frac{R}{R_\mathrm{ C}}\right)}, \end{eqnarray*}$$

(3)

where Σ₀ is the normalization constant set by the total disc mass (for a given R_C), and R_C is the scale radius, which sets the initial disc size. For all our simulations, we use an initial disc mass of 105 M_Jup (∼0.1M_*) and an initial scale radius R_C = 50 au. We also use a constant host star mass M_* = 1 M_⊙ and solve equation (1) over 2000 grid cells using an explicit finite difference scheme, adopting a non-uniform mesh grid, for which the grid spacing Δr scales with radius (Coleman & Nelson 2014). We set the local solids-to-gas ratio (Z₀) everywhere in the disc to 0.01 and note that the solids in the disc include both dust particles that contribute to the disc opacity and pebbles produced from dust coagulation (see Section 2.2).

We set an initial radial temperature profile as

$$\begin{eqnarray*} T_\textrm {{initial}} = T_{0}\left(\frac{R}{R_0}\right)^{-0.5}, \end{eqnarray*}$$

(4)

where R₀ = 1 au and T₀ is the temperature at R₀ (1 au), which is set to 280 K. The disc temperature after each time-step of the simulation is obtained by solving the following thermal equilibrium equation that describes the balance of heating from central star irradiation, heating from residual molecular cloud, viscous heating, and blackbody cooling using an iterative method:

$$\begin{eqnarray*} Q_\textrm {{irr}} + Q_{\nu } + Q_\textrm {{cloud}}-Q_\textrm {{cool}} = 0, \end{eqnarray*}$$

(5)

where Q_irr is the radiative heating rate due to the central star, Q_ν is the viscous heating rate per unit area of the disc, Q_cloud is the radiative heating due to the residual molecular cloud (with temperature 10 K), and Q_cool is the blackbody radiative cooling rate (see equations 4–10 in Coleman 2021, for further descriptions of the expressions in calculating the temperature). For reference, in Fig. A1 we show the surface density evolution for discs in different UV field strengths and with different shielding times.

2.1.1 External photoevaporation

For a disc in a dense stellar cluster, external photoevaporative winds can be launched by both EUV and FUV photons radiated from nearby massive stellar neighbours, though the FUV is generally expected to be dominant in setting the mass-loss rate (O’dell & Wen 1994; Johnstone et al. 1998; Adams et al. 2004; Facchini, Clarke & Bisbas 2016; Winter & Haworth 2022). In our simulations, the mass-loss rate due to external photoevaporation is calculated via interpolating over the FRIED grid (Haworth et al. 2018b). FRIED grid provides mass-loss rates for discs irradiated by FUV radiation as a function of the star/disc/FUV parameters. In our simulations, we determine the mass-loss rate at each time-step by linearly interpolating FRIED in three dimensions: disc size R_d, outer edge disc surface density Σ_out, and FUV field strength at this time F_FUV(t) (see Section 2.1.2 for details about the FUV field tracks used).

A cautionary note is that FRIED mass-loss rates are sensitive to the value of the disc outer radius R_d and the surface density there Σ_out. However, towards the disc outer edge the surface density can be a strong function of radius (for example, dropping off exponentially). Furthermore, if the disc outer radius is chosen at a point that would be optically thin to the external irradiation, the situation is unphysical since in reality radiation would propagate deeper into the disc and drive the wind from there. This means that the disc outer edge location should be at the optically thick/thin transition, where the interpolated mass-loss rate from FRIED is maximum. We therefore adopt that radius of maximum mass-loss R_max as the disc outer radius (see Sellek et al. 2020, particularly fig. 2 for more information).

We evaluate the FRIED mass-loss rate at each radius in the disc to determine R_max. We then extract the corresponding mass-loss from the disc via

$$\begin{eqnarray*} \dot{\Sigma }_{\textrm {ext,~FUV}}(r) = G_{\rm sm} \frac{\dot{M}_\rm {{ext}}(R_{\rm {max}})}{\pi (R^2_\rm {{d}}-{R_{\rm {\text {max}}}}^2)+A_{\rm sm}}, \end{eqnarray*}$$

(6)

where A_sm is a smoothing area equal to

$$\begin{eqnarray*} A_{\rm sm} = \dfrac{\pi (R_{\rm max}^{22}-(0.1 R_{\rm max})^{22})}{11R_{\rm max}^{20}} \end{eqnarray*}$$

(7)

and G_sm is a smoothing function

$$\begin{eqnarray*} G_{\rm sm} = \dfrac{r}{R_{\rm max}^{20}}. \end{eqnarray*}$$

(8)

Note that at this stage, to isolate the impact of external photoevaporation, we do not include internal photoevaporative winds (e.g. Ercolano & Pascucci 2017; Ercolano & Picogna 2022, for reviews). Internal photoevaporative winds can entrain small dust in the wind (Hutchison et al. 2016a; Hutchison, Laibe & Maddison 2016b; Hutchison & Clarke 2021) and eventually open up an inner hole on a time-scale that depends on the inward transport of gas and the strength of the radiation field (EUV, FUV, or X-ray) driving the wind (e.g. Owen et al. 2010; Owen, Clarke & Ercolano 2012; Picogna et al. 2019; Picogna, Ercolano & Espaillat 2021; Coleman & Haworth 2022; Sellek, Clarke & Ercolano 2022). Studying the combined interplay with these processes is left for future work.

2.1.2 Shielding time-scale prescription

We aim to investigate the effect of shielding time t_sh on planet formation by pebble accretion in a controlled way. Therefore, rather than follow Qiao et al. (2022) or Wilhelm et al. (2023) and use complicated time-varying FUV radiation fields from simulations of star formation, here we opt to use a parametrized time-varying FUV flux profile that transitions from a low FUV field F_FUV,0 to a high value F_FUV,max after a shielding period t_sh

$$\begin{eqnarray*} F_{\textrm {FUV}}(t) &=& F_{\textrm {FUV,0}} + \frac{1}{2}(F_{\textrm {FUV,max}}-F_{\textrm {FUV,0}}) \nonumber\\ &&\times \, \left(\tanh \left(\frac{t-t_\textrm {{shield}}}{t_{\textrm {trans}}}\right) + 1\right), \end{eqnarray*}$$

(9)

where t_trans = 5 × 10⁴ yr is a parameter that controls the time-scale over which the star/disc moves from a shielded to an unshielded environment. This rapid transition from shielded to irradiated represents the change in irradiation when stars/discs cease being embedded (Qiao et al. 2022). More complicated time variation of the UV could arise due to cluster dynamics (e.g. a star moving from a low UV region into a higher UV region), but including this in this study makes the parameter space unwieldy. We set F_FUV,0 = 10 G₀, and vary t_sh (0–3 Myr) and F_FUV,max (10–10 000 G₀) as parameters in our simulations. Fig. 1 illustrates some examples of our time-varying external FUV radiation parametrization.

$Examples of some of the time-varying FUV tracks used in the simulation described by equation (9); the colours indicate different shielding times, with blue, orange, green, and red representing tsh = 0.5, 1, 2, and 3 Myr, respectively, and the linestyles represent different values of FFUV,max, with the solid curves representing $F_{\textrm {FUV,~max}} = 10\,000\, \mathrm{ G}_0$, the dashed curves representing $F_{\textrm {FUV,~max}} = 1000\, \mathrm{ G}_0$, and dash–dotted lines representing $F_{\textrm {FUV,~max}} = 100\, \mathrm{ G}_0$.$

Figure 1.

Examples of some of the time-varying FUV tracks used in the simulation described by equation (9); the colours indicate different shielding times, with blue, orange, green, and red representing t_sh = 0.5, 1, 2, and 3 Myr, respectively, and the linestyles represent different values of F_FUV,max, with the solid curves representing |$F_{\textrm {FUV,~max}} = 10\,000\, \mathrm{ G}_0$|⁠, the dashed curves representing |$F_{\textrm {FUV,~max}} = 1000\, \mathrm{ G}_0$|⁠, and dash–dotted lines representing |$F_{\textrm {FUV,~max}} = 100\, \mathrm{ G}_0$|⁠.

Open in new tab Download slide

2.2 Pebble accretion

We use the model of pebble production and accretion implemented in Coleman (2021), which follows the models of Lambrechts & Johansen (2012, 2014). In this model of pebble production via dust growth, solids are separated into two populations: one of radially stationary small dust grains and one of larger inwardly drifting pebbles, formed from the coagulation of the vertically settling small dust grain population in the disc. At each disc radius, it takes a certain amount of time for the small ISM-like dust grains (about 0.1–1 |$\mu$|m; Weingartner & Draine 2001) to grow to a pebble size large enough to start drifting R_drift (⁠|$R_{\textrm {drift}}\, \approx$| 1–10 mm; Lambrechts & Johansen 2014), assuming a certain growth efficiency. In other words, at a certain time t, there is one particular disc radius r_g(t) that is called the pebble production front, where the dust particles have just grown to R_drift:

$$\begin{eqnarray*} r_{\mathrm{ g}}(t) = {\left(\frac{3}{16}\right)}^{1/3}{(G M_{*})}^{1/3}{(\epsilon _{\mathrm{ d}} Z_{0})}^{2/3}t^{2/3}, \end{eqnarray*}$$

(10)

where ϵ_d = 0.05 defines the size-dependent dust growth efficiency, and Z₀ = 0.01 is the solids-to-gas ratio, contributed by both pebbles (Z_peb) and dust (Z_dust):

$$\begin{eqnarray*} Z_{0} = Z_{\textrm {peb}} + Z_{\textrm {dust}}. \end{eqnarray*}$$

(11)

In this work, we assume that 90 per cent of the total solids are converted into pebbles, and this conversion rate remains constant throughout the simulation.

As the dust growth time-scale is shorter at smaller disc radii (i.e. dust at smaller radii grows faster), the pebble production front moves outwards from the disc inner edge over time [see Fig. 3 for the plot of r_g(t) for the disc configuration used in this paper]. This provides a pebble mass flux (flux of pebbles drifting inwards from r_g) defined as

$$\begin{eqnarray*} \dot{M}_{\textrm {flux}} = 2 \pi r_\mathrm{ g} \frac{\mathrm{ d}r_{\mathrm{ g}}}{\mathrm{ d}t} Z_{\textrm {peb}}(r_\mathrm{ g}) \Sigma _{\textrm {gas}}(r_\mathrm{ g}), \end{eqnarray*}$$

(12)

where Σ_gas(r_g) is the gas surface density at the pebble production front. From the mass flux, we define the pebble surface density profile Σ_peb the same way as Lambrechts & Johansen (2014):

$$\begin{eqnarray*} \Sigma _{\textrm {peb}} = \frac{\dot{M}_{\textrm {flux}}}{2 \pi r v_{r}}, \end{eqnarray*}$$

(13)

where r is the disc radius and v_r is the radial velocity of pebbles at r, defined as

$$\begin{eqnarray*} v_r = 2 \frac{{\, {\rm St}}}{{\, {\rm St}}^{2} + 1} \eta v_{k}-\frac{v_{\textrm {r,gas}}}{1+ {\, {\rm St}}^{2}} \end{eqnarray*}$$

(14)

(Weidenschilling 1977; Nakagawa, Sekiya & Hayashi 1986), where St is the Stokes number of the pebbles, v_k is the local Keplerian velocity, v_{r, gas} is the local gas radial velocity, and η is the dimensionless measure of gas pressure support:

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

(15)

where H is the disc scale height (Nakagawa et al. 1986). For the Stokes number, we assume it is equal to

$$\begin{eqnarray*} {\, {\rm St}}= \min ({\, {\rm St}}_{\rm drift}, {\, {\rm St}}_{\rm frag}) , \end{eqnarray*}$$

(16)

where |${\, {\rm St}}_{\rm drift}$| is the drift-limited Stokes number that is obtained through an equilibrium between the drift and growth of pebbles to fit constraints of observations of pebbles in protoplanetary discs and from advanced coagulation models (Birnstiel, Klahr & Ercolano 2012)

$$\begin{eqnarray*} {\, {\rm St}}_{\rm drift} = \dfrac{\sqrt{3}}{8}\dfrac{\epsilon _{\rm p}}{\eta }\dfrac{\Sigma _{\rm peb}}{\Sigma _{\rm gas}}, \end{eqnarray*}$$

(17)

where ϵ_p is the coagulation efficiency between pebbles that we assume is similar to the dust growth efficiency and is equal to 0.5. As well as the drift-limited Stokes number, we also include the fragmentation-limited Stokes number (⁠|${\, {\rm St}}_{\rm frag}$|⁠), which we follow (Ormel & Cuzzi 2007) and is equal to

$$\begin{eqnarray*} {\, {\rm St}}_{\rm frag} = \dfrac{v_{\rm frag}^2}{3\alpha c_{\rm s}^2} , \end{eqnarray*}$$

(18)

where v_frag is the impact velocity required for fragmentation, which we model as the smoothed function

$$\begin{eqnarray*} \dfrac{v_{\rm frag}}{1 \rm ms^{-1}} = 10^{0.5+0.5\tanh {((r-r_{\rm snow})/5H)}}, \end{eqnarray*}$$

(19)

where r_snow is the water snowline, which we assume to be where the disc mid-plane temperature is equal to 170 K. The fragmentation velocity therefore varies between 1 ms⁻¹ for rocky pebbles (Güttler et al. 2010) and 10 ms⁻¹ for icy pebbles, consistent with some results in the literature (though this is still an area of open research; Gundlach & Blum 2015; Musiolik & Wurm 2019).

For each simulation, we inject the planetary embryo core at its specified initial semimajor axis a_init as soon as r_g has reached a_init (i.e. we inject it at t when r_g = a_init) with the initial planetary embryo core mass M_core,init defined as 10 per cent of the transition mass M_trans for the core. The transition mass is a core mass threshold above which the pebble accretion on to the core transitions from the Bondi regime to the Hill regime (Johansen & Lambrechts 2017). The Bondi accretion regime happens initially for low embryo core mass, when the core accretes pebbles that pass through its Bondi radius R_B, which is smaller than its Hill radius R_H at low mass. As the core becomes more massive and reaches M_trans, its Bondi radius becomes similar to the Hill radius, and the accretion regime switches to the Hill regime where the accretion rate is limited by its Hill sphere. The transition mass is calculated as

$$\begin{eqnarray*} M_{\textrm {trans}} = \eta ^{3} M_{*}. \end{eqnarray*}$$

(20)

Coleman (2021) studied the sizes and distributions of planetary embryos and planetesimals formed in discs via pebble trapping in short-lived pressure bumps, and found that the masses of the planetary embryos formed are at least one order of magnitude lower than the transition mass, which is the reason we chose our initial embryo core masses as 0.1 M_trans.

As in Coleman (2021), we use two modes for core growth via pebble accretion: a 2D mode or a 3D mode, depending on the core’s Hill radius R_H and the scale height of the pebbles H_peb. If R_H < H_peb, we use the 3D accretion mode since there are pebbles crossing the core’s entire Hill sphere, but as the core mass increases and when R_H > H_peb, we switch to the 2D accretion mode as some regions of its Hill sphere are empty of pebbles. We use the 2D and 3D accretion rates in Johansen & Lambrechts (2017):

$$\begin{eqnarray*} \dot{M}_{\textrm {2D}} = 2 R_{\textrm {acc}} \Sigma _{\textrm {peb}} \delta v, \end{eqnarray*}$$

(21)

and

$$\begin{eqnarray*} \dot{M}_{\textrm {3D}} = \pi {R_{\textrm {acc}}}^{2} \rho _{\textrm {peb}} \delta v, \end{eqnarray*}$$

(22)

where ρ_peb is the mid-plane pebble density, and δv = Δv + ΩR_acc is the approach speed, with Δv as the gas sub-Keplerian velocity. R_acc is the accretion radius that depends on whether the accreting core is in the Hill or Bondi regime, and on the pebble friction time-scale. The dependence on friction time arises due to pebbles having to change directions significantly on time-scales shorter than the friction time-scale, which brings a criterion accretion radius |$\hat{R}_{\textrm {acc}}$| defined as

$$\begin{eqnarray*} \hat{R}_{\textrm {acc}} = {\left(\frac{4t_\mathrm{ f}}{t_\mathrm{ B}} \right)}^{1/2} R_\mathrm{ B } \end{eqnarray*}$$

(23)

for the Bondi regime, where |$t_\mathrm{ f} = {\, {\rm St}}/\Omega$| is the pebble friction time-scale, t_B is the Bondi sphere crossing time, and

$$\begin{eqnarray*} \hat{R}_{\textrm {acc}} = {\left(\frac{\Omega _{k} t_\mathrm{ f}}{0.1} \right)}^{1/2} R_\mathrm{ H } \end{eqnarray*}$$

(24)

for the Hill regime. The accretion radius is the defined as

$$\begin{eqnarray*} R_{\textrm {acc}} = \hat{R}_{\rm acc} \exp { \left[- \chi {(t_\mathrm{ f}/t_\mathrm{ p})}^{\gamma } \right] }, \end{eqnarray*}$$

(25)

where t_p = GM/(Δv + ΩR_H)³ is the characteristic passing time-scale with χ = 0.4 and γ = 0.65 (Ormel & Klahr 2010).

The planetary embryos grow by accreting pebbles until they reach the so-called pebble isolation mass, that is the mass required to perturb the gas pressure gradient in the disc: i.e. the gas velocity becomes super-Keplerian in a narrow ring outside of a planet’s orbit reversing the action of the gas drag. The pebbles are therefore pushed outwards rather than inwards and accumulate at the outer edge of this ring stopping the embryos from accreting solids (Paardekooper & Mellema 2006; Rice et al. 2006). Initial work found that the pebble isolation mass was proportional to the cube of the local gas aspect ratio (Lambrechts & Johansen 2014). More recent work, however, has examined what effects disc viscosity and the Stokes number of the pebbles have on the pebble isolation mass, finding that small pebbles that are well coupled to the gas are able to drift past the pressure bump exterior to the planet’s orbit (Ataiee et al. 2018; Bitsch et al. 2018). To account for the pebble isolation mass while including the effects of turbulence and stokes number, we follow Bitsch et al. (2018), and define a pebble isolation mass-to-star ratio,

$$\begin{eqnarray*} q_{\rm iso} = \left(q_{\rm iso}^{\dagger } + \frac{\Pi _{\rm crit}}{\lambda } \right) \frac{M_{\oplus }}{M_*} , \end{eqnarray*}$$

(26)

where |$q_{\rm iso}^\dagger = 25 f_{\rm fit}$|⁠, λ = 0.004 76/f_fit, |$\Pi _{\rm crit} = \frac{\alpha }{2\mathrm{ St}}$|⁠, and

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

(27)

with α₃ = 0.001. Once planetary embryos reach the pebble isolation mass, they no longer accrete pebbles from the discs in our simulations.

2.3 Gas envelope accretion

Once a planet reaches sufficient mass through pebble accretion, it is able to accrete a gaseous envelope from the surrounding disc. Ideally, we would incorporate 1D envelope structure models (e.g. Coleman, Papaloizou & Nelson 2017) into our simulations. However, these calculations are computationally expensive and would considerably increase the simulation run times; therefore, we opted to include gas accretion fits instead. Previous work (Coleman & Nelson 2016a) provided fits to the 1D calculations presented in Movshovitz et al. (2010), which was for a planet located at 5.2 |$\, {\rm au}$|⁠. More recently, Poon, Nelson & Coleman (2021) presented an updated gas accretion model based on fits to gas accretion rates obtained using a 1D envelope structure model (Papaloizou & Terquem 1999; Papaloizou & Nelson 2005; Coleman et al. 2017). To calculate these fits, Poon et al. (2021) performed numerous simulations, embedding planets with initial core masses between 2 and 15 |$\, {\rm M}_{\oplus }$| at orbital radii spanning 0.2–50 |$\, {\rm au}$|⁠, within gas discs of different masses. This allowed for the effects of varying local disc properties to be taken into account when calculating gas accretion fits, a significant improvement on fits from other works (e.g. Hellary & Nelson 2012; Coleman & Nelson 2016a). Using the 1D envelope structure model of Coleman et al. (2017), the embedded planets were then able to accrete gas from the surrounding gas disc until either the protoplanetary disc dispersed or the planets reached a critical state where they would then undergo runaway gas accretion. With the results of these growing planets, Poon et al. (2021) calculated fits to the gas accretion rates taking into account both the planet and the local disc properties.

Following the fits found in Poon et al. (2021), the gas accretion rate we use is equal to

$$\begin{eqnarray*} \left(\dfrac{\rm{d} M_{\mathrm{ge}}}{\rm{d} t} \right)_{\mathrm{local}}&=& 10^{-10.199} \left(\dfrac{\mathrm{M_{\oplus }}}{\mathrm{yr}}\right) f_{\mathrm{opa}}^{-0.963} \left(\dfrac{T_{\mathrm{local}}}{\mathrm{1\, K}}\right)^{-0.7049}\nonumber \\ &&\times \, \left(\dfrac{M_{\mathrm{core}}}{\mathrm{M}_{\oplus }} \right)^{5.6549} \left(\dfrac{M_{\mathrm{ge}}}{\mathrm{M}_{\oplus }} \right)^{-1.159}\nonumber \\ &&\times \, \left[\exp {\left(\dfrac{M_{\mathrm{ge}}}{M_{\mathrm{core}}} \right)} \right]^{3.6334}, \end{eqnarray*}$$

(28)

where T_local is the local disc temperature, f_opa is an opacity reduction factor (which reduces the grain opacity contribution, kept constant at f_opa = 0.01 in all simulations), and M_core and M_ge are the planet’s core and envelope masses, respectively. When comparing the masses of gas accreting planets calculated through equation (28) to the actual masses obtained using the 1D envelope structure model of Coleman et al. (2017), Poon et al. (2021) found excellent agreement, a considerable improvement on previous fits (e.g. Coleman & Nelson 2016a). In each of our simulations, we allow the planetary embryo to accrete a gaseous envelope once their mass exceeds an Earth mass, and all mass accreted by planets is removed from the disc to conserve mass.

2.4 Type I migration

In our simulations, planetary embryos that grow massive enough (significantly exceeding a Lunar mass) are subject to type I migration, which is implemented in the same way as Coleman & Nelson (2014), based on the torque formulae presented by Paardekooper et al. (2010) and Paardekooper, Baruteau & Kley (2011). These formulae specify how the different torques (that contribute to the total torque acting on the planetary embryo, expressed by equation 29) change with planet mass and local disc conditions (determined by changes in temperature, surface density, and metallicity/opacity). Corotation torques in particular vary sensitively with the ratio of the horseshoe libration time-scale to either the viscous or thermal diffusion time-scales.

The torque experienced by a low-mass planet embedded in a disc is contributed by the Lindblad torques and a weighted sum of the vorticity-related horseshoe drag, the entropy-related horseshoe drag, the vorticity-related linear corotation torque, and the entropy-related linear corotation torque. Combining equations (50)–(53) in Paardekooper et al. (2011), we express the total torque for type I migration Γ_tot acting on a planet as

$$\begin{eqnarray*} \Gamma _{\textrm {tot}} &=& \Gamma _{\mathrm{ LR}} + \lbrace \Gamma _{\textrm {VHS}} F_{\nu } G_{\nu } + \Gamma _{\textrm {EHS}} F_{\nu } F_{\mathrm{ d}} \sqrt{G_{\nu } G_{\mathrm{ d}}} \nonumber\\ &&+\, \Gamma _{\textrm {LVCT}}(1-K_{\nu }) +\Gamma _{\textrm {LECT}}(1-K_{\nu })(1-K_\mathrm{ d})\rbrace , \end{eqnarray*}$$

(29)

where Γ_LR is the Lindblad torque, Γ_VHS is the vorticity-related horseshoe drag, Γ_EHS is the entropy-related horseshoe drag, Γ_LVCT is the vorticity-related linear corotation torque, and Γ_LECT is the entropy-related linear corotation torque (as given by equations 3–7 in Paardekooper et al. 2011). The functions G_ν, G_d, F_ν, F_d, K_ν, and K_d are related to either the ratio between the viscous/thermal diffusion time-scales and the horseshoe liberation time-scale or the ratio of the viscous/thermal diffusion time-scales and the horseshoe U-turn time-scale (see equations 23, 30, and 31 in Paardekooper et al. 2011). As our simulations only consider circular coplanar orbits, the effects of a planet’s eccentricity and inclination on attenuating the Lindblad and corotation torques are not required, although these components are present in the code (see Coleman & Nelson 2014).¹

3 RESULTS AND DISCUSSION

We begin by studying the mass budget and duration of pebble flux through the disc as a function of external FUV field strength and shielding time-scale in Section 3.1. We then move on to discuss the impact on the planets forming in our models in Section 3.2. Finally, in Section 3.3 we discuss our understanding of the shielding time, how that links to the star formation process, and how that affects our ability to predict exoplanet populations.

3.1 Pebble cut-off time

Planet formation by pebble accretion relies on a flux of pebbles from the outer disc. Grains undergo radial drift once they reach sufficient size to dynamically decouple from the gas (Weidenschilling 1977). The time-scale for growth to such a size increases radially outwards through the disc (e.g. Birnstiel et al. 2012), resulting in a ‘pebble production front’ at radius r_g. Once the pebbles have grown to a size where they undergo radial drift, they are effectively safe from external photoevaporation, both because only small grains are entrained in a wind (Facchini et al. 2016) and because radial drift moves them in to smaller radii where external photoevaporation does not operate.

The pebble cut-off time, t_cut-off, is the time a disc has before its pebble production front r_g reaches the disc outer edge, cutting off the supply of pebbles. t_cut-off relative to the time that the pebble production front passed a planetary embryo’s location determines the time available for the embryo to grow via pebble accretion. If a disc is subject to external photoevaporation and the mass-loss rate is greater than the rate of viscous spreading, the disc will be truncated (Clarke 2007). Consequently, the disc outer edge moves inwards and meets the outward-moving pebble production front earlier than it otherwise would have. This reduces the pebble cut-off time compared to when the disc is not truncated. Therefore, how much the pebble cut-off time is shortened, i.e. how long the planetary embryos have to accrete pebbles, depends on the strength of the radiation field and the shielding time during which the disc is protected from truncation via external photoevaporation.

The dashed line in Fig. 2 shows the pebble production front r_g over time in our modelled disc. It compares the production front radius with the disc outer edge when subjected to a strong FUV irradiation field (F_FUV,max = 10⁴ G₀), for various shielding times t_sh. Note that here we define the disc outer edge (or disc radius) as the radius that contains 90 per cent of the gas mass, which we denote as R_d,90. It shows that as soon as the disc is exposed to the strong FUV field, the disc radius is rapidly truncated down to ∼40 au, and then undergoes a much slower phase of truncation until fully dispersed. A longer shielding time protects the disc from the rapid truncation phase and hence the pebble cut-off time, i.e. the time when r_g crosses R_d,90, gets longer as t_sh increases.

$The pebble production front rg as a function of time for the disc configuration simulated (dashed curve) and disc radii Rd,90 evolution under and FUV track with $F_{\textrm {FUV,~max}} =10^4\, \mathrm{ G}_0$ and various shielding times (solid lines). Longer shielding means that rg intersects the disc outer radius at a later time, resulting in a longer period of pebble flux through the disc.$

Figure 2.

The pebble production front r_g as a function of time for the disc configuration simulated (dashed curve) and disc radii R_d,90 evolution under and FUV track with |$F_{\textrm {FUV,~max}} =10^4\, \mathrm{ G}_0$| and various shielding times (solid lines). Longer shielding means that r_g intersects the disc outer radius at a later time, resulting in a longer period of pebble flux through the disc.

Open in new tab Download slide

The correlation of pebble cut-off time with shielding time depends on F_FUV,max, as illustrated in Fig. 3 that shows t_cut-off as a function of t_sh for discs evolved under various F_FUV,max. For the disc evolved with F_FUV,max = 10⁴ G₀ (the strongest FUV strength in Fig. 2), t_cut-off and t_sh are almost in a 1:1 relationship for t_sh < 2.5 Myr. This is due to the disc outer edge being immediately stripped to inside of r_g as soon as the shielding ends. By 2.5 Myr, r_g reaches the original disc outer edge, resulting in the maximum possible cut-off time, and so longer shielding beyond 2.5 Myr has no additional impact. For lower values of F_FUV,max in the range of t_sh < 2.5 Myr, t_cut-off initially stays almost constant for shorter lengths of t_sh (because weaker FUV radiation reduces the disc size less efficiently), before again becoming linearly correlated to t_sh. The weaker the F_FUV,max, the longer the t_cut-off remains constant before starting to increase with t_sh.

Figure 3.

The pebble cut-off time t_cut-off as a function of the shielding time t_sh (solid lines) for FUV tracks of different F_FUV,max values indicated by different colours. The dashed line highlights the 1:1 ratio as a reference.

Open in new tab Download slide

Overall, Fig. 3 shows that the cloud shielding time is important in determining t_cut-off over a wide range of UV field strengths. This affects the total mass budget in drifting pebbles and the time-scale over which there is a pebble flux through the disc. We next turn our attention to how this affects planet formation by pebble accretion.

3.2 Planet growth and migration

A planetary embryo begins accreting pebbles once the pebble production front reaches the orbital radius of the embryo (r_g = a). This continues until shortly after the pebble growth front reaches the disc outer edge (i.e. at t_cut-off). This in turn depends on the FUV radiation strength and shielding time, since those control the evolution of the disc outer radius.

Fig. 4 shows the evolution of planet masses M_p and semimajor axes a, for discs irradiated by a 10⁴ G₀ radiation field after some shielding time. The circles represent the initial mass/semimajor axis of the planetary embryos, and the arrowed tracks show the subsequent evolution. Each panel represents the evolution of sets of planetary embryos for different shielding time-scales, from no shielding in the upper left to 3 Myr of shielding in the lower right panel.

$Evolution of planetary embryo mass Mp against semimajor axis a in discs evolved with FUV tracks with FFUV,max = $10^4\, \mathrm{ G}_0$. In each subplot, simulations with one specific shielding time tsh are included, and each coloured curve represents the evolution trajectory of one planetary embryo in one disc with arrows indicating the direction of evolution. The initial planetary embryo core mass and semimajor axis (ainit) when injected are indicated by the dots. Note that the plots only include the simulations with ainit < 60 au though the full ainit range simulated is 10–100 au.$

Figure 4.

Evolution of planetary embryo mass M_p against semimajor axis a in discs evolved with FUV tracks with F_FUV,max = |$10^4\, \mathrm{ G}_0$|⁠. In each subplot, simulations with one specific shielding time t_sh are included, and each coloured curve represents the evolution trajectory of one planetary embryo in one disc with arrows indicating the direction of evolution. The initial planetary embryo core mass and semimajor axis (a_init) when injected are indicated by the dots. Note that the plots only include the simulations with a_init < 60 au though the full a_init range simulated is 10–100 au.

Open in new tab Download slide

This figure illustrates how the cloud shielding time affects the growth and migration of planetary embryos at different semimajor axes. First, the top left subplot with t_sh = 0 demonstrates the effects of strong external photoevaporation in suppressing the core growth by the fast truncation of the disc outer edge and hence quickly cutting off the supply of pebble flux (shortening t_cut-off) for planetary embryo cores injected at almost all semimajor axes. For r > 35 au, the disc is truncated so efficiently that the disc outer edge has already moved inside the pebble production front r_g, before r_g reaches 35 au, and so no embryo cores were injected beyond 35 au. Even within 35 au, there is hardly any growth and migration for cores injected at >20 au due to short t_cut-off, and only the innermost three cores (injected at 10, 11, and 13 au) had enough time to grow and accrete meagre gaseous envelopes but only reached masses between 1 and 10 M_⊕.

With such strong irradiation, even a short shielding time makes a difference in delaying rapid truncation and increases t_cut-off, allowing the cores more time to receive the pebble flux and grow. The top right subplot with t_sh = 0.3 Myr, for example, already shows a difference in the final masses that the cores were able to reach, especially for those with a_init between 14 and 25 au, whose final masses are at least one order of magnitude larger than the case with t_sh = 0. As t_sh increases, more cores injected at larger a_init are able to accrete for longer and so grow more massive and migrate. This behaviour is demonstrated further in the top plot of Fig. 5, which plots the final masses M_final that the planetary embryos injected at various a_init were able to reach as a function of t_sh.

$The final mass Mfinal (upper panel) and final semimajor axis afinal (lower panel) that the injected planetary embryos reached as a function of shielding time tsh in discs evolved with the FUV radiation track of $F_{\textrm {FUV,~max}} = 10^4\, \mathrm{ G}_{0}$, with different colours indicating the initial semimajor axis ainit of the planetary embryos.$

Figure 5.

The final mass M_final (upper panel) and final semimajor axis a_final (lower panel) that the injected planetary embryos reached as a function of shielding time t_sh in discs evolved with the FUV radiation track of |$F_{\textrm {FUV,~max}} = 10^4\, \mathrm{ G}_{0}$|⁠, with different colours indicating the initial semimajor axis a_init of the planetary embryos.

Open in new tab Download slide

The top plot of Fig. 5 shows that for a planetary embryo injected at some radial distance a_init, the M_final it is able to reach increases with shielding time initially until a certain value of t_sh, after which M_final remains roughly the same. This is due to the shielding time being long enough for the pebble production front to reach the disc outer edge before being exposed to external photoevaporation. In that case, the disc retains the maximal pebble mass budget and so extending the shielding time-scale further has no impact on planet formation.

The shielding time makes the most difference for the embryos injected between a_init = 14 and 31 au in terms of the steepness of the correlation between M_final and t_sh before plateauing. Embryos in this radial range reach masses of around a lunar mass (∼0.01 M_⊕) when not shielded, but were able to reach ∼5 M_⊕ when shielded for longer than ∼1.5 Myr. With long enough t_sh, all the embryos with a_init < = 31 au are able to reach at least 0.1 M_⊕, with those of a_init < = 25 au able to achieve |$1\!-\!10\, \mathrm{ M}_{\oplus }$|⁠. For the embryos injected with a_init > 31 au, as it takes longer for r_g to reach their location so that they can start accreting pebbles, the time left for their pebble accretion before r_g crosses disc outer edge is not long enough for them to grow larger than 0.1 M_⊕ even with relatively long t_sh. The opposite occurs for the embryos with a_init < 14 au; as r_g reaches their locations quite early on, they had enough time to grow beyond 1 M_⊕ and migrate inward rapidly even without shielding.

Figs 4 and 5 demonstrate that the shielding time makes a significant difference to the growth and migration of planets in 10⁴ G₀ environments. In Section 3.1, we showed that t_sh also affects t_cut-off in lower FUV radiation environments. To explore this, Fig. 6 shows how the final mass of an object resulting from an embryo injected at a_init = 25 au scales as a function of the FUV radiation environment once the disc is unshielded. The final planet mass scales non-linearly with t_sh up to around 1.5 Myr in the high UV field cases (10³–10⁴ G₀). Even in the case of a 100 G₀ radiation field, the difference in mass between no shielding and 1.5 Myr of shielding is an order of magnitude. Fig. 7 also shows a contour plot that summarizes how the final planet mass/semimajor axes vary as a function of the initial semimajor axis and shielding time-scale.

Figure 6.

The final mass M_final that the planetary embryos injected at a_init = 25 au reached as a function of shielding time t_sh in discs evolved with FUV tracks of different F_FUV,max values, indicated by different colours. There is a non-linear impact of shielding time on the planet’s final mass up to 1 Myr for F_FUV,max ≥ 100 G₀.

Open in new tab Download slide

The colour scale shows the final mass of a planet as a function of the shielding time and initial orbital radius of the planetary seed, in a 104 G0 environment. The contours show the final semimajor axis (afinal) of the planets in au.

Figure 7.

The colour scale shows the final mass of a planet as a function of the shielding time and initial orbital radius of the planetary seed, in a 10⁴ G₀ environment. The contours show the final semimajor axis (a_final) of the planets in au.

Open in new tab Download slide

Beyond a certain shielding time-scale (⁠|${\sim} 1.5\,$|Myr in this case), the mass reservoir for pebble accretion has reached the disc outer edge and further shielding from external photoevaporation has no additional impact on planet properties. Beyond such a time, other processes, such as planet migration and N-body interactions, will have a greater impact on the types of planets that form in the simulations.

3.3 The boundedness of clouds and shielding time-scales of discs

We have demonstrated that planet formation by pebble accretion, a widely used planet formation paradigm, is sensitive to external photoevaporation. Furthermore, it is sensitive to the time-scale that a disc can be shielded from external photoevaporation up to shielding time-scales of ∼1.5 Myr. We now turn our attention to the expected shielding time-scales of discs.

It is well established in numerical simulations that the extent and rapidity with which a star-forming cloud is dispersed by feedback is a function of the cloud boundedness (virial ratio; see, for example, Dale et al. 2012a, b, 2013a, b), whose simulations demonstrate that globally bound clouds are dispersed more slowly. Since the rate that the molecular gas is dispersed decreases for more bound clouds, we therefore expect that the shielding time-scale will also correlate with the boundedness of the cloud. Other recent studies of feedback in clouds find that protostellar jets and other forms of feedback can play a role in supporting against very low virial ratios, preventing it from dropping significantly below unity on a global scale of the cloud (Grudić et al. 2022; Guszejnov et al. 2022).

The expectation that shielding will be affected by cloud boundedness is supported by recent simulations that study the expected UV irradiation of discs in star formation simulations (Qiao et al. 2022; Wilhelm et al. 2023). These papers studied star formation from clouds with initial global virial ratios of α = 2 and 0.25, respectively, which effectively correspond to an upper limit on the virial ratio for collapse to occur and lower extreme of observed virial ratio, respectively. In accordance with our expectations from the above discussion, Qiao et al. (2022) and Wilhelm et al. (2023) find shielding time-scales of <0.5 Myr and of at least 0.5 Myr, respectively. There will be a continuum of behaviours between these two examples, which is yet to be explored in detail. It will be important to do so, since when coupled with our results, this range of shielding times could lead to either planet formation being dramatically influenced by the environment for a shielding time <0.5 Myr or it being a much more marginal (or even negligible) effect for shielding |${\gt} 1.5\,$|Myr.

We note that measured virial ratios of clouds vary broadly (see e.g. fig. 1 of Kauffmann, Pillai & Goldsmith 2013, where they vary from α ∼ 3 × 10⁻² to almost 10³). A large marginally bound cloud like that of Qiao et al. (2022) will end up with very subvirial clumps at the sites of massive star formation (see e.g. Camacho et al. 2020, for simulations following the virial ratio of substructures from larger scale turbulent clouds). We also note that the virial ratio is also easy to be underestimated through systematics, e.g. as highlighted by Singh et al. (2021) for Gould Belt clumps (see also Traficante et al. 2018).

The discovery of proplyds in the very young (0.2–1 Myr) star-forming region NGC 2024 demonstrates that external photoevaporation certainly does happen at a time that would influence planet formation by pebble accretion (according to our models; Haworth et al. 2021). At the other extreme, Ginski et al. (2022) observed a disc with signs of planet formation emerging from an irradiated cloud in a region that is |${\sim} \, 5$|Myr, raising the possibility of a very long shielding time-scale. Observationally then, it is clear that both very short and long shielding time-scales are happening, but the actual probability distribution for real systems across the wide range of star forming scenarios is yet to be established theoretically or observationally.

3.4 Implications for planet populations

We have demonstrated here that planet formation by pebble accretion can be affected by the radiation environment, and that the shielding time-scale can have a non-linear impact on the resulting planets. However, the uncertainty of the shielding time-scale distribution and the likely coupling of that distribution to the particular star formation scenario could provide a major source of uncertainty in planetary population synthesis models.

While we have shown that the shielding time-scales can have non-linear effects on the resulting planets, our results also show that the radiation environment itself also has a large effect. Indeed, from Fig. 6 when there is little shielding, there is a large difference in the final planet masses for those forming at 25|$\, {\rm au}$| in weak environments (e.g. 10 G₀) to those in strong environments (e.g. 10⁴ G₀). The impact on the planetary population is therefore likely to vary from cluster to cluster. In the strong UV environments of massive star-forming regions, such as NGC 2024, the ONC, and Carina (Winter et al. 2020a), this could indicate that the final planet population might be significantly different from the one that would form in a weak environment, e.g. Taurus, mainly due to the reduction in pebbles and solids that are able to be accreted, as shown earlier. The differences in the environment will also affect the scale of migration that planets can undergo, which will influence the final populations that can be observed, as seen both here and in Winter et al. (2022).

There could also be additional processes that are affected by the early truncation of the protoplanetary disc through the removal of shielding from the parent cloud. For example, a significant simplification of our model is that we do not consider the effect of substructures in the disc on the dust evolution. These can cause pressure bumps that control the distribution of dust in a disc. Here, we assumed that the cores of planets can only grow through the accretion of pebbles. However, recent work that examined the population of planetesimals and planetary embryos that arise through the gravitational collapse of pebbles trapped at specific locations in the disc (Lenz, Klahr & Birnstiel 2019; Coleman 2021) could also be affected by the efficient truncation of the disc. Since the formation of planetesimals, and subsequently planetary embryos, may be dependent on the pebble production front, it is reasonable to assume that the effects on the front discussed in this paper will also affect their formation as well. This will further reduce the amount of solids available for accretion by planets, while simultaneously reducing the number and/or mass of the initial planetary embryos. With fewer, less massive planets, this could result in large discrepancies in the populations from one cluster to another, depending purely on the amount of shielding, which we will explore in future work.

An additional effect that would affect planetesimal and embryo formation is that the removal of gas would lead to higher dust-to-gas ratios in the outer disc. This may lead to an increase in the planetesimal formation rate via the streaming instability (Carrera et al. 2017). This could address some of the issues around a reduced planetesimal population through the pebble production front reaching the disc outer edge.

4 SUMMARY AND CONCLUSIONS

In this paper, we explored how planet formation by pebble accretion is affected by external photoevaporation, and how the impact of external photoevaporation can be mitigated by shielding the disc by embedding it in a molecular cloud for some time. We use a detailed and well-established disc evolution and pebble accretion model coupled with the FRIED grid of mass-loss rates to implement external photoevaporation. We draw the following main conclusions from this work:

It is already well established that a growth front of pebbles propagates from the inner disc to the disc outer edge. Those pebbles then undergo radial drift inwards, providing a mass flux for planet formation. External photoevaporation can truncate the outer disc dramatically (and rapidly). We show that if that happens before the pebble growth front reaches the disc outer edge, it acts to reduce the mass reservoir of pebbles and the time-scale over which there will be a pebble flux throughout the disc. Because the pebbles undergo radial drift, this impact on the outer disc can in principle affect planet formation anywhere in the disc.
The shielding time-scale has a non-linear impact on the resulting planet mass for planetary embryos originating beyond about 14 au. Note that planets formed originally at these locations in the disc do migrate into short-period orbits such as those observed for many exoplanet systems. A |$10^{-3}\, \mathrm{ M}_{\oplus }$| planetary embryo at 25 au, for example, stays at 25 au with a lunar mass if the disc is immediately irradiated by a 10³ G₀ field, but grows and migrates to be approximately Earth-like in both mass and orbital radius if the disc is shielded for just 1 Myr. Even relatively brief periods of shielding (<0.5 Myr) can drastically affect the final planetary system. Conversely, shielding for time-scales >1 Myr essentially nullifies the impact of external photoevaporation on the resulting planets.
The distribution of shielding time-scales that all stars experience is not well understood, but does have some constraints. Feedback clears less bound clouds more rapidly, exposing discs to UV radiation on a shorter time-scale. Qiao et al. (2022) simulations of marginally bound clouds find that shielding is <0.5 Myr, which is a time-scale that is impactful for planet formation by pebble accretion according to our models here. Conversely, at an extremely low global virial ratio Wilhelm et al. (2023) find shielding times of at least 0.5 Myr. Observed cloud virial ratios are uncertain, and span a range of values (discussed more in Section 3.3). Furthermore, externally irradiated discs are observed at both very young (perhaps as young as 0.2 Myr in NGC 2024) and much older ages. We hence do expect that planet formation by pebble accretion can be affected by the environment in practice. However, how widespread this impact is remains to be determined, which provides significant uncertainty for planetary population synthesis models.

In summary, we demonstrated that external photoevaporation and cloud shielding can have potentially significant impacts on planets formed via pebble accretion. Here, we used a simplified model of one planet per disc and a parametrized FUV track with a constant field strength (F_FUV,max) after a shielding period. However, in a dynamically evolving star-forming region, the FUV field radiated upon a disc constantly varies with time (Qiao et al. 2022; Wilhelm et al. 2023). In order for a realistic study of to what extent the formation of planets is affected by the cluster they are born in, planet formation models should be coupled to the FUV tracks traced from stellar cluster and feedback simulations. Including multiple planets in one disc could result in higher mass cores via collisions, which could also act to enhance the pebble accretion rate. These more complicated models will be explored in future work.

ACKNOWLEDGEMENTS

TJH funded by a Royal Society Dorothy Hodgkin Fellowship. GALC was funded by the Leverhulme Trust through grant RPG-2018-418. This work was performed using the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure. This research utilized Queen Mary’s Apocrita HPC facility, supported by QMUL Research-IT (http://doi.org/10.5281/zenodo.438045). We thank the anonymous reviewer for their comments that significantly improved the clarity of the paper.

DATA AVAILABILITY

The simulation results are available upon request to the corresponding author.

Footnotes

Note that although our code includes the component for type II migration, none of the planets in our simulations grew massive enough to open a gap in the disc and switch to type II migration, according to the gap opening criterion in Crida, Morbidelli & Masset (2006).

References

Adams

F. C.

Hollenbach

Laughlin

Gorti

2004

ApJ

611

360

Month:	Total Views:
April 2023	68
May 2023	86
June 2023	30
July 2023	45
August 2023	41
September 2023	39
October 2023	41
November 2023	28
December 2023	26
January 2024	49
February 2024	33
March 2024	39
April 2024	100
May 2024	49
June 2024	21
July 2024	47
August 2024	45
September 2024	38
October 2024	29
November 2024	40
December 2024	29
January 2025	14
February 2025	20
March 2025	36
April 2025	29

Article Contents

Planet formation via pebble accretion in externally photoevaporating discs

ABSTRACT

1 INTRODUCTION

2 NUMERICAL METHOD

2.1 Disc evolution

2.1.1 External photoevaporation

2.1.2 Shielding time-scale prescription

2.2 Pebble accretion

2.3 Gas envelope accretion

2.4 Type I migration

3 RESULTS AND DISCUSSION

3.1 Pebble cut-off time

3.2 Planet growth and migration

3.3 The boundedness of clouds and shielding time-scales of discs

3.4 Implications for planet populations

4 SUMMARY AND CONCLUSIONS

ACKNOWLEDGEMENTS

DATA AVAILABILITY

Footnotes

References

APPENDIX A: EXAMPLES OF THE DISC SURFACE DENSITY EVOLUTION

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only