© The Authors 2024

1 Introduction

Protoplanetary discs, the dusty and gaseous birthplaces of planets, are not stable and static environments. During the few million years in which the growth from dust to planets takes place, many physical parameters change. Very young stars experience strong outbursts that heat the disc, and these outbursts are followed by a phase of cooling. In this paper, we investigate how these ubiquitous heating events in young stellar systems influence the particle growth by shifting the water ice line outwards to up to tens of astronomical units.

In the quiescent stage of a typical protoplanetary disc, the water ice line is situated at 1–3 au. This corresponds to approximately 170 K, depending on the evolutionary stage of the disc (Martin & Livio 2012; Bitsch et al. 2015). Beyond the ice line, water transitions from vapour to a solid form, allowing the solid density to double at this disc location (Hayashi 1981; Abod et al. 2019). We thus expect the water ice line to potentially be an important location for growth towards planetesimals and planets. Therefore, understanding how particle growth takes place around the ice line is important in order to map the general growth from dust to larger bodies in the disc. Most works have so far been focused on a quiescent disc with a static ice line (e.g. Stevenson & Lunine 1988; Cuzzi & Zahnle 2004; Ros & Johansen 2013; Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017). However, the shifting of the ice line due to stellar outbursts will effectively reset the distribution of icy particles, and therefore, the consequences need to be investigated.

Early works that specifically modelled vapour condensation in the context of planet formation include Stevenson & Lunine (1988), who investigated outwards vapour diffusion across the ice line and found significant growth of icy bodies there. However, this work ignored dynamical processes, such as inwards diffusion and radial transport of the resulting icy bodies, and they also assumed the instantaneous formation of new icy bodies without the need of a seed particle, followed by a rapid growth stage into planetesimals. Cuzzi & Zahnle (2004) modelled the inwards drift of particles followed by sublimation and re-condensation of vapour on immobile sink particles. These works showed that the process of sublimation and condensation of vapour at the ice line has the potential to facilitate planetesimal growth there.

During the past decade, the dynamics of both vapour and particles were taken into account in modelling condensation and sublimation at the ice line by for example Ros & Johansen (2013); Schoonenberg & Ormel (2017); Ida & Guillot (2016). These works showed dust growth to pebble sizes and pile-ups of particles due to sublimation and re-condensation of vapour at the ice line: As icy particles drift inwards past the ice line, they sublimate, and the vapour that is released can diffuse back across the ice line and re-condense onto particles there. It is nevertheless still debated whether the conditions for triggering planetesimal formation by the streaming instability, in terms of particle Stokes number and metallicity, can really be reached around the water ice line (Estrada & Umurhan 2023).

The outcome of condensation at the ice line is strongly affected by the number of particles onto which vapour can condense. A large particle number leads to vapour being distributed on a large surface area and thus a small net growth per particle, whereas fewer particles means that each individual particle experiences a larger growth. However, not the total number of solid particles around the ice line alone matters: The fraction of these particles that have an icy mantle is even more important. It is known from atmospheric chemistry that the condensation of vapour onto bare rocky particles, so-called heterogeneous nucleation, is more difficult than condensation onto particles that are already covered by ice, so-called deposition (Vehkamäki 2006). The difference between heterogeneous nucleation and deposition was also confirmed in laboratory experiments (see Hoose & Möhler 2012, and references therein). This bottleneck for heterogeneous nucleation was not taken into account in some previous dust evolution models that included FU Orionis outbursts (Schoonenberg et al. 2017; Houge & Krijt 2023). The presence of particles prone to heterogeneous nucleation and condensation growth was first suggested to be important by Hubbard (2017). Assuming a set difference in the saturation ratio, that is, in the ratio of vapour pressure and saturated vapour pressure needed for condensation, between icy versus non-icy particles, the author found that this could lead to the growth of these icy particles to pebble sizes after the FU Orionis outburst at the expense of the non-icy particles, which stayed small.

In Ros et al. (2019), we modelled nucleation and deposition onto particles at a static water ice line using experimental temperature-dependent data for the saturation ratio required for condensation (Iraci et al. 2010). We found that the differentiation between icy and non-icy particles is indeed important in order to understand particle growth at the ice line. Growth takes place predominantly for already ice-covered particles because the vapour pressure outside the ice line is too low for heterogeneous nucleation (i.e. the formation of the first ice layer on an existing substrate) on bare silicate grains. This therefore leads to a dichotomy in which icy particles can grow to centimetre-sized pebbles in a narrow region around the ice line, while silicate particles stay dust-sized and diffuse out over the disc.

In this paper, we extend the work of Ros et al. (2019) to the cooling phase after the disc has experienced a thermal heating event. The FU Orionis cooling timescale of 100–1000 yr is short enough for the dynamics of vapour and pebbles to be negligible (Schoonenberg et al. 2017), but it is long enough for heterogeneous nucleation to become a rare event that only benefits a small subset of the silicate particles. In this way, we demonstrate the formation of centimetre-sized particles outside of the water ice line. These particles are large enough to trigger planetesimal formation by the streaming instability. An advantage to forming planetesimals in this way is that the pebble growth does not depend on the relatively poorly constrained diffusion rate of vapour across the ice line, in contrast to the prevalent model for planetesimal formation at the water ice line (Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017; Ros et al. 2019; Estrada & Umurhan 2023).

