© ESO, 2013

1. Introduction

Planets form in protoplanetary discs of gas and dust that surround young stars. In the classical planet formation scenario dust grains collide, stick together, and form larger and larger bodies (Safronov 1969). Particles of up to millimetre sizes stick together owing to contact forces and kilometre-sized and larger bodies are held together by gravity. How particles grow from millimetre to kilometre-sized planetesimals is difficult to explain by collisions, since these particles tend to fragment or bounce off each other when they collide, instead of sticking (Blum & Wurm 2008). Furthermore, solids on Keplerian orbits experience a headwind from the slow-orbiting pressure-supported gas, and therefore lose angular momentum and drift in towards the central star (Weidenschilling 1977). This radial drift velocity peaks for metre-sized particles, which drift into the star from 1 AU on a time scale of only 100 orbits (Brauer et al. 2008).

Once particles grow to centimetre-sized to decimetre-sized pebbles, further growth into planetesimals is possible via particle concentration and gravitational collapse. The streaming instability causes particles to clump together owing to the velocity difference between the particles and the gas in the disc, and planetesimals form through subsequent gravitational collapse (Youdin & Goodman 2005; Johansen et al. 2007). Particles also concentrate in long-lived vortices (Barge & Sommeria 1995; Klahr & Bodenheimer 2003) and in pressure bumps (Johansen et al. 2009) excited in the turbulent flow. However, an efficient mechanism for growth from millimetre-sized dust grains to pebbles is needed to explain the formation of particles that are large enough to take part in such concentration events.

Sticking, or coagulation, as a growth mechanism has been extensively studied, both experimentally (Blum & Wurm 2008; Zsom et al. 2010) and numerically (Brauer et al. 2008; Windmark et al. 2012). However, an often overlooked concept in planet formation models is that of growth via condensation near ice lines. Stevenson & Lunine (1988) investigated material enhancement at the ice line due to diffusive redistribution and subsequent condensation of water vapour from the inner part of the disc. With the assumption of efficient condensation by homogeneous nucleation and by ignoring inward transport of material, the authors found a significant solid density enhancement. In contrast, Cuzzi & Zahnle (2004) considered material drifting inwards across the ice line, which would enhance the vapour density in the inner disc, followed by a phase of accumulation of solids by condensation onto immobile planetesimals formed outside of the ice line.

In this work the dynamical behaviour and growth of small (~1 mm) seed particles via ice condensation is studied in computer simulations. Ice particles moving with the turbulent gas encounter the water ice line where ice sublimates and becomes water vapour. We consider both the radial ice line, which separates the inner, hotter part of the disc from the colder region farther away from the star, and the atmospheric ice line, which separates the hot atmosphere of the disc from the cold midplane layer (Dullemond & Dominik 2004; Meijerink et al. 2009). Owing to the turbulent motion of the gas in the protoplanetary disc, water vapour diffuses across the ice line into the condensation region, where it condenses onto existing ice particles. Many ice particles will move across the ice line and sublimate. However, since small particles are coupled to the gas, thereby effectively moving in a random walk due to turbulence, a significant fraction will be lucky enough to stay within the condensation region of the disc, growing to larger sizes as diffusing water vapour condenses onto them. The total mass in the ice component is maintained during this process, because the loss of ice particles to sublimation is completely compensated for by the growth of the remaining particles.

In a moderately turbulent disc, we find particle growth to beyond decimetre sizes on a time scale of a few thousand orbits. The most significant growth takes place locally, close to the radial ice line. However, efficient radial mixing supplies the midplane with large particles, even at distances of several scale heights outside of the ice line. The atmospheric ice line, located at z ≳ 3 H where disc material is heated directly by the central star (Chiang & Goldreich 1997), contributes less to the particle growth, generally of the order of 1 μm.

The radial water ice line, or condensation front, is found at r ≈ 3 AU around a solar type star (Lecar et al. 2006). However, growth by condensation is not only applicable to water ice, but also to any other volatile found in the protoplanetary disc. Other important condensation fronts relevant for planet formation include those of ammonia, methane, carbon monoxide and molecular nitrogen at much larger radii than the water ice line, and of silicates at r ≪ 1 AU (Lodders 2003; Qi et al. 2011). These condensation fronts will be the topic of a future study.

The paper is organised as follows. In Sect. 2 the numerical model used for particle dynamics and growth by condensation is described, including the units used in this paper. The Monte Carlo scheme for condensation is further explained in Sect. 3. In Sect. 4 we describe the test problems used to validate the code. We present our results in Sect. 5, and investigate the influence of the position of the atmospheric ice line, the turbulence strength and the effects of radial mixing. In a moderately turbulent disc (α = 10^-2), pebbles of centimetres to decimetres form within a few thousand years. We discuss general assumptions and simplifications made in Sect. 6. Finally, in Sect. 7 we conclude that particle growth by condensation is an important mechanism for forming pebbles in protoplanetary discs, which should be taken into account in future studies of planetesimal formation.

2. Numerical model

2.1. Simulation box, units and boundary conditions

We simulate a two-dimensional region around the water ice line. The simulation domain is set in the radial r and vertical z direction, whereas the azimuthal direction is ignored, for simplicity and under the assumption of axisymmetry. The fundamental units are the gas scale height H as a length unit, sound speed c_s as a velocity unit and inverse Keplerian orbits Ω^-1 as a time unit. The relation between these quantities is H = c_s/Ω. Setting H = c_s = Ω = 1 gives a system that is scalable to any region of the protoplanetary disc. Hence our results apply equally well to other condensation fronts, such as those of CO or silicates, although the transition and scaling between drag regimes and the fractional ice abundance depend on the choice of radial location in the disc. Here we interpret the results in terms of the water ice line at around 3 AU from the star, using a water ice mass fraction of 0.571% (Lodders 2003).

Condensation is considered as a neighbour interaction. A linear scaling between number of particles and number of calculations is obtained by the use of a mapping scheme, where the simulation domain is divided into a number of grid cells, with particle interactions possible only between particles within the same grid cell.

Initially, both vapour and ice particles have a Gaussian distribution in the vertical direction, following the gas distribution. At the beginning of a simulation particles are set to be ice or vapour, depending on whether they are located within or outside of the condensation region.

2.2. Particle sizes

Particles are coupled to the gas via drag forces and can be characterised by their dimensionless friction time (Weidenschilling 1977), or Stokes number, $\begin{array}{ccc}& \mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}^{\mathrm{\left(}\mathrm{Ep}\mathrm{\right)}}\mathrm{=}\frac{\mathit{R}{\mathit{\rho }}_{\mathrm{•}}}{\mathit{H}{\mathit{\rho }}_{\mathrm{g}}}& \\ & \mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}^{\mathrm{\left(}\mathrm{St}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{4}}{\mathrm{9}}\frac{{\mathit{R}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{•}}}{\mathit{H}{\mathit{\rho }}_{\mathrm{g}}\mathit{\lambda }}& \end{array}$where λ is the mean free path of the gas, ρ_• is the material density and ρ_g is the gas density. The superscripts (Ep) and (St) denote Epstein and Stokes drag regime, the two drag regimes relevant for the particle sizes and gas density considered in this paper. We formulate the particle radius R in units of R₁, where $\mathit{R}\mathit{/}{\mathit{R}}_{\mathrm{1}}\mathrm{=}\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}^{\mathrm{\left(}\mathrm{Ep}\mathrm{\right)}}$ so that $\begin{array}{ccc}{\mathit{R}}_{\mathrm{1}}\mathrm{=}\frac{\mathit{H}{\mathit{\rho }}_{\mathrm{g}}}{{\mathit{\rho }}_{\mathrm{•}}}\mathrm{·}& & \end{array}$(3)To find the corresponding particle size in metres, we use the material density of ice ρ_• ≈ 1 g cm^-3, the midplane gas density ρ_g and column density Σ from the minimum mass solar nebula (MMSN) model of Hayashi (1981), $\begin{array}{ccc}& & {\mathit{\rho }}_{\mathrm{g}}\mathrm{=}\frac{\mathit{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\sqrt{\mathrm{2}\mathit{\pi }}\mathit{H}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathit{,}\\ & & \end{array}$giving the scaling of $\begin{array}{ccc}{\mathit{R}}_{\mathrm{1}}\mathrm{\approx }\mathrm{1.3}\hspace{0.17em}\mathrm{m}{\left(\frac{\mathit{r}}{\mathrm{3}\hspace{0.17em}\mathrm{AU}}\right)}^{-1.5}& & \end{array}$(6)close to the water ice line.

2.3. Condensation scheme

The rate of change in particle mass caused by condensation and sublimation is $\begin{array}{ccc}\frac{\mathrm{d}\mathit{m}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{R}}^{\mathrm{2}}{\mathit{v}}_{\mathrm{th}}{\mathit{\rho }}_{\mathrm{v}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{P}}_{\mathrm{sat}}}{{\mathit{P}}_{\mathrm{v}}}\right)\hspace{0.17em}& & \end{array}$(7)(Supulver & Lin 2000). Here, v_th is the thermal velocity of vapour, ρ_v is the vapour density and P_sat and P_v is the saturated vapour pressure and vapour pressure, respectively. Both the vapour pressure and the saturated vapour pressure can be expressed as densities using the ideal gas law. The vapour pressure can be written as $\begin{array}{ccc}{\mathit{P}}_{\mathrm{v}}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{B}}\mathit{T}}{{\mathit{m}}_{\mathrm{v}}}{\mathit{\rho }}_{\mathrm{v}}\mathit{,}& & \end{array}$(8)and the saturated vapour pressure can be expressed in an equivalent way. Here, k_B is the Boltzmann constant, T is the temperature, ρ_v is the vapour density and m_v is the mass of one vapour particle. Assuming spherical particles and using vapour densities instead of pressures, Eq. (7) can be rewritten as $\begin{array}{ccc}\frac{\mathrm{d}\mathit{R}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{th}}}{{\mathit{\rho }}_{\mathrm{•}}}\left({\mathit{\rho }}_{\mathrm{v}}\mathrm{-}{\mathit{\rho }}_{\mathrm{sat}}\right)\mathit{,}& & \end{array}$(9)in the limit where ρ_v ≪ ρ_sat sublimation dominates the time evolution. The sublimation time scale in terms of particle radius is then $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{s}\mathit{,}\hspace{0.17em}\mathrm{R}}\mathrm{=}\frac{\mathit{R}{\mathit{\rho }}_{\mathrm{•}}}{{\mathit{v}}_{\mathrm{th}}{\mathit{\rho }}_{\mathrm{sat}}}\mathit{,}& & \end{array}$(10)and in terms of mass $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{s}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\frac{\mathit{R}{\mathit{\rho }}_{\mathrm{•}}}{{\mathit{v}}_{\mathrm{th}}{\mathit{\rho }}_{\mathrm{sat}}}\mathrm{·}& & \end{array}$(11)Correspondingly, when ρ_sat ≪ ρ_v condensation dominates. The condensation time scale in terms of radius is $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{c}\mathit{,}\hspace{0.17em}\mathrm{R}}\mathrm{=}\frac{\mathit{R}{\mathit{\rho }}_{\mathrm{•}}}{{\mathit{v}}_{\mathrm{th}}{\mathit{\rho }}_{\mathrm{v}}}\mathit{,}& & \end{array}$(12)and in terms of mass $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{c}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\frac{\mathit{R}{\mathit{\rho }}_{\mathrm{•}}}{{\mathit{v}}_{\mathrm{th}}{\mathit{\rho }}_{\mathrm{v}}}\mathrm{·}& & \end{array}$(13)In equilibrium ρ_v = ρ_sat with τ_s = τ_c. By definition ρ_v = f_H2Oρ_g, where f_H2O is the water mass fraction, at an ice line. The rapid increase of temperature into the condensation region implies ρ_sat ≫ ρ_g and hence τ_s ≪ τ_f, allowing for modelling of sublimation as an instantaneous process as long as the particle friction time, τ_f, is shorter than the characteristic time scale of the turbulent eddies, Ω^-1.