Our paper is organised as follows. In Sect. 2 we describe how we calculated the heterogeneous nucleation rate of ice on silicate dust, and we present our model for a protoplanetary disc that undergoes an FU Orionis outburst. We present the results of our calculation in Sec. 3. We discuss here that the condensation process is so fast that we can ignore physical effects such as collisions and sedimentation during the growth. In Sec. 4 we demonstrate that the ice particles grow large enough to potentially trigger planetesimal formation by the streaming instability. We briefly summarise our results and conclude in Sec. 5.

2 Method

2.1 Nucleation calculations

We discuss here the concept of heterogeneous nucleation, which means the formation of the first ice layer on a silicate particle, and how to calculate the nucleation rate. Nucleation onto silicate particles requires a significantly higher saturation than deposition (i.e. vapour-to-ice condensation) on grains that have already acquired an ice mantle. This is a key concept. A low nucleation rate therefore results in only a few but large icy particles, whereas a high nucleation rate allows many particles to grow, but only to comparatively smaller sizes. In order to understand the resulting size distribution, it is thus crucial to be able to model the nucleation rate correctly.

We calculated the nucleation rate using classical nucleation theory (CNT), which agrees well with experimental results when the typical prefactors are adopted (Hoose & Möhler 2012; Barahona 2012). In the following brief derivation, we differentiate between homogeneous nucleation, where water vapour nucleates into solid ice without a seed particle, and heterogeneous nucleation, where nucleation takes place on an existing solid seed particle. For water vapour at the ice line, only heterogeneous nucleation is relevant, but for simplicity, we here start from homogeneous nucleation in order to arrive at the expression for heterogeneous nucleation. For a more in-depth discussion of the various forms of nucleation, we refer to Hoose & Möhler (2012).

The energy barrier for the phase transition between vapour and ice is expressed in terms of the change in Gibbs free energy by formation of an ice particle of size R, $$\mathrm{\Delta }G=\frac{4\pi }{3}{R}^3\mathrm{\Delta }{g}_{\mathrm{v}}+4\pi R^2\sigma ,$$(1)

where σ is the surface tension, and Δg_v is the change in free energy per volume for the formation of the solid phase.

The critical cluster size beyond which an increase in R leads to a reduction in free energy is found from setting δΔG/δR = 0 to obtain $$R_{\mathrm{c}}=\frac{2\sigma }{-\mathrm{\Delta }{g}_{\mathrm{v}}}.$$(2)

The free energy difference can be expressed through the supersaturation factor S = P_v/P_sat, where P_v is the vapour pressure, and P_sat is the saturated vapour pressure, $$\mathrm{\Delta }{g}_{\mathrm{v}}=\left({\rho }_{•}/\mu \right)k_{\mathrm{B}}T\ln S,$$(3)

where ρ_•. is the solid density, μ is the molecular mass of the water molecule, k_B is the Boltzmann constant, and T is the temperature. The vapour pressure can be expressed using the ideal gas law, $$P_{\mathrm{v}}=\frac{k_{\mathrm{B}}T}{m_{\mathrm{v}}}{\rho }_{\mathrm{v}},$$(4)

where m_v is the mass of a water molecule, and ρ_v is the density of vapour. The saturated vapour pressure is given by the Clausius-Clapeyron equation, and it yields $$P_{\text{sat }}=6.034\times {10}^{11}\mathrm{P}\mathrm{a}\times {\mathrm{e}}^{-5938\mathrm{K}/T}$$(5)

using experimentally determined material constants for water (Haynes et al. 1992). The energy barrier for homogeneous nucleation to reach the critical size is thus $$\mathrm{\Delta }{G}_{\mathrm{h}\mathrm{o}\mathrm{m}}=\frac{16\pi {\sigma }^3}{3{\left(\mathrm{\Delta }{g}_{\mathrm{v}}\right)}^2}=\frac{16\pi {\left(\mu /{\rho }_{•}\right)}^2{\sigma }^3}{3{(k_{\mathrm{B}}T\ln S)}^2}.$$(6)

The rate at which clusters of the critical size form (expressed as per cubic metre per second) is $$J=J_0{\mathrm{e}}^{-\mathrm{\Delta }{G}_{\mathrm{h}\mathrm{o}\mathrm{m}}/\left(k_{\mathrm{B}}T\right)},$$(7)

where J₀ is a kinetic pre-factor, k_B is the Boltzmann constant, and T is the temperature (Kozasa & Hasegawa 1987). This is the homogeneous nucleation rate, which gives the number of individual nucleation events per volume and the time for the case of the formation of new ice particles without a seed. In Johansen & Dorn (2022), we discussed the various forms that have been suggested for the pre-factor J₀.

For heterogeneous nucleation, we followed Hoose & Möhler (2012) and reformulated the change in Gibbs free energy as $$\mathrm{\Delta }{G}_{\text{het }}=\mathrm{\Delta }{G}_{\text{hom}}f\text{, }$$(8)

where the parameter f(θ) describes the wetting of the surface through $$f(\theta )=\frac{(2+\cos \theta )(1-\cos \theta {)}^2}{4}\text{. }$$(9)

Here, θ is the contact angle, which describes the flattening of the nucleating particle due to its coupling with the surface energy of the substrate.

We can now generalise Eq. (7) to the heterogeneous nucleation rate through $$j=j_0{\mathrm{e}}^{-\mathrm{\Delta }{G}_{\text{het }}/\left(k_{\mathrm{B}}T\right)}.$$(10)

Here, j defines the number of nucleation events per particle surface area per second. For the heterogeneous kinetic pre-factor of Eq. (10), we used $$j_0=\frac{{\alpha }_{\mathrm{c}}}{\sqrt{f}}\frac{P_{\mathrm{v}}^2(\mu /\rho )}{\mu k_{\mathrm{B}}T{v}_{\mathrm{s}}}\sqrt{\frac{\sigma }{k_{\mathrm{B}}T}}\exp \left(\frac{\mathrm{\Delta }{G}_{\mathrm{d}}}{k_{\mathrm{B}}T}\right).$$(11)

Here, α_c is the sticking probability (assumed to be unity), P_v = n_vk_BT is the vapour pressure, ν_s is the lattice vibration frequency of ice, and ΔG_d is the desorption energy per molecule (Barahona 2012). The numerical values for the physical parameters needed to evaluate Eq. (11) are given in Barahona (2012). Through Eq. (10), we calculated the nucleation rate for a given temperature T and super-saturation level S.

The contact angle θ, introduced in Eq. (9), governs the ease with which nucleation takes place. For θ = π (no flattening), we recover homogeneous nucleation (f = 1), where vapour does not benefit from nucleating on existing seed particles, but instead forms new ice particles. The other limit case where cos θ = 1, giving f = 0, corresponds to a scenario without any barrier to nucleation. Practically, this is equivalent to neglecting the bottleneck for heterogeneous nucleation, as was assumed in many earlier studies of vapour deposition in an astronomical context. Neglecting this bottleneck results in a small overall growth of each particle. We calculated the growth in radius, ΔR, that each particle undergoes in a scenario where we compare the total particle mass including ice (m_ice+rock) to the total rock mass alone (m_rock) through $$\frac{m_{\text{ice }+\text{ rock }}}{m_{\text{rock }}}=\frac{\sum _{\mathrm{j}=1}^{\mathrm{n}}\frac{4\pi }{3}(R_{\mathrm{j}}^3{\rho }_{\text{rock }}+({(R_{\mathrm{j}}+\mathrm{\Delta }R)}^3-R_{\mathrm{j}}^3){\rho }_{\text{ice }})}{\sum _{\mathrm{j}=1}^{\mathrm{n}}\frac{4\pi }{3}{R}_{\mathrm{j}}^3{\rho }_{\text{rock }}}.$$(12)

Here, n is the number of particles, and ρ_rock and ρ_ice are the material densities of rock and ice, respectively. For m_ice+rock/m_rock = 2, and using the particle size distribution described in Sec. 2.3, we find that when heterogeneous nucleation is neglected, each particle grows equally in terms of radius by a very small amount, ΔR = 5.8 μm.

In our simulations with heterogeneous nucleation, we used a contact angle of cos θ = 0.83, corresponding to f = 0.0085, which is relevant for nucleation on silicate dust at a low temperature (Iraci et al. 2010). We performed additional experiments with cos θ = 0.9, corresponding to f = 0.0027, and cos θ = 0.95, corresponding to f = 0.00066. The first case simply reduced the region of large pebbles outside of the ice line, while the second case substantially suppressed the overall particle growth as the rate of heterogeneous nucleation increased dramatically. This agrees with the findings from the limit case where f = 0.

2.2 Disc model and cooling

We implemented a simple 1D model of a very young protoplanetary disc that accretes gas onto the star at a rate of $\dot{M}={10}^{-7}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$. We assumed that the accretion rate is steady throughout the inner regions of the protoplanetary disc. This allowed us to calculate the gas column density ∼ from the flux expression $\dot{M}=3\pi v\mathrm{\Sigma }$. Here, ν = αc_sH is the viscosity of the disc, c_s denotes the sound speed and H = c_s/Ω the gas scale height (Ω is the Keplerian frequency at distance r). We followed Ida et al. (2016) and set the disc temperature to $$T=150\mathrm{K}{\left(\frac{L_{\star }}{L_{\odot }}\right)}^{2/7}{\left(\frac{r}{\mathrm{a}\mathrm{u}}\right)}^{-3/7}.$$(13)

Here, L_* is the luminosity of the star. We ignored viscous heating here because this is likely inefficient if the angular momentum transport is dominated by disc winds (Mori et al. 2021).

During FU Orionis outbursts, young stars reach luminosities of up to 10²−10³ L⊙ (Hartmann & Kenyon 1996). These outbursts likely originate in the innermost regions of the protoplanetary disc (Zhu et al. 2007, 2009; Hillenbrand & Findeisen 2015; Pignatale et al. 2018) and decay on timescales of decades to centuries. We focused on large outbursts, with L = 10³ L_⊙, and on two values of the cooling timescale (100 and 1000 yr). The protoplanetary disc model was initialised in its high-luminosity state and was subsequently exponentially cooled down to its pre-outburst equilibrium temperature.

We define the nucleation timescale for a particle as $$t_{\text{nuc }}=\frac{1}{j{A}_{\text{sur }}},$$(14)