Condensation can not be assumed to occur instantaneously, since ρ_v never increases much beyond f_H2Oρ_g. The condensation time scale is proportional to the size of the particle, with the dimensionless condensation time scale being $\begin{array}{ccc}\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{c}}\mathrm{=}\frac{\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}^{\mathrm{\left(}\mathrm{Ep}\mathrm{\right)}}}{\mathrm{3}}\frac{{\mathit{\rho }}_{\mathrm{g}}}{{\mathit{\rho }}_{\mathrm{v}}}\mathit{,}& & \end{array}$(14)assuming ${\mathit{c}}_{\mathrm{s}}\mathrm{\approx }\sqrt{\frac{{\mathit{\mu }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{0}}}{\mathit{\mu }}}{\mathit{v}}_{\mathrm{th}}$, where μ is the mean molecular weight of molecular hydrogen and μ_H2O is the corresponding quantity for water. The condensation process is modelled by a Monte Carlo approach, where the probability of condensation of vapour onto an ice particle in the condensation zone is proportional to the size of the available ice surface. This is further described in Sect. 3.

2.4. Superparticle approach

The ice and vapour components are modelled using a superparticle approach, related to the coagulation algorithm of Zsom & Dullemond (2008), where a superparticle is a numerical representation of a large number of physical particles with identical properties. Here the most important properties are the physical state (ice or vapour), size of constituent particles if ice, and material density. The number of superparticles N is much smaller than the number of physical particles, so that superparticle i represents n_i physical particles. In these simulations, typically N = 10 000 for a simulation area of 12 H in the vertical direction and 1.5 H in the radial direction, a number chosen to be computationally easy to handle. The number of superparticles is fixed throughout a simulation, and each superparticle represents an equal and fixed mass M that does not change. The mass m of the physical particles represented by the superparticle does however change, and so does n_i, the number of physical particles represented by a superparticle, since the total mass represented by the superparticle is always M = m_i n_i. Also when an ice particle undergoes a phase transition and becomes vapour the total mass M of the superparticle stays constant. The superparticle approach to condensation and sublimation is sketched in Fig. 1.

Fig. 1 Sketch of sublimation and condensation for superparticles. Blue represents ice and red vapour. Left panel: ice superparticle, representing a mass M in the form of 4 physical ice particles, each of mass m, sublimates and turns into a vapour superparticle of mass M. Right panel: vapour superparticle of mass M condenses onto an ice superparticle of mass M, representing the mass 4 m. The physical ice particles are after the event represented by two ice superparticles, each representing a total mass M, but now in the form of 2 physical ice particles each with mass 2 m. Both mass and the number of superparticles are conserved between t0 and t1. In the case of condensation, also the number of physical particles is conserved.

2.5. Random walk of particles

In this section the modelling of particle motions is outlined. We start by introducing the random walk used to describe the turbulent motions of the gas, sufficient for modelling small particle motions, and then move on to include the more complex behaviour of larger particles. When particles grow larger, additional effects need to be taken into account. Firstly, the random step length becomes smaller, since larger particles follow the turbulent eddies a shorter length each time step. Secondly, for larger particles gravity towards the midplane has an increasingly important effect, and causes particles to sediment towards the midplane. Finally, radial drift towards the central star must be taken into account.

With the dimensionless constant α introduced by Shakura & Sunyaev (1973), the turbulent diffusion coefficient of gas and small particles in an accretion disc can be written as $\mathit{D}\mathrm{=}\mathit{\alpha }{\mathit{c}}_{\mathrm{s}}\mathit{H}\mathit{.}$(15)We here assume that particles move with velocity $\mathit{v}\mathrm{~}\sqrt{\mathit{\alpha }}{\mathit{c}}_{\mathrm{s}}$, and that the length a particle moves coherently is $\mathit{l}\mathrm{~}\sqrt{\mathit{\alpha }}\mathit{H}$, which gives rise to D ~ vl ~ αc_sH by mixing length arguments. The correlation time is defined as $\mathit{\tau }\mathrm{=}\frac{\mathit{l}}{\mathit{v}}\mathrm{=}{\mathrm{\Omega }}^{-1}\mathit{.}$(16)This means that the effect of turbulence on a small particle is to move it in a random direction with a characteristic speed v, during the correlation time τ. The characteristic length scale, velocity scale and time scale are all set by the turbulence. We use a default α-value of 10^-2, motivated by observations of accretion rates in young stars (Hartmann et al. 1998). However, the presence of a dead zone can give a lower α-value (Fleming & Stone 2003; Oishi & Mac Low 2009), and there are also observations suggesting α = 10^-4 in protoplanetary discs (Mulders & Dominik 2012). Therefore we also run simulations with α = 10^-3 and α = 10^-4.

The distance a particle moves each time step, where one time step is Δt = Ω^-1, is found by equating the generic expression for the diffusion coefficient of a random walk in 2D, D = l²/(4τ), with that of diffusion described by Eq. (15), giving $\mathit{l}\mathrm{=}\mathrm{2}\sqrt{\mathit{\alpha }}\mathit{H}\mathit{.}$(17)Practically, the random walk for particles coupled to the gas is modelled by setting the length the particle moves in the radial direction Δr and in the vertical direction Δz to $\begin{array}{ccc}\mathrm{(}\mathrm{\Delta }\mathit{r}{\mathrm{)}}_{\mathrm{rw}}& \mathrm{=}& \\ \mathrm{(}\mathrm{\Delta }\mathit{z}{\mathrm{)}}_{\mathrm{rw}}& \mathrm{=}& \end{array}$Here 0 ≤ θ < 2π is chosen randomly for each particle and time step, so that a particle moves in a random direction, but a fixed length, each time step. This forms the basis of the particle dynamics used in this code, and is valid for small particles perfectly coupled to the turbulent gas. However, for larger particles this simple approximation breaks down. A number of corrections to the random walk model will now be introduced, one at a time, before arriving at the final model for the particle dynamics.