where A_sur is the surface area of the seed particle, and j is given by Eq. (10). This shows that the nucleation timescale scales with particle size as t_nuc ∝ R⁻².

To illustrate the rate of heterogeneous nucleation, we show in Fig. 1 the nucleation timescale for a particle of size R = 100 μm over a range of super-saturations and temperatures, calculated for the conditions at the water ice line. The figure clearly shows that the nucleation rate is a very steep factor of supersaturation. For all temperatures shown here, the nucleation timescale is extremely long for S ≤ 2. This means that nucleation effectively does not occur until S ≥ 2, at which point the cooling path of a particle crosses a steep nucleation cliff.

The cooling path of a particle at typical conditions at the water ice line is schematically illustrated with a green line. The cooling path is the same regardless of the cooling timescale, but the outcome in terms of the resulting number of ice-nucleated particles and their size changes depends on how fast the disc cools. When a disc cools quickly, the nucleation cliff is also climbed quickly, meaning that many particles are able to nucleate. The resulting particle sizes are therefore smaller. In the opposite case, where the cooling is slow, fewer particles nucleate and grow to larger sizes because the particles climb the nucleation cliff more slowly.

Fig. 1 Nucleation timescale for a single 100-micrometre-sized grain for a range of super-saturations (S) and temperatures (T). The dotted line marks the 1 Myr nucleation timescale, and the dashed line shows the 1 Gyr nucleation timescale. We illustrate a schematic cooling path at the conditions of the water ice line with the green line.

2.3 Particle growth

In Johansen & Dorn (2022), we calculated the homogeneous nucleation of refractory minerals in the innermost regions of the protoplanetary disc after an FU Orionis outburst. Here, we focus instead on heterogeneous nucleation of water ice farther away from the star. Heterogeneous nucleation has the added complexity that only a very small fraction of the silicate particles nucleate ice. We therefore followed two populations of particles: (i) the silicate particles that may experience heterogeneous nucleation, and (ii) the rare ice-nucleated particles that grow directly by vapour deposition.

We assumed that the disc is populated by silicate particles with a minimum size of 1 μm and a maximum size of 1 mm. We ignored smaller dust sizes because their high curvature will likely prevent the formation of a stable ice layer (Kuroiwa & Sirono 2011). The maximum size is motivated by models of dust coagulation (Güttler et al. 2010; Zsom et al. 2010) and by observational constraints on pebble sizes in protoplanetary discs (Carrasco-González et al. 2019). Any solid ice in the modelled part of the disc was assumed to have sublimated during the outburst and to have left behind water vapour and either intact or disrupted silicate aggregates (Spadaccia et al. 2022; Aumatell & Wurm 2011). The dust number N as a function of size is distributed as a power law dN/dR ∝ R^−q with index q = 3.5. These particles, as well as the water vapour, were assumed to be mixed with the gas and not sedimented to the mid-plane. We discretised the silicate size range into 12 logarithmically spaced super-particles populating each of the 400 radial annuli that span from 0.4 to 30 au. During the cooling of the disc, we kept track of the temperature T and super-saturation level S in each annulus. In every time step, we then calculated the number of physical particles in each super-particle, $$\mathrm{\Delta }{N}_i^{(\text{ice })}$$, that underwent heterogeneous nucleation.

Equation (10) gives us the rate of nucleation per unit area. This yields $\mathrm{\Delta }{N}_i^{(\mathrm{i}\mathrm{c}\mathrm{e})}=j4\pi R_i^2{N}_i\mathrm{\Delta }t$ new ice particles, with R_i denoting the radius of physical particles in a super-particle of index i, ΔN_i denoting the total number of physical particles in the super-particle in each time step, and Δt denoting the time step. These ice-covered particles were transferred to a new super-particle whose physical particles, initialised at the critical cluster size, were allowed to undergo direct depositional growth. In this way, the number of super-particles increased until the vapour reservoir was depleted.

After a particle has nucleated, it will grow purely by deposition of vapour. The growth rate by deposition for a spherical particle is $$\dot{R}=\frac{v_{\perp }}{{\rho }_{•}}({\rho }_{\mathrm{v}}-{\rho }_s).$$(15)

Here, ρ_•. is the material density of water ice, ρ_v is the water vapour mass density, and ρ_s is the saturated vapour density of water at the environmental temperature. We did not include dust smaller than 1 micron, and we therefore ignored that ρ_s increases for very small dust sizes (Kuroiwa & Sirono 2011). The water vapour flux speed is $$v_{\perp }=\sqrt{\frac{k_{\mathrm{B}}T}{2\pi m_{\mathrm{v}}}},$$(16)

where k_B is the Boltzmann constant, T is the temperature, and m_v is the water molecule mass.

3 Results

3.1 Pebble growth after outburst

In Fig. 2 we compare the mean particle sizes after an FU Orionis outburst for two different decay times, 100 and 1000 yr. For both decay times, growth to the largest sizes is achieved just outside of the quiescent ice line, where the cooling stops at a super-saturation just above unity. However, the cooling timescale does have a significant impact on the resulting particle sizes. For a fast cooling timescale, t_cool = 100 yr, a large number of grains nucleates. The vapour is therefore spread out over a large number of particles, resulting in a smaller growth for each individual particle. For a slow cooling timescale, t_cool = 1000 yr, we observe fewer nucleation events. The icy particles that do form by nucleation can grow to much larger sizes in this case.

The size distribution of icy particles is shown in Fig. 3 for representative locations in the disc. At each location (1.8, 2.1, and 3.6 au), the mass is concentrated in a narrow size range with particle sizes of about 10 cm just outside the quiescent ice line at 1.8 au, 1 cm at 2.1 au, and 1 mm at 3.6 au for the short cooling time of t_cool = 100 yr. For a longer cooling time of t_cool = 1000 yr, the particle sizes are even larger. The figure also shows the initial size distribution of silicate particles, as a comparison.

The evolution of a representative particle nucleating close to the ice line is shown in Fig. 4. Particles predominantly nucleate during the short time when supersaturation (shown in panel a) is about S = 2. Nucleation can be seen as an instantaneous event. As shown in panels b and c, the continued growth to pebbles with St ≈ 0.01 by vapour deposition also occurs quickly, on a timescale of about ten years. We can therefore assume that nucleation and depositional growth in this early growth phase take place without being influenced by other processes such as sedimentation and collisions. We examine this assumption more closely in the following sections.

Fig. 2 Mean sizes and Stokes number for particles after an FU Orionis outburst with a decay time of 100 yr in blue and 1000 yr in red. The sizes of silicate particles at the start of the simulation are shown as dotted lines, with a maximum size of one millimetre. Growth to the largest particles occurs just outside of the quiescent ice line.

Fig. 3 Mass distribution of particles (normalised such that the integral over the logarithmic radius is always unity) just outside of the ice line at 1.8 (yellow line), 2.1 (red line), and 3.6 au (blue line) after the ice line has returned to its quiescent state in a scenario with a cooling time of tcool = 100 yr. We indicate the particle size for the longer cooling time of 1000 yr with circles. The black line shows the mass distribution of bare silicate particles.

3.2 Particle collisions

In our model, we ignored collisions that could take place during the period of condensation growth. The condensed water molecules form layers of solid ice on the ice-nucleated particles. Macroscopic particles of solid ice, as opposed to aggregates made of microscopic monomers, are not expected to stick to each other at the relative speeds experienced in the protoplanetary disc (see Johansen et al. 2014, and references therein). However, solid ice particles may experience cracking and destruction in relatively low-speed collisions.

In Fig. 5 we plot the collision velocity as a function of particle size for three different turbulence levels, δ = 10⁻³, δ = 10⁻⁴, and δ = 10⁻⁵ (δ denotes the dimensionless turbulent diffusion coefficient, defined equivalently to the α-viscosity coefficient as the turbulent diffusion coefficient normalised by c_sH). We calculated the collision speed from the expression valid for two equal-sized particles, $$v_{\mathrm{c}}=\sqrt{3}\sqrt{\mathrm{S}\mathrm{t}}\sqrt{\delta }{c}_{\mathrm{s}}$$(17)

(Ormel & Cuzzi 2007). Here, c_s is the sounds speed of the gas, and St is the dimensionless Stokes number, which quantifies the strength of coupling between gas and dust and is defined as $$\mathrm{S}\mathrm{t}=\frac{\pi }{2}\frac{R{\rho }_{•}}{\mathrm{\Sigma }}$$(18)

in the Epstein drag regime (Weidenschilling 1977). Here, ∼ denotes the gas surface density. From Eq. (18), we can describe how the Stokes number scales with particle size, $$\text{ St }\sim 0.0015\left(\frac{R}{\mathrm{c}\mathrm{m}}\right){\left(\frac{r}{2\mathrm{a}\mathrm{u}}\right)}^{-15/14},$$(19)

where r is the distance from the central star. We assumed a gas accretion rate ${\dot{M}}_{\star }=3\pi \mathrm{\Sigma }\alpha c_{\mathrm{s}}H={10}^{-7}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$ when calculating the gas surface density ∼.

We also show in Fig. 5 the critical fragmentation velocity, where icy particles start to fragment in collisions (Hauk et al. 2015). A comparison of the collision speeds to the critical fragmentation velocity shows that icy particles reach at least centimetre-sizes without being hindered by fragmentation for all turbulence levels. Any growth beyond these sizes could nevertheless be impacted by fragmenting collisions. To understand the potential impact of collisions, we now compare the condensation growth timescale to the collision timescale.

The collision timescale can be written as τ_col = 1/(nσυ_col), where n is the number density, and σ is the collision cross section. Inserting the collision speed from Eq. (17), we arrive at the expression $${\tau }_{\mathrm{c}\mathrm{o}\mathrm{l}}=\left(\frac{300}{\mathit{\Omega }}\right){\left(\frac{\mathrm{S}\mathrm{t}}{0.01}\right)}^{0.5}{\left(\frac{{\rho }_{\mathrm{p}}/{\rho }_{\mathrm{g}}}{0.01}\right)}^{-1}{\left(\frac{\delta }{{10}^{-4}}\right)}^{-1/2},$$(20)

where ρ_p and ρ_g are the mid-plane dust and gas density. For centimetre-sized particles with St = 0.001, we then obtain a collision timescale of τ_col ≈ 40 yr and a timescale of τ_col ≈ 135 yr for decimetre-sized particles with St = 0.01.