Firstly, the random step length must be adjusted, since larger particles move shorter distances with the turbulent eddies each time step. The diffusion coefficient D = vl = l²/τ = l²Ω is modified to ${\mathit{D}}^{\mathrm{\ast }}\mathrm{=}\frac{\mathit{D}}{\mathrm{1}\mathrm{+}\mathrm{(}\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}{\mathrm{)}}^{\mathrm{2}}}\hspace{0.17em}$(20)following the analytical theory developed and tested by Youdin & Lithwick (2007). This gives a modified size-dependent step length of ${\mathit{l}}^{\mathrm{\ast }}\mathrm{=}\frac{\mathrm{2}\sqrt{\mathit{\alpha }}\mathit{H}}{\sqrt{\mathrm{1}\mathrm{+}\mathrm{(}\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}{\mathrm{)}}^{\mathrm{2}}}}\mathrm{=}\frac{\mathit{l}}{\sqrt{\mathrm{1}\mathrm{+}\mathrm{(}\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}{\mathrm{)}}^{\mathrm{2}}}}\mathrm{·}$(21)For small particles l^∗ ≈ l, whereas for larger particles l^∗ < l. Separating the radial and vertical directions, Eqs. (18) and (19) are replaced by $\begin{array}{ccc}\mathrm{(}\mathrm{\Delta }\mathit{r}{\mathrm{)}}_{\mathrm{rw}}& \mathrm{=}& \frac{\mathit{l}\hspace{0.17em}\cos \mathit{\theta }}{\sqrt{\mathrm{1}\mathrm{+}\mathrm{(}\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}{\mathrm{)}}^{\mathrm{2}}}}\mathit{,}\\ \mathrm{(}\mathrm{\Delta }\mathit{z}{\mathrm{)}}_{\mathrm{rw}}& \mathrm{=}& \end{array}$When particles grow larger, and drag forces decrease in strength, the gravity towards the midplane influences the particle motion more, so that large particles sediment to the midplane. Since this only affects the vertical direction, it is only necessary to modify Δz. This is done by changing the step size in the vertical direction to (Δz)_tot = (Δz)_rw + (Δz)_sed, where (Δz)_sed is an increasing function of particle size. The sedimentation step length is set by the terminal velocity of the particles v_t =  −τ_fΩ²z. However, in order to modify the step length correctly, we must ensure that the resulting particle scale height is the same as the analytically expected scale height. The analytically expected particle scale height, for a given particle size and turbulence strength, resulting from diffusion and sedimentation is $\frac{{\mathit{H}}_{\mathrm{p}}}{\mathit{H}}\mathrm{=}\sqrt{\frac{\mathit{\alpha }}{\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}\mathrm{+}\mathit{\alpha }}}\mathit{,}$(24)as determined in computer simulations of particles in turbulent protoplanetary discs (Carballido et al. 2006) and explained in the analytical diffusion framework of Youdin & Lithwick (2007). This can now be compared with the particle scale height that arises from considering a combination of a random walk and terminal velocity sedimentation, which is what is used for modelling. The random walk of particles mimics turbulent diffusion with coefficient D = αc_sH, and sedimentation-diffusion equilibrium occurs at particle scale height $\frac{{\mathit{H}}_{\mathrm{p}}}{\mathit{H}}\mathrm{=}\sqrt{\frac{\mathit{\alpha }}{\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}}}\mathrm{·}$(25)To construct a model that gives the correct scale height, as given by Eq. (24), we therefore need to modify the dimensionless friction time in the terminal velocity expression to $\mathrm{(}\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}{\mathrm{)}}^{\mathrm{\ast }}\mathrm{=}\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}\mathrm{+}\mathit{\alpha }\mathit{.}$(26)Effectively we let even tiny particles sediment in order to follow the Gaussian distribution of the gas, rather than the uniform distribution that results from a pure random walk. Changing Ωτ_f → (Ωτ_f)^∗ in Eq. (25) reproduces the correct particle scale height in Eq. (24). With the modified dimensionless friction time (Ωτ_f)^∗ in the terminal velocity expression, the step size in the vertical direction can be written as $\mathrm{(}\mathrm{\Delta }\mathit{z}{\mathrm{)}}_{\mathrm{tot}}\mathrm{=}\mathrm{(}\mathrm{\Delta }\mathit{z}{\mathrm{)}}_{\mathrm{rw}}\mathrm{-}\mathrm{(}\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}\mathrm{+}\mathit{\alpha }\mathrm{)}\mathrm{\Omega }\mathit{z}\mathrm{\Delta }\mathit{t}\mathit{.}$(27)However, this expression leads to a numerical instability in which large particles overshoot the midplane while they sediment, since for Ωτ_f > 1 and Δt = Ω^-1 we get $\mathrm{(}\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}{\mathrm{)}}^{\mathrm{\ast }}\mathrm{\Omega }\mathit{z}\mathrm{\Delta }\mathit{t}\mathit{>}\mathit{z}\mathit{.}$(28)This gives particle oscillations that are amplified in time. Equation (27) is therefore modified to $\mathrm{(}\mathrm{\Delta }\mathit{z}{\mathrm{)}}_{\mathrm{tot}}\mathrm{=}\mathrm{(}\mathrm{\Delta }\mathit{z}{\mathrm{)}}_{\mathrm{rw}}\mathrm{-}\frac{\mathrm{(}\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}\mathrm{+}\mathit{\alpha }\mathrm{)}\mathrm{\Omega }\mathit{z}\mathrm{\Delta }\mathit{t}}{\mathrm{1}\mathrm{+}\mathrm{(}\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}{\mathrm{)}}^{\mathrm{2}}}\mathit{,}$(29)following the scheme of Youdin & Lithwick (2007), which eliminates the amplifying of the vertical particle oscillations. Instead it gives larger particles a longer settling time, corresponding to the time they would have spent oscillating before settling. For a comparison between the particle scale height resulting from Eq. (29) and the theoretical scale height from Eq. (24), see Fig. 2. The correlation between the analytical and modelled curve in this figure shows that the particle dynamics in the vertical direction are correctly modelled by Eq. (29).

The step in the radial direction needs to be modified because of the radial drift towards the central star. This is also particle size dependent, and (Δr)_tot is modified to $\mathrm{(}\mathrm{\Delta }\mathit{r}{\mathrm{)}}_{\mathrm{tot}}\mathrm{=}\mathrm{(}\mathrm{\Delta }\mathit{r}{\mathrm{)}}_{\mathrm{rw}}\mathrm{+}\mathrm{(}\mathrm{\Delta }\mathit{r}{\mathrm{)}}_{\mathrm{rd}}\mathrm{=}\mathrm{(}\mathrm{\Delta }\mathit{r}{\mathrm{)}}_{\mathrm{rw}}\mathrm{-}\frac{\mathrm{2}\mathit{\eta }{\mathit{v}}_{\mathrm{K}}}{\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}\mathrm{+}\mathrm{(}\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}{\mathrm{)}}^{-1}}\mathrm{\Delta }\mathit{t,}$(30)where v_K is the Keplerian velocity and ηv_K is the velocity difference between the gas and the particles, following Weidenschilling (1977). This gives a radial drift velocity that peaks for Ωτ_f = 1 particles. We use ηv_K = 0.05c_s, a reasonable value at 3 AU in the MMSN (Cuzzi et al. 1993; Chiang & Youdin 2010).

The friction time of a particle is inversely proportional to the gas density, $\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}\mathrm{\propto }{\mathit{\rho }}_{\mathrm{g}}^{-1}$, which decreases in the vertical direction with distance from the midplane. Taking into account the Gaussian distribution of the gas around the midplane, the dimensionless friction time of a particle is modified according to $\begin{array}{ccc}\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}\mathrm{\to }\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}\mathit{ϵ}\hspace{0.17em}\mathit{,}& & \end{array}$(31)where $\begin{array}{ccc}\mathit{ϵ}\mathrm{=}\exp \mathrm{[}{\mathit{z}}^{\mathrm{2}}\mathit{/}{\mathrm{(}\mathrm{2}{\mathit{H}}^{\mathrm{2}}\mathrm{)}}^{\mathrm{]}}\mathit{.}& & \end{array}$(32)Simulations have shown that the heterogeneous vertical gas distribution leads to a turbulent α-value that increases with height above the midplane (Fromang & Nelson 2009). We do however not take this effect into account in our simulations.

For large simulation domains radial gas density gradients are also of importance, and the random walk would then have to be adjusted in order to mimic the motion of especially small particles, well coupled to the gas (Hughes & Armitage 2010). Since the model used in this paper is local, we ignore this effect.

The equations used to describe particle dynamics in the code are thus $\begin{array}{ccc}\mathrm{\Delta }\mathit{z}& \mathrm{=}& \\ \mathrm{\Delta }\mathit{r}& \mathrm{=}& \end{array}$where the first term for both the radial and vertical direction describe the random walk, with a step length modified according to particle size, and the second term describes sedimentation and radial drift, for the vertical and radial direction, respectively.

2.6. Transition between drag regimes

For small particles with a radius of less than (9/4) of the gas mean free path, Epstein drag applies, whereas for larger particles Stokes drag is the relevant regime. The transition size given by Eqs. (1) and (2) can be rewritten in dimensionless friction time as $\begin{array}{ccc}{\left(\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}\right)}^{\mathrm{\left(}\mathrm{Ep}\mathrm{\right)}}\mathit{<}\frac{\mathrm{9}}{\mathrm{4}}\frac{{\mathit{\rho }}_{\mathrm{•}}}{{\mathit{\rho }}_{\mathrm{g}}}\frac{\mathit{\lambda }}{\mathit{H}}\mathrm{=}\frac{\mathrm{9}\mathit{\pi }}{\mathrm{2}}\frac{{\mathit{\rho }}_{\mathrm{•}}\mathit{\mu H}}{{\mathit{\Sigma }}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{mol}}}\mathrm{\approx }\mathrm{0.5}\mathit{,}& & \end{array}$(35)when applied to the midplane of the disc. For molecular hydrogen, the mean molecular weight is μ = 3.9 × 10^-24 g and the molecular cross-section is σ_mol = 2 × 10^-15 cm² (Nakagawa et al. 1986). The column density Σ at 3 AU is given by Eq. (5).

Converting dimensionless friction time in the Epstein drag regime to the Stokes drag regime is done using $\begin{array}{ccc}\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}^{\mathrm{\left(}\mathrm{St}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{9}\mathit{\pi }}\frac{{\mathit{\Sigma }}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{mol}}}{{\mathit{\rho }}_{\mathrm{•}}\mathit{\mu H}}{\mathrm{(}\mathrm{\Omega }{{\mathit{\tau }}_{\mathrm{f}}^{\mathrm{\left(}\mathrm{Ep}\mathrm{\right)}}}^{\mathrm{)}}}^{\mathrm{2}}\mathrm{\approx }\mathrm{2.0}{\mathrm{(}\mathrm{\Omega }{{\mathit{\tau }}_{\mathrm{f}}^{\mathrm{\left(}\mathrm{Ep}\mathrm{\right)}}}^{\mathrm{)}}}^{\mathrm{2}}\mathit{,}& & \end{array}$(36)where the scale height is given by the MMSN model (Hayashi 1981) as $\begin{array}{ccc}\mathit{H}\mathrm{=}\mathrm{0.033}{\left(\frac{\mathit{r}}{\mathrm{AU}}\right)}^{\mathrm{1.25}}\mathrm{·}& & \end{array}$(37)