We compared this to the growth timescale of a particle, defined as $${\tau }_{\mathrm{g}\mathrm{r}\mathrm{o}}=\frac{R}{\dot{R}},$$(21)

where the growth rate by deposition for spherical particles is given by Eq. (15).

The saturated vapour pressure is ρ_s ≪ ρ_v during nucleation because ice nucleation occurs when the supersaturation is S > 2. We therefore write the condensation rate as $$\dot{R}=3.5\mathrm{m}\mathrm{m}{\mathrm{y}\mathrm{r}}^{-1}\left(\frac{{\rho }_{\mathrm{v}}}{{10}^{-9}\mathrm{k}\mathrm{g}{\mathrm{m}}^{-3}}\right){\left(\frac{T}{170\mathrm{K}}\right)}^{1/2}.$$(22)

Here, ρ_v = 10⁻⁹ kg m⁻³ represents the characteristic vapour density at 2 au from the distance from the star and an accretion rate of $\dot{M}={10}^{-7}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$. This estimate corresponds well to the actual deposition rate reported in panel c of Fig. 4. The growth timescale follows from $R/\dot{R}$ as $${\tau }_{\mathrm{g}\mathrm{r}\mathrm{o}}\approx 0.3\mathrm{y}\mathrm{r}\left(\frac{R}{1\mathrm{m}\mathrm{m}}\right){\left(\frac{{\rho }_{\mathrm{v}}}{{10}^{-9}\mathrm{k}\mathrm{g}{\mathrm{m}}^{-3}}\right)}^{-1}{\left(\frac{T}{170\mathrm{K}}\right)}^{-1/2}.$$(23)

The collision timescale is thus several times longer than the growth timescale for millimetre-sized and centimetre-sized particles. This comparison thus shows that the frequency with which these collisions can occur during the growth phase is low, and it would therefore not impact our results significantly. However, collisions will be important during the subsequent sedimentation phase.

Fig. 4 Evolution of a representative pebble at 1.8 au, nucleating when S = 2, corresponding to a time of t ≈ 840 yr. Panel a shows the super-saturation S at 1.8 au. Panel b shows the particle size in metres and Stokes number. We mark St = 0.01 with a dotted line. Panel c shows the growth rate in millimetres per year, and panel d shows the sedimentation and growth timescales. The cooling time of the disc is tcool =100 yr.

Fig. 5 Ice particle collision speed as a function of the mean ice particle size in each radial bin for a turbulent diffusion coefficient of δ = 10−3, δ = 10−4, and δ = 10−5. The critical velocity for ice fragmentation is over-plotted as a dotted line. Ice particles can grow to a few centimetres before they risk fragmentation in collisions.

3.3 Sedimentation

We compared the growth timescale to the sedimentation timescale to show that our assumption that the particles are located in the full gas scale-height layer rather than in a dense sedimented mid-plane layer is valid for the duration of their growth process.

The timescale to sediment to the mid-plane of the protoplanetary disc is $${\tau }_{\text{sed }}\sim \frac{1}{\mathit{\Omega }\mathrm{S}\mathrm{t}},$$(24)

(Youdin & Lithwick 2007) where Ω is the orbital frequency. This gives $${\tau }_{\text{sed }}=450\mathrm{y}\mathrm{r}{\left(\frac{\mathrm{S}\mathrm{t}}{0.001}\right)}^{-1}{\left(\frac{r}{2\mathrm{a}\mathrm{u}}\right)}^{3/2}.$$(25)

Comparing this to the growth timescale in Eq. (23), we find that with a difference of two orders of magnitude at 1 cm (St = 0.001) and a similar magnitude at 10 cm (St = 0.01), we can conclude that the sedimentation timescale is too long for sedimentation to be important for growth by deposition. We therefore did not include sedimentation of particles in our simulations. However, settling will become important shortly after the deposition growth phase.

4 Implications for planetesimal formation

An important question for understanding the role of heterogeneous ice nucleation for planet formation is whether and how icy pebbles that formed by nucleation and deposition of vapour can grow further into planetesimals after they have sedimented to the mid-plane. Here, we place our results in the context of the leading mechanism of planetesimal formation, the streaming instability, in which particles clump together due to gas drag felt by the solids to form clumps that are dense enough for them to gravitationally collapse into solid planetesimals (Youdin & Goodman 2005; Johansen et al. 2007).

In Fig. 6 we plot the threshold for clumping via the streaming instability for a range of particle Stokes numbers and metallicities using the fit to computer simulations from Li & Youdin (2021). A combination of a high enough metallicity and particle Stokes numbers gives the parameter space within which clumping via streaming instability is possible. This area is the upper region in the figure, and the lower region corresponds to the parameter space within which clumping via streaming instability is not possible.

For a disc with with solar metallicity (Z = 0.005), where we only allowed water ice to contribute to the clumping, the particles are not large enough to trigger the streaming instability in our simulations. However, for the scenario in which the disc, either globally or locally, has a higher metallicity than Z = 0.01 or where both rock and ice are included as solids that can contribute to clumping (e.g. if the silicate dust aggregates collide with and adhere to the ice particles) the outcome is more promising. In this case, the largest particles that can reach centimetre-sizes just outside of the ice line fall in the range in which clumping via streaming instability is possible. This would also change the chemical composition of the planetesimals from pure ice (if formed only out of the ice particles that grew by condensation) to a mixture of silicates and ice.

We emphasise that we did not model the step from particle growth through sedimentation and to clumping by the streaming instability, but rather compared the particle sizes resulting from our model to typical theoretical limits for planetesimal formation through the streaming instability (Li & Youdin 2021; Lim et al. 2023).

Clumping by streaming instability is expected to take place in the mid-plane. The ratio of the solids to gas density here can be calculated for a specific particle size from the overall metallicity, Z, and the dimensionless turbulent diffusion coefficient, δ, as $$\begin{array}{rl}{\left({\rho }_{\mathrm{p}}/{\rho }_{\mathrm{g}}\right)}_{\text{mid }}& =Z\sqrt{\frac{\mathrm{S}\mathrm{t}+\delta }{\delta }}\sim Z\sqrt{\frac{\mathrm{S}\mathrm{t}}{\delta }}\\ & =0.3\left(\frac{Z}{0.01}\right){\left(\frac{\mathrm{S}\mathrm{t}}{0.01}\right)}^{0.5}{\left(\frac{\delta }{{10}^{-5}}\right)}^{-0.5}\end{array}$$(26)

(Johansen et al. 2014). For a disc with a low turbulence of δ = 10⁻⁵ and solar metallicity in which both rock and ice contribute to clumping, we find a mid-plane ratio of solids to gas of 0.3, corresponding to the value in the simulations by Li & Youdin (2021). Newer simulations by Lim et al. (2023) explored the threshold for streaming instabilities for a range of values of δ only down to 10⁻⁴, but the trend of their simulations indicates that δ = 10⁻⁵ would lead to clumping for values corresponding to our resulting particle sizes. A low turbulence like this would also decrease the effect of collisional fragmentation, as is clear from Fig. 5.

We note that the turbulence level of δ = 10⁻⁵ that we adopted here is lower than in many simulations, and that a high background diffusion would restrict the effectiveness of clumping by the streaming instability (Umurhan et al. 2020). However, there is observational support for a low level of turbulence (Villenave et al. 2022). Theoretically, we also expect the mid-plane layer to be less turbulent than the upper disc layers, as it has been shown that the vertical shear instability stirs the upper layers of the disc, but not the mid-plane (Schäfer et al. 2020; Schäfer & Johansen 2022).

Finally, we need to assess the amount of material available for clumping into planetesimals at the ice line. The mass in the annulus where particles grow above the streaming instability threshold can be estimated in Earth masses, M_E, as $$\begin{array}{rl}M& =2\pi rW{\mathrm{\Sigma }}_{\text{peb }}\\ & \approx 0.4M_{\mathrm{E}}\left(\frac{r}{1.8\mathrm{a}\mathrm{u}}\right)\left(\frac{W}{0.1\mathrm{a}\mathrm{u}}\right)\left(\frac{{\mathrm{\Sigma }}_{\mathrm{g}}}{{10}^4\mathrm{k}\mathrm{g}{\mathrm{m}}^{-2}}\right)\left(\frac{Z}{0.01}\right),\end{array}$$(27)

where ∼_peb is the surface density of pebbles, and W is the width of the annulus. Estimating the width of the annulus to W = 0.2 au, and assuming that both ice and rock contributes to the clumping, we thus determine a mass of the material available for clumping of M = 0.8 M_E. This confirms the potential for a pathway towards planets around the ice line.

Fig. 6 Threshold for particle clumping via the streaming instability for different particle stopping times and metallicities. The green area in the upper right part denotes the parameter space within which clumping via streaming instability is possible. The filled coloured circles show our data for different disc locations with a cooling time of tcool = 100 yr. The lower row corresponds to a disc with solar metallicity within which the interacting solids are made of ice, and the upper row is either a high-metallicity disc or a disc with solar metallicity, where both ice and rock contribute to clumping.

5 Conclusions

We have investigated the effects of stellar outbursts, specifically, FU Orionis outbursts, and the subsequent cooling of the disc on heterogeneous ice nucleation onto silicate particles and particle growth by deposition. During an outburst, the temperature in the disc increases dramatically, causing the ice line to shift outwards and to sublimate all the ice in the affected region. The subsequent cooling phase allows the ice line to move back to its quiescent position, and vapour is then redistributed onto solid grains. Since nucleation of ice from vapour onto bare grains requires a higher supersaturation than the deposition of vapour onto existing ice grains, this results in the condensation of ice on a few particles that can then grow to be very large. Our main findings are summarised below.

• 1.

  Particle growth: a stellar outburst resets the particle size distribution so that particle growth starts over from bare silicate grains. Beginning from a particle distribution with a maximum size of R = 1 mm, we find resulting sizes of several centimetres just outside the quiescent ice line. Growth nevertheless affects a broad region between r_in ≈ 2 au and r_out ≈ 10 au. Importantly, growth by condensation following heterogeneous nucleation is several times faster than growth by coagulation, meaning that it can be less dependent on typical collisional growth barriers.