3. Condensation in a Monte Carlo scheme

We treat condensation as a neighbour interaction, and require vapour to be present near an ice particle for interaction to occur. In this model “near” means being in the same grid cell, so that it is possible for an ice particle to interact only with vapour within the same grid cell. Condensation in each time step is thus decided one grid cell at a time. Any grid cell can, and most often does, harbour more than one ice particle. Therefore, whether or not condensation takes place in a given grid cell and a given time step, is decided in a two-step process for each vapour particle present in the grid cell.

The two-step procedure is here described for one grid cell and one time step Δt, with one vapour particle present in this grid cell. In the case where more than one vapour particle is present in the grid cell, the two-step procedure is repeated for each vapour particle present. The first step is to decide whether this vapour particle condenses onto any of the ice particles. If it does, the second step is to decide which of the ice particles it condenses onto. In this scheme, a vapour particle can only condense onto one single ice particle, and thus can not be shared amongst several particles. One ice particle can on the other hand experience several growth events during one time step, provided that several vapour particles are present.

As a first step it is decided whether or not the vapour particle is involved in a condensation event during the time step Δt. In order to do this the total interaction probability for all ice particles in the grid cell is needed. The interaction probability for ice particle i is ${\mathit{p}}_{\mathit{i}}\mathrm{=}\mathrm{1}\mathrm{-}\exp \mathrm{(}\mathrm{-}\mathrm{\Delta }\mathit{t}\mathit{/}{\mathit{\tau }}_{\mathit{i}}\mathrm{)}\mathit{,}$(38)where τ_i is the condensation time scale for particle i, given by Eq. (14). The expression for interaction probability is chosen to always give 0 ≤ p_i < 1, for any Δt, also for Δt ≫ τ_i. The algorithm proceeds by calculating the probability that no interaction occurs. For n ice particles, the total probability that no interaction occurs is ${\mathit{p}}_{\mathrm{0}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{p}}_{\mathrm{1}}\mathrm{)}\mathrm{·}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{p}}_{\mathrm{2}}\mathrm{)}\mathrm{·}\hspace{0.17em}\mathrm{·}\mathrm{·}\mathrm{·}\hspace{0.17em}\mathrm{·}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{p}}_{\mathit{n}}\mathrm{)}\mathit{.}$(39)To decide whether condensation occurs, a random number r₁ ∈ ] 0,1 [ is generated. If r₁ is smaller than p₀ nothing happens, whereas if r₁ is larger than p₀ the vapour particle condenses onto one of the ice particles.

If condensation occurs, a second step is needed in order to decide onto which of the ice particles in the grid cell the vapour particle condenses. In this step it is no longer the absolute probability that is of interest, but the relative probabilities of the ice particles in the grid cell. This can be expressed as ${\mathit{p}}_{\mathit{i}}^{\mathrm{\ast }}\mathrm{=}\frac{\mathrm{\Delta }\mathit{t}}{{\mathit{\tau }}_{\mathit{i}}}\mathrm{·}$(40)The relative probability for each of the ice particles in the grid cell are placed in a sequence, $\mathrm{0}\mathit{,}\hspace{0.17em}{\mathit{p}}_{\mathrm{1}}^{\mathrm{\ast }}\mathit{,}\hspace{0.17em}{\mathit{p}}_{\mathrm{1}}^{\mathrm{\ast }}\mathrm{+}{\mathit{p}}_{\mathrm{2}}^{\mathrm{\ast }}\mathit{,}\hspace{0.17em}{\mathit{p}}_{\mathrm{1}}^{\mathrm{\ast }}\mathrm{+}{\mathit{p}}_{\mathrm{2}}^{\mathrm{\ast }}\mathrm{+}{\mathit{p}}_{\mathrm{3}}^{\mathrm{\ast }}\mathit{,}\hspace{0.17em}\mathit{...}\mathit{,}\sum _{\mathit{i}}{\mathit{p}}_{\mathit{i}}^{\mathrm{\ast }}\mathit{,}$(41)and a new random number is generated, ${\mathit{r}}_{\mathrm{2}}\mathrm{\in }\mathrm{\left]}\mathrm{0}\mathit{,}{\sum }_{\mathit{i}}{\mathit{p}}_{\mathit{i}}^{\mathrm{\ast }}\mathrm{\right[}$ . Which interval in the sequence r₂ falls in decides which ice particle the vapour particle condenses onto, such that r₂ falling in the interval $\mathrm{\left]}{\mathit{p}}_{\mathit{i}\mathrm{-}\mathrm{1}}^{\mathrm{\ast }}\mathit{,}{\mathit{p}}_{\mathit{i}}^{\mathrm{\ast }}\mathrm{\right]}$ implies condensation onto ice particle i.

The condensation time scale in Eq. (13) is proportional to the radius of the particle, so that a smaller particles has a shorter, and a larger particle a longer, condensation time scale. This corresponds to a larger condensation probability for smaller particles and vice versa, which may seem counter-intuitive, but is explained by the fact that the simulation handles superparticles, each representing the same mass M. Therefore a superparticle representing small physical particles represents a larger combined surface area than one representing large particles. Condensation is thus more likely to happen to a superparticle representing small particles. Looking at a single-physical-particle-basis, the result is correct averaged over many time steps and particles. A condensation event means doubling the mass of the physical particles involved, implying that a small particle experiences many condensation events, but the absolute growth in radius each event is small, whereas a single large physical particle experiences few condensation events, but with a large growth in absolute radius each event.

4. Test problems

In order to understand the results of the computer simulations and to make sure that the code is functioning correctly, two major tests of the algorithms used were made. Firstly, the dynamical behaviour of the code was tested without including particle growth. Secondly, the particle growth algorithm was tested, without taking spatial dimensions into account.

4.1. Test of the dynamical behaviour of the particles