• 2.

  Cooling time dependence: the timescale on which the disc cools is important for the resulting particle sizes. A slower cooling leads to larger particles. However, for both the long and short cooling timescales that we investigated, we find growth to at least centimetre-sized pebbles outside of the ice line. The width of the growth zone and the maximum particle size are nevertheless both larger in the case of slower cooling. For the longer cooling time investigated in this paper, t_cool = 1000 yr, centimetre-sized particles form out to 3 au, whereas the shorter cooling time of t_cool = 100 yr results in growth to at least centimetre-sized particles out to 2 au.

• 3.

  Collisions and sedimentation: we find that both collisions and sedimentation can be ignored in the modelling of nucleation and deposition in the aftermath of an FU Orionis outburst because of the rapid condensational growth that follows heterogeneous nucleation. Collisions between particles of a few centimetres and larger could nevertheless lead to fragmentation for typical turbulence levels after the phase of condensational growth.

• 4.

  Relevance for the streaming instability: after an outburst, we find particles of several centimetres outside the ice line. This corresponds to Stokes numbers as high as St = 0.01–0.02, which is sufficient to trigger the streaming instability for a metallicity of Z = 10⁻² in a low-turbulence disc (δ = 10⁻⁵). We did not model the actual step from particle growth through sedimentation and to the streaming instability, but we tentatively find that stellar outbursts lead to favourable conditions for growth into planetesimals.

• 5.

  Importance of nucleation: heterogeneous nucleation is a crucial component in particle growth by ice condensation. The final particle sizes depend on how many particles are available for vapour deposition (typically referred to as condensation in astronomical literature), with larger particles being favoured in a scenario in which nucleation is rare. We have demonstrated that the aftermath of a stellar outburst with relatively slow cooling (100 yr ≲ t_cool ≲ 1000 yr) provides just such favourable conditions as needed for particle growth. Importantly, the condensation growth that follows nucleation is so rapid that dynamical processes such as diffusion, sedimentation, and drift have little or no influence on the ensuing particle sizes.

We did not model the collisional evolution of the size distribution of silicate particles and solid-ice particles after the outburst. The largest ice particles will collide fast enough to experience cracking collisions, which will reduce their size to a few centimetres at most. At the same time, silicate dust aggregates colliding with solid-ice particles may adhere their monomers to the icy surface. Laboratory experiments on the outcome of collisions between dust aggregates and macroscopic solid water ice particles are lacking in the literature so far. Additional modelling that takes the coexistence of dust aggregates and solid ice into account as well as laboratory experiments are therefore needed to better understand particle growth and planetesimal formation outside the water ice line after FU Orionis outbursts.

Acknowledgements

We thank the referee for a thorough reading of the manuscript and for providing comments and questions that helped improve the quality of this paper. A.J. acknowledges funding from the European Research Foundation (ERC Consolidator Grant 724687-PLANETESYS), the Knut and Alice Wallenberg Fountdation (Wallenberg Scholar Grant 2019.0442), the Swedish Research Council (Project Grant 2018-04867), the Danish National Research Foundation (DNRF Chair Grant DNRF159) and the Göran Gustafsson Foundation.

References

All Figures

	Fig. 1 Nucleation timescale for a single 100-micrometre-sized grain for a range of super-saturations (S) and temperatures (T). The dotted line marks the 1 Myr nucleation timescale, and the dashed line shows the 1 Gyr nucleation timescale. We illustrate a schematic cooling path at the conditions of the water ice line with the green line.
In the text

	Fig. 2 Mean sizes and Stokes number for particles after an FU Orionis outburst with a decay time of 100 yr in blue and 1000 yr in red. The sizes of silicate particles at the start of the simulation are shown as dotted lines, with a maximum size of one millimetre. Growth to the largest particles occurs just outside of the quiescent ice line.
In the text

Fig. 3 Mass distribution of particles (normalised such that the integral over the logarithmic radius is always unity) just outside of the ice line at 1.8 (yellow line), 2.1 (red line), and 3.6 au (blue line) after the ice line has returned to its quiescent state in a scenario with a cooling time of tcool = 100 yr. We indicate the particle size for the longer cooling time of 1000 yr with circles. The black line shows the mass distribution of bare silicate particles.

In the text

Fig. 4 Evolution of a representative pebble at 1.8 au, nucleating when S = 2, corresponding to a time of t ≈ 840 yr. Panel a shows the super-saturation S at 1.8 au. Panel b shows the particle size in metres and Stokes number. We mark St = 0.01 with a dotted line. Panel c shows the growth rate in millimetres per year, and panel d shows the sedimentation and growth timescales. The cooling time of the disc is tcool =100 yr.

In the text

	Fig. 5 Ice particle collision speed as a function of the mean ice particle size in each radial bin for a turbulent diffusion coefficient of δ = 10−3, δ = 10−4, and δ = 10−5. The critical velocity for ice fragmentation is over-plotted as a dotted line. Ice particles can grow to a few centimetres before they risk fragmentation in collisions.
In the text

Fig. 6 Threshold for particle clumping via the streaming instability for different particle stopping times and metallicities. The green area in the upper right part denotes the parameter space within which clumping via streaming instability is possible. The filled coloured circles show our data for different disc locations with a cooling time of tcool = 100 yr. The lower row corresponds to a disc with solar metallicity within which the interacting solids are made of ice, and the upper row is either a high-metallicity disc or a disc with solar metallicity, where both ice and rock contribute to clumping.

In the text