The dynamical behaviour of the particles is tested by excluding condensation and sublimation from the simulation. This means that particles move because of turbulence that stirs them up in a turbulent diffusion, and because of gravity directed towards the midplane, but no particle growth is included. We also ignore radial drift. The particles settle to an equilibrium, where the particle scale height depends on the size of the particles, since the strength of the coupling to the turbulent gas is a function of particle size. Figure 2 shows the particle scale height relative to the gas scale height H_p/H as a function of particle size R. The particles settle into a Gaussian distribution around the midplane, so the particle scale height in equilibrium can be retrieved as the root mean square of the particle positions, ${\left(\frac{{\mathit{H}}_{\mathrm{p}}}{\mathit{H}}\right)}_{\mathrm{mod}}\mathrm{=}\sqrt{\frac{\mathrm{1}}{\mathit{N}}\mathrm{(}{\mathit{z}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{+}\mathit{...}\mathrm{+}{{\mathit{z}}_{\mathit{n}}^{\mathrm{2}}}^{\mathrm{)}}}\mathit{,}$(42)where n is the number of particles and z_i is the vertical position of particles. This is compared to the theoretically expected particle scale height, given by an equilibrium between sedimentation and turbulent diffusion, ${\left(\frac{{\mathit{H}}_{\mathrm{p}}}{\mathit{H}}\right)}_{\mathrm{theo}}\mathrm{=}\sqrt{\frac{\mathit{\alpha }}{\mathrm{\Omega }{\mathit{\tau }}_{\mathit{f}}\mathrm{+}\mathit{\alpha }}}\mathrm{·}$(43)Figure 2 shows how the particle scale height depends on the size of the particles, with the modelled values represented by the black full line and the expected analytical values as a red dotted line. The modelled line clearly follows the analytical curve, showing that the dynamical behaviour of the particles is as expected. The particle size is given in units of R₁, where R₁ ≈ 1 m at r = 3 AU. From the figure it can be seen that small particles are fully coupled to the gas, as H_p/H ≈ 1 for particles of R ≲ 10^-3 R₁, a size that corresponds to mm-sized ice particles near the water ice line. Furthermore, a slight change of the slope can be seen at R ≈ 0.5 R₁. This is caused by the change from Epstein to Stokes drag regime, which makes the particles decouple more quickly from the gas and hence increases sedimentation.

4.2. Test of the condensation algorithm

We test the condensation algorithm by letting particles grow via condensation, without including the spatial dimensions. The ice line is therefore not included, so sublimation is not taken into account. The number of ice superparticles is known, as is the number of vapour superparticles. The total mass available for growth is known since we run the simulation until all vapour has condensed onto the ice particles. By tracking the initial and final mass of all superparticles the growth in radius is followed, and can be compared to a theoretical expectation.

In reality, condensation is a continuous process on macroscopic scales, where single water molecules condense onto ice particles so that growth happens little by little all the time. In the Monte Carlo scheme used in this code this continuous behaviour is modelled by a discrete growth process, where a condensation event involves one ice superparticle and one vapour superparticle. The mass of the physical ice particles represented by the ice superparticle doubles as the vapour superparticle completely condenses onto the ice, implying that there is no vapour superparticle left after the event. To conserve the number of ice particles and superparticles (for easier coding purposes) this implies that the physical ice particles that before the event were represented by the one ice superparticle involved in the event, now are represented by two superparticles, i.e. the number of physical ice particles before the event are after the event split between two ice superparticles. The mass of the physical particles involved in an event thus doubles when a condensation event occurs. Since particle radius is connected to mass as m ∝ R³, for any number of condensation events n_ce the radius increases as $\mathit{R}\mathrm{\to }{\left({\mathrm{2}}^{{\mathit{n}}_{\mathrm{ce}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{R}\mathit{.}$(44)The final radius R_final is plotted as a function of R_init as black crosses in Fig. 3. In this test 2000 ice superparticles were distributed evenly in 20 different initial size bins over the range R/R₁ =  [0,1], and 20 000 vapour superparticles were added. The possible growth in radius from mass doubling events is plotted as red lines in Fig. 3, for a different number of condensation events 0 ≤ n_ce ≤ 10. From the figure it is clear that all modelled values indeed fall on these lines, and never in-between, as expected because of the discrete nature of the Monte Carlo method used.

Fig. 2 Comparison between modelled and analytical particle scale height as a function of particle size, for a turbulence strength of α = 10-2. The black full line denotes modelled values and the red dotted line analytical values. Small particles R ≲ 10-3 are perfectly coupled to the gas, whereas larger particles sediment towards the midplane. The slope change at R ≈ 0.5 is caused by the change from Epstein to Stokes drag regime. The particle size is given in units of R1, where R1 ≈ 1 m at r = 3 AU. Overall, modelled and analytical values are in good agreement.

The modelled particle growth can be compared to the analytically expected ΔR. For a physical ice particle i the total mass growth due to condensation during the time Δt is $\mathrm{\Delta }{\mathit{m}}_{\mathit{i}}\mathrm{=}{\mathit{m}}_{\mathit{i}}\mathrm{-}{\mathit{m}}_{\mathrm{0}\mathit{,i}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}\mathrm{(}{\mathit{R}}_{\mathit{i}}^{\mathrm{3}}\mathrm{-}{{\mathit{R}}_{\mathrm{0}\mathit{,i}}^{\mathrm{3}}}^{\mathrm{)}}{\mathit{\rho }}_{\mathrm{•}}\mathit{,}$(45)where ρ_• is the material density of ice, R_0,i is the initial and R_i the final radius of the ice particle. However, the physical particles are represented by superparticles in the code, and the mass growth must therefore be given for superparticles and not physical particles. One superparticle represents the mass M_i = n_im_i, where n_i is the number of physical particles represented by the superparticles.

The total mass of vapour that condenses onto the ice particles is given by the difference between the sum of final and initial ice superparticles, ${\mathit{M}}_{\mathrm{v}}\mathrm{=}{\mathit{\rho }}_{\mathrm{v}}\mathit{V}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}{\mathit{\rho }}_{\mathrm{•}}\left(\sum _{\mathit{i}}{\mathit{n}}_{\mathit{i}}{\mathit{R}}_{\mathit{i}}^{\mathrm{3}}\mathrm{-}\sum _{\mathit{i}}{\mathit{n}}_{\mathrm{0}\mathit{,i}}{\mathit{R}}_{\mathrm{0}\mathit{,i}}^{\mathrm{3}}\right)\mathit{.}$(46)The total number of physical ice particles is constant, as is the total number of ice and vapour superparticles, but the number of ice superparticles is not. While the physical particles represented by one superparticle grow in mass, the number of physical particles represented by the superparticle decreases since the total mass of a superparticle is constant. In order to conserve the number of physical particles, total number of superparticles and the total mass of each superparticle, the number of ice superparticles increases as the physical particles grow in mass, so that each superparticle represents fewer and fewer physical particles, which is why the two sums in Eq. (46) are over different numbers of superparticles.

Fig. 3 Comparison between average modelled growth and analytically expected growth of ice particles. The final particle size as a function of initial particle size is shown for 2000 ice superparticles distributed evenly in 20 different initial size bins over the range R/R1 =  [0,1] , and with 20 000 vapour superparticles added. Black crosses show modelled values and blue stars show the average growth in each size bin. The modelled average size roughly agrees with the analytical value, shown as a blue line. Red lines denote expected size for a number of condensation events between 0 (lower line) and 10 (upper line).

The right handside of Eq. (46) can be rewritten using ${\mathit{n}}_{\mathit{i}}\mathrm{=}\frac{{\mathit{M}}_{\mathit{i}}}{{\mathit{m}}_{\mathit{i}}}\mathrm{=}\frac{{\mathit{\rho }}_{\mathit{i}}\mathit{V}}{{\mathit{\rho }}_{\mathrm{•}}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}{\mathit{R}}_{\mathit{i}}^{\mathrm{3}}}\mathrm{·}$(47)Here V is the volume represented by a superparticle, which in this one-grid-cell test problem is equal to the total volume of the box. By cancelling terms and gathering the densities on the left-hand side Eq. (46) can be rewritten as $\frac{{\mathit{\rho }}_{\mathrm{v}}}{{\mathit{\rho }}_{\mathit{i}}}\mathrm{=}{\mathit{N}}_{\mathrm{p}}\mathrm{-}{\mathit{N}}_{\mathrm{p}\mathit{,}\mathrm{0}}\hspace{0.17em}$(48)where N_p and N_p,0 is the final and initial number of ice superparticles, respectively. All superparticles represent the same density, and therefore the ratio between the two densities corresponds to the ratio between the number of vapour and ice superparticles represented, so that $\frac{{\mathit{\rho }}_{\mathrm{v}}}{{\mathit{\rho }}_{\mathit{i}}}\mathrm{=}\frac{{\mathit{N}}_{\mathrm{v}}}{\mathrm{1}}\mathit{,}$(49)with N_v being the number of vapour superparticles added to the simulation. Equation (46) is therefore equivalent to ${\mathit{N}}_{\mathit{v}}\mathrm{=}{\mathit{N}}_{\mathrm{p}}\mathrm{-}{\mathit{N}}_{\mathrm{p}\mathit{,}\mathrm{0}}\mathit{,}$(50)which the code can be checked against, since both the number of vapour superparticles put into the system and the initial number of ice superparticles is specified, and the final number of superparticles is given as an output from the program. This analysis shows that the code is mass-conserving, by construction.

From Eq. (46) the theoretically expected ΔR, equal for all particles, can be found. The final radius can be rewritten as a_i = ΔR + R_0,i, where R_0,i is the initial radius of the particle. The two sums are not necessarily over the same number of superparticles, as each condensation event introduces a new ice superparticle in the system, at the expense of the vapour particle representing the condensing vapour. The physical particles originally represented by the first ice superparticle are now split between the new and the old ice superparticle. The initial size of a physical particle represented by a converted vapour superparticle can therefore be found as the initial size represented by the first ice superparticle. Equation (46) is thereby rewritten as ${\mathit{M}}_{\mathrm{v}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}{\mathit{\rho }}_{\mathrm{•}}\left[\sum _{\mathit{i}}{\mathit{n}}_{\mathit{i}}\mathrm{(}\mathrm{\Delta }\mathit{R}\mathrm{+}{\mathit{R}}_{\mathrm{0}\mathit{,i}}{\mathrm{)}}^{\mathrm{3}}\mathrm{-}\sum _{\mathrm{0}\mathit{,i}}{\mathit{n}}_{\mathit{i}}{\mathit{R}}_{\mathrm{0}\mathit{,i}}^{\mathrm{3}}\right]\mathit{,}$(51)which, by using Eq. (47), cancelling terms and noting that (as in Eq. (46)) ρ_v/ρ_i = N_v, gives ${\mathit{N}}_{\mathrm{v}}\mathrm{=}\sum _{\mathit{i}}\frac{\mathrm{1}}{{\mathit{R}}_{\mathit{i}}^{\mathrm{3}}}\mathrm{(}\mathrm{\Delta }\mathit{R}\mathrm{+}{\mathit{R}}_{\mathrm{0}\mathit{,i}}{\mathrm{)}}^{\mathrm{3}}\mathrm{-}{\mathit{N}}_{\mathrm{p}\mathit{,}\mathrm{0}}\mathit{.}$(52)The initial number of vapour and ice superparticles N_v and N_p,0 are known, as they are input parameters to the simulation. The radius represented by a superparticle in the end of the simulation, a_i, and when it is created, R_c,i, are both given as output from the code. The radius growth for each particle ΔR can therefore be found. For N_v = 20 000 and N_p,0 = 2000, a value of ΔR ≈ 0.202 is found. In Fig. 3 the blue full line represents the expected growth of the 2000 initial ice superparticles, (R_0,i + Δa), as a function of initial radius R_0,i. Comparing the expected final particle radii with the actual final radii shows that for small particle sizes the Monte Carlo method works fairly well, whereas for larger particle sizes most particles experience too little (zero) growth and a few grow several orders of magnitude more than expected. This means that one must be careful not to draw conclusions from the growth of individual particles. However, statistically speaking, i.e. by averaging over many particles, the result can be considered to be correct (Zsom & Dullemond 2008).

5. Results

Fig. 4 Pressure ice line for a vertically isothermal disc. The full black line denotes the ice line, whereas the coloured dotted lines show the location of 1 H, 2 H and 3 H. As particles are preferably found within a few scale heights from the midplane, growth due to crossing of the atmospheric ice line is only important within less than 1 AU of the radial ice line.

The results are divided into two parts: runs including only the atmospheric ice line (Table 1: runs 1 and 2) and runs including both the radial and atmospheric ice lines (Table 1: runs 3 and 4). We assess how the growth depends on turbulence strength in runs 2 and 4, its dependence on atmospheric ice line position in runs 1, 2 and 3 and the dependence on box size in run 4. In all runs particle collisions are ignored, an assumption that we discuss in Sect. 6.4.

Fig. 5 State of simulation with only the atmospheric ice line included, for t = 0, 200 Ω and 2000 Ω-1, from left to right. Upper panels: distribution of blue ice particles and red vapour particles, with the grey area representing the condensation zone. The size of the blue dots are proportional to the size of the ice particles and the number of particles shown is inversely scaled with size for visibility. Lower panels: evolution of the size distribution of ice particles at corresponding times. N is the number of ice superparticles, and the size is given as R/R1. At the water ice line, R1 ≈ 1.3 m.

5.1. Atmospheric ice line

We start by exploring the atmospheric ice line, not including the radial ice line. This means that the simulation domain is located just outside of the radial ice line, with a colder midplane and hotter outer layers (Dullemond & Dominik 2004).

The thermal atmospheric ice line in a typical protoplanetary disc is located at z ≳ 3 H (Chiang & Goldreich 1997). However, the decrease of pressure with height above the midplane leads to a narrow radial zone just outside of the radial ice line where the atmospheric ice line is located at z ≤ 3 H. The pressure ice line in a vertically isothermal disc is shown in Fig. 4.

The height of the atmospheric ice line above the midplane z_ice is varied from z_ice = 0.4 H to z_ice = 3 H. The lower limit is the lowest value for which not all ice particles immediately sublimate. A lower value gives a distance between the midplane and the ice line that is comparable to the length a particle moves during one time step (see Eq. (17)), making it possible for all ice particles to escape from the condensation region before any growth has taken place. Periodic boundary conditions are used in both the radial and vertical direction.

Figure 5 shows the time evolution of a simulation with the atmospheric ice line at z_ice = 1 H. The distribution of red vapour and blue ice particles is shown together with the size distribution of ice particles, for t = 0, t = 200 Ω^-1 and t = 2000 Ω^-1. Vapour diffuse into the grey condensation region, condensing onto already existing ice particles. As small ice particles tend to diffuse out of the condensation region and sublimate, whereas larger ones stay in the midplane and experience continued growth, the result is a narrow size distribution with decimetre-sized ice pebbles residing in a midplane layer.

5.1.1. Varying the atmospheric ice line position

Fig. 6 Mean ice particle size as a function of time for different positions of the atmospheric ice line. The different lines show the growth for different distances from the atmospheric ice line to the midplane, zice. From top to bottom zice =  0.4 H, 0.6 H, 1.0 H, 1.4 H,1.8 H and 3.0 H. Ice particles grow faster, and to larger sizes, the closer to the midplane the atmospheric ice line is. Full lines denote simulations with 10 000 particles, dashed lines 1000 particles and dotted lines 100 000 particles.

Fig. 7 Mean ice particle size as a function of the position of the atmospheric ice line and the strength of the turbulence. The mean particle size ⟨ R ⟩ is shown as a function of the distance from the atmospheric ice line to the midplane zice, in gas scale heights H. The system is shown when the particle sizes have approximately reached an equilibrium, at t = 5000 Ω-1 for α = 10-2, shown in blue, at t = 25 000 Ω-1 for α = 10-3, shown in red, and at t = 100 000 Ω-1 for α = 10-4, shown in yellow. The initial particle size is 10-4 R1, corresponding to approximately 0.1mm at the ice line. Solid lines are modelled values and dotted lines are analytical values, set by Eq. (53). The modelled and analytical values follow the same slope in the applicable range.

Simulations were run with the atmospheric ice line placed at different heights z_ice above the midplane. Decreasing z_ice corresponds in a physical protoplanetary disc to placing the simulation box radially closer in towards the radial ice line, where a hypothetical z_ice = 0 would correspond to the position of the radial ice line, and thus assessing the narrow radial zone where the pressure ice line dominates.

In Fig. 6 particle growth for different z_ice is shown. The mean ice particle size ⟨ R ⟩ in units of R₁ is shown as a function of time in units of Ω^-1. As is clear from the figure, particles grow to larger sizes as the ice line is moved closer to the midplane. This is effectively due to the fact that for ice lines closer to the midplane, i.e. a narrower condensation region, vapour particles can easier penetrate to the larger ice particles residing in the midplane. Figure 2 shows how the particle scale height H_p decreases with increasing particle radius R. Small particles have a scale height comparable to the gas scale height H_p,small ≈ H, whereas for larger particles H_p,large ≪ H. Placing the ice line at a height larger than 1 H is thus equivalent to saying that even the smallest particles tend to stay within the condensation zone. There is therefore very little material available for particle growth. The few particles that do move across the ice line and sublimate will mostly redistribute material among the particles at the border of the condensation zone. Growth in runs with z_ice ≳ 1 H is both slow, and to modest particle sizes. For ice line positions of z_ice ≥ 3 H the particle growth is of the order of ΔR ~ 10^-4 R₁ or less in 1000 Ω^-1. For the larger extent of the disc, where z_ice ≳ 3H, the atmospheric ice line is thus of little importance for particle growth by condensation.

In Fig. 6 we show runs with 1000, 10 000 and 100 000 particles for each ice line position. The negligible difference between runs with 10 000 and 100 000 particles in combination with a much lower computational cost motivates the use of 10 000 particles as a reference resolution.

5.1.2. Varying the turbulent α-value

Simulations were run with different values of α in order to test the influence of turbulence strength on the results. In Fig. 7 the mean particle size ⟨ R ⟩ is shown as a function of the height of the ice line above the midplane z_ice, starting from an initial particle size of 10^-4  R₁. The system is shown when it has approximately reached an equilibrium, at t = 5000 Ω^-1 for α = 10^-2, shown in blue, at t = 25 000 Ω^-1 for α = 10^-3, shown in red, and at t = 100 000 Ω^-1 for α = 10^-4, shown in yellow. Full lines represent modelled values and dotted lines are analytical estimates. As can be seen, lowering the strength of turbulence gives less growth. This is to be expected since the main effect of decreasing the level of turbulence is to lower the effective particle scale height, set by an equilibrium between sedimentation and turbulent diffusion. As a comparison, dotted lines show the analytical relation between mean particle size and particle scale height in a system where vertical gravity and turbulent diffusion is dominating. The scale height of particles in terms of the gas scale height is $\frac{{\mathit{H}}_{\mathrm{p}}}{\mathit{H}}\mathrm{=}\sqrt{\frac{\mathit{\alpha }}{\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}\mathrm{+}\mathit{\alpha }}}\mathrm{·}$(53)Setting H_p ~ (1/c) z_ice gives an analytical estimate of the maximum particle size achievable by condensation, shown in Fig. 7, as $\mathrm{\Omega }{\mathit{\tau }}_{\mathrm{f}}\mathrm{=}\mathit{\alpha }\left\{\frac{\mathrm{1}}{\mathrm{[}{\mathit{z}}_{\mathrm{ice}}\mathit{/}\mathrm{(}\mathit{c}\hspace{0.17em}\mathit{H}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}}\mathrm{-}\mathrm{1}\right\}\mathit{.}$(54)We chose c = 3, but even in this case the analytical expression underestimates the resulting R slightly. Nevertheless our simulations show that particles grow approximately to a size where their scale height is 1/3 of the ice line height.

Fig. 8 Mean ice particle size ⟨ R ⟩ as a function of time in inverse Keplerian orbits Ω-1 for simulations including both radial and atmospheric ice lines. The height of the atmospheric ice line above the midplane zice is denoted as: black 0.4 H, violet 0.6 H, blue 1.0 H, green 1.4 H, red 1.8 H and yellow 3.0 H. For the runs with zice = 3.0 H, the dashed line corresponds to 1000 particles and the dotted line to 100 000 particles. Full lines denote runs with 10 000 particles. The position of the radial ice line is rice =  −0.3 H, and the strength of turbulence is α = 10-2.

5.2. Atmospheric and radial ice line

Including both the atmospheric and radial ice lines is done by letting the simulation domain represent a region around the radial ice line. This means that ice particles can sublimate both by moving away from the midplane, and by moving closer to the central star. For a realistic atmospheric ice line position at z_ice ≳ 3 H, sublimation by moving across the radial ice line is however much more important than by moving across the atmospheric ice line. The simulation domain is 12 H in the vertical direction. In the radial direction the box size is varied from 1.5 H to 5.25 H, with the number of particles and grid cells changed accordingly in order to keep the particle density and the effective resolution fixed. Periodic boundary conditions are used in the vertical direction, and reflective boundary conditions in the radial direction. In the vertical direction there is in effect no difference between + z and − z, because the conditions above the midplane mirror those below it. In the radial direction, when including the radial ice line, there is no such symmetry. At r = r_min conditions such that ice sublimates prevail at all z, whereas at r = r_max the midplane is part of the condensation zone, where ice particles can exist without sublimating. Using periodic boundary conditions in the radial direction would therefore artificially introduce vapour in the system and thereby cause a too large particle growth.

5.2.1. Varying the atmospheric ice line position

Including the radial ice line greatly increases the growth efficiency. As shown in Sect. 5.1, models with only the atmospheric ice line gives negligible growth for ice line positions of z_ice ≳ 3 H. However, also taking the radial ice line into account leads to particle growth beyond 0.1 R₁, corresponding to decimetre-sized pebbles at the ice line, within 1000 Ω^-1.

Figure 8 shows how the mean particle size ⟨ R ⟩ as a function of time varies with atmospheric ice line position z_ice when the radial ice line is included. Growth is somewhat faster for z_ice close to the midplane and slower for z_ice farther away from it. As previously discussed, this is due to the larger supply of material available for growth in a narrow condensation region, as compared to a wider condensation region. For ice line positions above 3 H, a growth similar to that for z_ice = 3.0 H is found.

Growth stops at ⟨ R ⟩  ≈ R₁ for all z_ice. This is due to radial drift towards the star that eventually removes all ice particles from the simulation. From Eq. (30) it can be seen that the radial drift inwards peaks for particles of ⟨ R ⟩  ≈ R₁. When particles reach this size they thus drift inwards, past the radial ice line, and sublimate.

We also show the growth for ice line position of z_ice = 3.0 H for a smaller (n_p = 1000) and a larger (n_p = 100 000) number of particles in Fig. 8, in addition to n_p = 10 000 particles. For particle sizes R ≲ 0.1 R₁ we find converging results, whereas for larger particle sizes the growth is dominated by statistical fluctuations caused by runaway growth of a few lucky particles. Our condensation model is constructed to give correct results in a regime where many particles compete for vapour, i.e. up to pebbles of a few decimetres, but may be unphysical in the regime of metre-sized boulders and beyond.

Fig. 9 State of an extended simulation box with both the radial and the atmospheric ice line included, for t = 0, 100 and 1000 Ω-1, from top to bottom. The grey area represents the condensation zone, and ice and vapour is shown in blue and red, respectively. The sizes of the blue dots are proportional to the size of the ice particles and the number of particles shown is inversely scaled with size for visibility. The turbulent α-value is 10-2.

5.2.2. Varying the radial extent of the simulation and the turbulent α-value

To assess the importance of the size of the simulation domain, we ran simulations where the box size was varied in the radial direction. Figure 9 shows the time evolution of a run where the radial extent of the simulation box has been increased, so that r_max = 5.25 H. The most significant growth takes place close to the radial ice line, but owing to radial mixing large particles can be found throughout the entire radial extent of the simulation box. In Fig. 10 the partial growth in four equally sized subdomains can be compared to the total growth in the whole simulation domain, as shown in Fig. 9. This illustrates that the total particle growth is dominated by the growth near to the radial ice line.

The left panel of Fig. 11 shows the particle growth for different radial extents of simulation domains for a strongly turbulent disc, α = 10^-2. The growth time scale increases with larger simulation domains, as the amount of ice particles is increased while the supply of water vapour across the ice line remains unchanged. Results from simulations with weak turbulence are presented in the middle (α = 10^-3) and right (α = 10^-4) panels of Fig. 11. These low-turbulence simulations show that particles can grow via condensation also in dead zones, although slowly.

6. Discussion

Our results show that ice condensation is an efficient way to form centimetre-sized and decimetre-sized pebbles. Ice condensation does not suffer from bouncing and fragmentation barriers and could be the dominant mode of growth in protoplanetary discs, to sizes where particles can concentrate strongly in the turbulent gas and continue towards planetesimal sizes by gravitational collapse. However, our results rely on a number of assumptions and simplifications that we discuss here.

6.1. Other ice lines

This model has been developed primarily to investigate the water ice line at r ≈ 3 AU. Similar ice lines, or condensation fronts, exist also for other chemical species, and as this model is scalable to any r it can easily be adapted to include any condensation front. Of importance for planet formation purposes are in particular the ammonia, methane, carbon monoxide and molecular nitrogen ice lines farther out in the protoplanetary disc, and the numerous silicate condensation fronts farther in towards the central star. A global disc model, which includes all condensation fronts of the most important chemical species in the protoplanetary disc, will be an important extension of this work.

6.2. Disc evolution

The systematic motion of gas accretion onto the central star has been ignored in this work. However, compared to the time scales of particle growth by condensation the accretion process is slow enough to be safely neglected (Hartmann et al. 1998).

The ice line has throughout this work been considered as fixed. Seen over the life time of the disc this is not true, since the ice line probably moves both in a systematic way over long time scales and because of shorter heating events (Armitage 2011). However, as the time scale for growth by condensation found in this work is of the order of τ_c ≈ 1000 Ω^-1, which is very short compared to the life time of the disc τ_disc ~ 3 Myr, the assumption of a fixed ice line is reasonable.

We also ignore the potential effect the processes of condensation and sublimation have on the ice line position, since the resulting release and absorption of energy is relatively small compared to the total thermal energy of the disc (Stevenson & Lunine 1988).

Fig. 10 Growth in different radial subdomains of the simulation box shown in Fig. 9. Black shows the average particle size in the total box, whereas yellow, red, green and blue give the average particle sizes in the different subdomains, from left to right. The total particle size at any time is dominated by the growth near the radial ice line.

Fig. 11 Particle growth for different radial extents of the simulation box and turbulence strengths. The mean particle size ⟨ R ⟩ is shown as a function of time in inverse Keplerian orbits Ω-1 for α = 10-2 (left panel), α = 10-3 (middle panel) and α = 10-4 (right panel). The size of the box has been extended in the radial direction, with rmax colour coded as: black 0.75 H, blue 2.25 H, red 3.75 H and yellow 5.25 H. For all runs zmin =  −6.0 H, zmax = 6.0 H and rmin =  −0.75 H. The inlay figure shows the simulation domain, with the considered radial extents marked.

6.3. Ice nuclei

We have modelled ice particles as one-component particles; homogeneous water ice particles. In reality these particles would have (at least) a two-component composition, that of a rocky core and an icy mantle. This is due to the fact that supersaturated water vapour is needed for homogeneous nucleation, i.e. for the spontaneous formation of a new ice particle. As this is typically not the case in the conditions prevailing in a protoplanetary disc, water vapour instead condenses onto an existing grain, either a rocky dust particle, or a particle with an already existing ice layer on top of the dust grain. Whether vapour preferably condenses onto a bare dust grain or an ice particle depends on material properties and on their respective sizes.

Ignoring the small refractory dust particles present in the disc possibly leads to an over-estimation of the growth of the larger ice particles. For a fixed amount of mass, smaller dust particles have a larger combined surface area than larger particles, and therefore there is a higher probability that a vapour particle condenses onto a smaller dust superparticle than onto a larger superparticle. This is for the case when the material properties of dust and ice are such that two equally large dust and ice particles have equal probabilities of having a vapour particle condensing onto them. Instead of growth of a few large ice particles that sediment towards the midplane and thereby are protected against sublimation, there would thus be a continuous sublimation and condensation of the small ice particles at the ice line. This effect can already be seen in our condensation model, as small particles are more likely to cross the ice line and sublimate, and conversely, to grow by condensation, but could be amplified by the introduction of refractory grains.

Material effects might make this problem more severe. The dust grain composition is assumed to be similar to that of interstellar dust grains, either being pure interstellar dust grains or conglomerates thereof. This gives us the most important dust species as silicates and carbonaceous material, with silicates typically assumed to be the dominating one (Draine 2003). Silicates have material properties that makes them more efficient as ice nuclei than pure water ice particles (Goumans et al. 2009), implying that the effect of small particle growth at the expense of the larger ones is amplified when taking the material properties into account.

There are however a number of possible solutions to this problem, in terms of reasons for vapour to condense onto ice particles instead of the small dust grains. Firstly, the dust grains are not likely to all have the same size. Rather, they would have a size distribution where the smallest grains can be nanometre sized (Draine 2003). The largest dust particles of relevance are the ones of the largest size that is strongly enough coupled to the gas to be present at the relevant vertical distance from the midplane. This is around micrometre-scale (see Fig. 3). The equilibrium vapour density for a curved surface is higher than for a flat surface, so that vapour more easily condenses onto a flat surface than a curved one. This is typically referred to as the Kelvin curvature effect, and is most important for the smallest particles, around nanometre sizes (Sirono 2011). Therefore vapour will effectively prefer condensing onto the larger particles, leaving the smallest ones as bare dust grains.

Secondly, small grains are likely to be hotter than larger particles, moving the ice line of small particles closer to the midplane than the ice line of large particles. The temperature of the grains is set by the balance between absorption and emission efficiency. Since grains absorb and emit radiation most efficiently in wavelengths smaller than and up to their own size, there is a size-dependent effect (Krivov et al. 2008; Vitense et al. 2012). The spectral energy distribution of a solar-type star peaks at wavelengths of about 500 nm. All grains of this size and larger therefore absorbs incoming radiation with approximately equal efficiency. However, grains emit in infra-red wavelengths (larger than 1 μm), and thus the larger grains can emit more energy than the smaller ones. The resulting temperature for the smaller grains is therefore higher than for the larger grains. Also the material affects the temperature of the grains. As water ice is nearly transparent in visible light, where a solar-type star peaks, an ice particle is heated less by the stellar radiation than a particle of a more absorbing material (Lecar et al. 2006). Both mechanisms are valid for grains residing in the atmosphere of the disc, where the disc is optically thin to the stellar radiation, but not for grains in the midplane, where the received stellar radiation is mainly the re-emitted infrared radiation from neighbouring grains.

Furthermore, there is a possibility that material effects might prevent water vapour from condensing onto the dust grains. As mentioned above, bare silicate grains are more efficient as ice nuclei than ice-covered ones. However, for carbonaceous grains the opposite is true (Papoular 2005), meaning that vapour condenses more easily onto a water ice particle than onto a bare carbonaceous grain. For a grain population with ice coated particles and bare carbonaceous grains, the carbonaceous grain would therefore be left bare, whereas the ice grains would continue growing.

6.4. Particle collisions

We have neglected collisions between particles throughout this work. Depending on material properties and relative velocities of the particles, collisions can lead to increased growth via coagulation, to fragmentation or to bouncing. In general, small particles tend to stick together to a larger extent than large particles, if colliding, whereas larger particles are more prone to bouncing or fragmenting. Collisions amongst large particles therefore have a tendency to be destructive, whereas small-particle collisions are more likely to favour growth (or at least not counteract it). For silicate particles growth via collisions involving equal-sized particles is very difficult beyond millimetres, however collisional growth involving small particles colliding with a large target is possible, although slow (Wurm et al. 2005; Johansen et al. 2008; Windmark et al. 2012).

Whereas extensive experimental data exists on collisions between rocky particles, the outcome of ice particle collisions is not equally well-known. Bridges et al. (1996) performed low-velocity collisions (v_rel < 5 cm s^-1) between ice pebbles and found that ice particles covered with a frost layer have an increased stickiness compared to rocky particles. A mechanism of collisional fusion, in which ice particles undergo a phase change during collision, has been suggested as a way for particles to stick also in collisions with higher velocities (1  m s^-1 < v_rel < 100 m s^-1) (Wettlaufer 2010).

Collision experiment for higher velocities have been used to derive criticial relative velocities above which fragmentation can occur, v_crit. Higa et al. (1996) found a critical velocity of v_rel ≈ 1 m s^-1, but other groups have found significantly larger values (Arakawa 1999; Arakawa et al. 2002). Computer simulations suggest critical velocities of up to v_rel ≈ 100 m s^-1 (Benz & Asphaug 1999).

Collision velocities for millimetre-sized particles are set by the turbulent velocities and are expected to be of the order of a few m s^-1, whereas collisions involving larger particles approach the radial drift velocity, v_rel ≈ 50 m s^-1 (Brauer et al. 2008). Both coagulation and fragmentation are therefore expected to occur as a result of collisions. A likely scenario is that collisions between large particles lead to fragmentation, where the resulting small ice fragments are later swept up in collisions with other particles.

An important implication of collisions is that it provides a natural means of removing small dust grains released from the ice particles when sublimating. As these dust grains are very efficient ice nuclei, they might prevent growth of already large particles by “stealing” all the vapour, as discussed in Sect. 6.3. However, the small dust grains are likely to be swept up by larger ice particles in collisions. This does not increase the growth significantly, but has the important benefit that the small dust grains are removed so that vapour has to condense onto growing ice particles.

We plan to include the effects of multiple condensable species, refractory ice nuclei and particle collisions in future work. The present work highlights that simple turbulent dynamics can cause significant particle growth by condensation of volatiles, motivating follow-up studies with increasingly realistic models for the condensation and collision processes.

7. Conclusions

In this paper we have shown that ice condensation is an important particle growth mechanism that must to be taken into account in models of early planet formation. As the more thoroughly investigated mechanism of coagulation is inefficient in forming particles larger than millimetres, growth by condensation is an important mechanism that could complement coagulation or even be the dominant particle growth mode.

Our results show that growth by condensation near ice lines is rapid and results in large particle sizes. The growth time scale from millimetre-sized dust grains to decimetre-sized pebbles is τ_c ≈ 1000 years. Significant growth is obtained for a range of turbulent α-values from 10^-4 to 10^-2, where the higher value corresponds to the strength of turbulence that has been inferred from observations of protoplanetary discs, and the lower value has been suggested in a dead zone with weak turbulence stirred by the narrow active layers.

An implication of our model is a lower column density of vapour in the entire region outside the radial ice line compared to the inner part of the disc, caused by condensation onto existing grains and subsequent sedimentation of particles. This is in agreement with observations of CO (Qi et al. 2011), where a drop in abundance was found outside of the CO ice line. The far outer disc regions have also been observationally inferred to be depleted in water vapour (Hogerheijde et al. 2011). Observed remnant water and CO vapour in the outer disc shows the importance of the coexistence of a cold midplane and an upper atmospheric layer where ice can sublimate or photodesorb into vapour.

Ice condensation can also explain observations of large quantities of pebbles in protoplanetary discs (Wilner et al. 2005). Such pebbles are crucial in planet formation models, as they are the preferred starting size for planetesimal formation by clumping via streaming instabilities, in pressure bumps and in vortices, followed by gravitational collapse. Once large planetesimals have formed, continued pebble accretion is a very efficient path to formation of the cores of gas and ice giants (Ormel & Klahr 2010; Lambrechts & Johansen 2012; Morbidelli & Nesvorny 2012).

Ice condensation is a natural consequence of temperature gradients in protoplanetary discs. With our simulations we have shown that condensation is an efficient growth mechanism with the potential to explain the formation of decimetre-sized pebbles in protoplanetary discs, thereby providing the missing link to further growth into planetesimals and planets.

Acknowledgments

We thank Jürgen Blum, Melvyn B. Davies, Ewine van Dishoeck, Kees Dullemond, Michiel Lambrechts, Klaus Pontoppidan, Erik Swietlicki, John Wettlaufer and Andrew Youdin for discussions and helpful suggestions. We also thank the referee, Phil Armitage, for comments that helped improve the quality of the paper. This research was partially funded by the European Research Council under ERC Starting Grant agreement 278675-PEBBLE2PLANET.

References

All Tables

All Figures

Fig. 1 Sketch of sublimation and condensation for superparticles. Blue represents ice and red vapour. Left panel: ice superparticle, representing a mass M in the form of 4 physical ice particles, each of mass m, sublimates and turns into a vapour superparticle of mass M. Right panel: vapour superparticle of mass M condenses onto an ice superparticle of mass M, representing the mass 4 m. The physical ice particles are after the event represented by two ice superparticles, each representing a total mass M, but now in the form of 2 physical ice particles each with mass 2 m. Both mass and the number of superparticles are conserved between t0 and t1. In the case of condensation, also the number of physical particles is conserved.

In the text

Fig. 2 Comparison between modelled and analytical particle scale height as a function of particle size, for a turbulence strength of α = 10-2. The black full line denotes modelled values and the red dotted line analytical values. Small particles R ≲ 10-3 are perfectly coupled to the gas, whereas larger particles sediment towards the midplane. The slope change at R ≈ 0.5 is caused by the change from Epstein to Stokes drag regime. The particle size is given in units of R1, where R1 ≈ 1 m at r = 3 AU. Overall, modelled and analytical values are in good agreement.

In the text

Fig. 3 Comparison between average modelled growth and analytically expected growth of ice particles. The final particle size as a function of initial particle size is shown for 2000 ice superparticles distributed evenly in 20 different initial size bins over the range R/R1 =  [0,1] , and with 20 000 vapour superparticles added. Black crosses show modelled values and blue stars show the average growth in each size bin. The modelled average size roughly agrees with the analytical value, shown as a blue line. Red lines denote expected size for a number of condensation events between 0 (lower line) and 10 (upper line).

In the text

	Fig. 4 Pressure ice line for a vertically isothermal disc. The full black line denotes the ice line, whereas the coloured dotted lines show the location of 1 H, 2 H and 3 H. As particles are preferably found within a few scale heights from the midplane, growth due to crossing of the atmospheric ice line is only important within less than 1 AU of the radial ice line.
In the text

Fig. 5 State of simulation with only the atmospheric ice line included, for t = 0, 200 Ω and 2000 Ω-1, from left to right. Upper panels: distribution of blue ice particles and red vapour particles, with the grey area representing the condensation zone. The size of the blue dots are proportional to the size of the ice particles and the number of particles shown is inversely scaled with size for visibility. Lower panels: evolution of the size distribution of ice particles at corresponding times. N is the number of ice superparticles, and the size is given as R/R1. At the water ice line, R1 ≈ 1.3 m.

In the text

Fig. 6 Mean ice particle size as a function of time for different positions of the atmospheric ice line. The different lines show the growth for different distances from the atmospheric ice line to the midplane, zice. From top to bottom zice =  0.4 H, 0.6 H, 1.0 H, 1.4 H,1.8 H and 3.0 H. Ice particles grow faster, and to larger sizes, the closer to the midplane the atmospheric ice line is. Full lines denote simulations with 10 000 particles, dashed lines 1000 particles and dotted lines 100 000 particles.

In the text

Fig. 7 Mean ice particle size as a function of the position of the atmospheric ice line and the strength of the turbulence. The mean particle size ⟨ R ⟩ is shown as a function of the distance from the atmospheric ice line to the midplane zice, in gas scale heights H. The system is shown when the particle sizes have approximately reached an equilibrium, at t = 5000 Ω-1 for α = 10-2, shown in blue, at t = 25 000 Ω-1 for α = 10-3, shown in red, and at t = 100 000 Ω-1 for α = 10-4, shown in yellow. The initial particle size is 10-4 R1, corresponding to approximately 0.1mm at the ice line. Solid lines are modelled values and dotted lines are analytical values, set by Eq. (53). The modelled and analytical values follow the same slope in the applicable range.

In the text

Fig. 8 Mean ice particle size ⟨ R ⟩ as a function of time in inverse Keplerian orbits Ω-1 for simulations including both radial and atmospheric ice lines. The height of the atmospheric ice line above the midplane zice is denoted as: black 0.4 H, violet 0.6 H, blue 1.0 H, green 1.4 H, red 1.8 H and yellow 3.0 H. For the runs with zice = 3.0 H, the dashed line corresponds to 1000 particles and the dotted line to 100 000 particles. Full lines denote runs with 10 000 particles. The position of the radial ice line is rice =  −0.3 H, and the strength of turbulence is α = 10-2.

In the text

Fig. 9 State of an extended simulation box with both the radial and the atmospheric ice line included, for t = 0, 100 and 1000 Ω-1, from top to bottom. The grey area represents the condensation zone, and ice and vapour is shown in blue and red, respectively. The sizes of the blue dots are proportional to the size of the ice particles and the number of particles shown is inversely scaled with size for visibility. The turbulent α-value is 10-2.

In the text

	Fig. 10 Growth in different radial subdomains of the simulation box shown in Fig. 9. Black shows the average particle size in the total box, whereas yellow, red, green and blue give the average particle sizes in the different subdomains, from left to right. The total particle size at any time is dominated by the growth near the radial ice line.
In the text

Fig. 11 Particle growth for different radial extents of the simulation box and turbulence strengths. The mean particle size ⟨ R ⟩ is shown as a function of time in inverse Keplerian orbits Ω-1 for α = 10-2 (left panel), α = 10-3 (middle panel) and α = 10-4 (right panel). The size of the box has been extended in the radial direction, with rmax colour coded as: black 0.75 H, blue 2.25 H, red 3.75 H and yellow 5.25 H. For all runs zmin =  −6.0 H, zmax = 6.0 H and rmin =  −0.75 H. The inlay figure shows the simulation domain, with the considered radial extents marked.

In the text