Ice condensation as a planet formation mechanism
Key Words.:Accretion, accretion disks - Turbulence - Methods: numerical - Planets and satellites: formation - Protoplanetary disks
We show that condensation is an efficient particle growth mechanism, leading to growth beyond decimeter-sized pebbles close to an ice line in protoplanetary discs. As coagulation of dust particles is frustrated by bouncing and fragmentation, condensation could be a complementary, or even dominant, growth mode in the early stages of planet formation. Ice particles diffuse across the ice line and sublimate, and vapour diffusing back across the ice line recondenses onto already existing particles, causing them to grow. We develop a numerical model of the dynamical behaviour of ice particles close to the water ice line, approximately 3 AU from the host star. Particles move with the turbulent gas, modelled as a random walk. They also sediment towards the midplane and drift radially towards the central star. Condensation and sublimation are calculated using a Monte Carlo approach. Our results indicate that, with a turbulent -value of 0.01, growth from millimeter to at least decimeter-sized pebbles is possible on a time scale of 1000 years. We find that particle growth is dominated by ice and vapour transport across the radial ice line, with growth due to transport across the atmospheric ice line being negligible. Ice particles mix outwards by turbulent diffusion, leading to net growth across the entire cold region. The resulting particles are large enough to be sensitive to concentration by streaming instabilities, and in pressure bumps and vortices, which can cause further growth into planetesimals. In our model, particles are considered to be homogeneous ice particles. Taking into account the more realistic composition of ice condensed onto rocky ice nuclei might affect the growth time scales, by release of refractory ice nuclei after sublimation. We also ignore sticking and fragmentation in particle collisions. These effects will be the subject of future investigations.
Planets form in protoplanetary discs of gas and dust, surrounding 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 millimeter-sizes stick together due to contact forces and kilometer-sized and larger bodies are held together by gravity. How particles grow from millimeter to kilometer-sized planetesimals is difficult to explain by collisions, as these particles tend to fragment or bounce off each other as they collide, instead of sticking (Blum & Wurm, 2008). Further, 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 meter-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 centimeter-sized to decimeter-sized pebbles, further growth into planetesimals is possible via particle concentration and gravitational collapse. The streaming instability causes particles to clump together due 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 pressure bumps (Johansen et al., 2009) excited in the turbulent flow. However, an efficient mechanism for growth from millimeter-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 assumptions of efficient condensation by homogeneous nucleation, and by ignoring inwards transport of material, the authors found a significant solid density enhancement. In contrast, Cuzzi & Zahnle (2004) considered inwards material drift 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 behavior 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, separating the inner, hotter part of the disc from the colder region further 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). Due 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 and thereby effectively move 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, as 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 decimeter-sizes on a time-scale of a few 1000 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 several scale heights outside of the ice line. The atmospheric ice line, located at where disc material is heated directly by the central star (Chiang & Goldreich, 1997), contributes less to the particle growth, generally on the order of .
The radial water ice line, or condensation front, is found at 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 (Lodders, 2003; Qi et al., 2011). These condensation fronts will be the topic of a future study.
The paper is organised as follows. In Section 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 Section 3. In Section 4 we describe the test problems used to validate the code. We present our results in Section 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 (), pebbles of centimeters to decimeters form within a few years. We discuss general assumptions and simplifications made in Section 6. Finally, in Section 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 and vertical direction, whereas the azimuthal direction is ignored, for simplicity and under the assumption of axisymmetry. The fundamental units are the gas scale height as a length unit, sound speed as a velocity unit and inverse Keplerian orbits as a time unit. The relation between these quantities is . Setting gives a system which 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 (Lodders, 2003).
Condensation is considered as a neighbour interaction. A linear scaling between number of particles and number of calculations is obtained by use of a mapping scheme, where the simulation domain is divided into a number of grid cells, with particle interactions possible between particles within the same grid cell only.
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 characterized by their dimensionless friction time (Weidenschilling, 1977), or Stokes number,
where is the mean free path of the gas and and are the material and gas density. The superscripts and 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 in units of , where so that
To find the corresponding particle size in meters, we use the material density of ice , the midplane gas density and column density from the Minimum Mass Solar Nebula (MMSN) model of Hayashi (1981),
giving the scaling of
close to the water ice line.
2.3 Condensation scheme
The rate of change in particle mass due to condensation and sublimation is
(Supulver & Lin, 2000). Here, is the thermal velocity of vapour, is the vapour density and and 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
and the saturated vapour pressure can be expressed in an equivalent way. Here, is the Boltzmann constant, is the temperature, is the vapour density and is the mass of one vapour particle. Assuming spherical particles Eq. 7 can be rewritten in terms of particle radius and using vapour densities instead of pressures,
In the limit where sublimation dominates the time evolution. The sublimation time scale in terms of particle radius is then
and in terms of mass
Correspondingly, when condensation dominates, with condensation time scale in terms of radius
and in terms of mass
In equilibrium with . By definition , where is the water mass fraction, at an ice line. The rapid increase of the temperature into the condensation region implies and hence , allowing for modelling of sublimation as an instantaneous process as long as the particle friction time, , is shorter than the characteristic time-scale of the turbulent eddies, .
Condensation can not be assumed to occur instantaneously, as never increases much beyond . The condensation time scale is proportional to the size of the particle, with the dimensionless condensation time scale being
assuming , where and are the mean molecular weights of water and molecular hydrogen. The condensation process is modelled by a Monte Carlo approach, where the probability of condensation for vapour onto an ice particle in the condensation zone is proportional to the size of the available ice surface. This is further described in Section 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 is much smaller than the number of physical particles, so that superparticle represents physical particles. In these simulations, typically for a simulation area of in the vertical direction and 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 that does not change. The mass of the physical particles represented by the superparticle does however change, and so does , the number of physical particles represented by a superparticle, as the total mass represented by the superparticle is always . Also when an ice particle undergoes a phase transition and becomes vapour the total mass of the superparticle stays constant. The superparticle approach to condensation and sublimation is sketched in Fig. 1.
2.5 Random walk of particles
In this section the modelling of particle motions is described. 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. As particles grow larger, additional effects need to be taken into account. Firstly, the random step length becomes smaller, as 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, causing 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
We here assume that particles move with velocity , and that the length a particle moves coherently is , which gives rise to by mixing length arguments. The correlation time is defined as
This means that the effect of turbulence on a small particle is to move it in a random direction with a characteristic speed , 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 , motivated by observations of accretion rates in young stars (Hartmann et al., 1998). However, the prescence of a dead zone can give a lower -value (Fleming & Stone, 2003; Oishi & Mac Low, 2009), and there are also observations suggesting in protoplanetary discs (Mulders & Dominik, 2012). Therefore we also run simulations with and .
The distance a particle moves each time step, where one time step is , is found by equating the generic expression for the diffusion coefficient of a random walk in 2D, ), with that of diffusion described by Eq. 15, giving
Practically, the random walk for particles coupled to the gas is modelled by setting the length the particle moves in the radial direction and in the vertical direction to
Here is chosen randomly for each particle and time step, giving a random direction, but a fixed length, for how a particle moves in one 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, as larger particles move shorter distances with the turbulent eddies each time step. The diffusion coefficient is modified to
following the analytical theory developed and tested by Youdin & Lithwick (2007). This gives a modified size-dependent step length of
As 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 affects the vertical direction only, it is only necessary to modify . This is done by changing the step size in the vertical direction to , where is an increasing function of particle size. The sedimentation step length is set by the terminal velocity of the particles . 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
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 arising 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 , and sedimentation-diffusion equilibrium occurs at particle scale height
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
Effectively we let even tiny particles sediment in order to follow the Gaussian distribution of the gas, rather than the uniform distribution resulting from a pure random walk. Changing in Eq. 25 reproduces the correct particle scale height in Eq. 24. With the modified dimensionless friction time in the terminal velocity expression, the step size in the vertical direction can be written as
However, this expression leads to a numerical instability in which large particles overshoot the midplane as they sediment, since for and we get
This gives particle oscillations that are amplified in time. To eliminate this problem, Eq. 27 is modified to
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 due to the radial drift towards the central star. This is also particle size dependent, and is modified to
where is the Keplerian velocity and is the velocity difference between the gas and the particles, following Weidenschilling (1977). This gives a radial drift velocity that peaks for particles. We use , a reasonable value at in the MMSN (Cuzzi et al., 1993; Chiang & Youdin, 2010).
The friction time of a particle is inversely proportional to the gas density, , 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
Simulations have shown that the heterogeneous vertical gas distribution leads to a turbulent -value increasing 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
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 particle Stokes drag is the relevant regime. The transition size given by Eqs. 1 and 2 can be rewritten in dimensionless friction time as
Converting dimensionless friction time in the Epstein drag regime to the Stokes drag regime is done using
where the scale height is given by the MMSN model (Hayashi, 1981) as
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 , 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, any vapour particle can only condense onto one single ice particle, so it can not be shared amongst several particles. One ice particle can on the other hand experience several growth events during one time step, if 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 . 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 is
where is the condensation time scale for particle , given by Eq. 14. The expression for interaction probability is chosen to always give , for any , also for . The algorithm proceeds by calculating the probability that no interaction occurs. For ice particles, the total probability that no interaction occurs is
To decide whether condensation occurs, a random number is generated. If is smaller than nothing happens, whereas if is larger than vapour 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. 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
The relative probability for each of the ice particles in the grid cell are placed in a sequence,
and a new random number is generated, . Which interval in the sequence falls in decides which ice particle the vapour particle condenses onto, such that falling in the interval implies condensation onto ice particle .
The condensation time scale in Eq. 13 is proportional to the radius of the particle, so that smaller particles have a shorter, and larger particles 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 . 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 timesteps 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 due to turbulence, stirring them up in a turbulent diffusion, and due to 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. Fig. 2 shows the particle scale height relative to the gas scale height as a function of particle size . 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,
where is the number of particles and 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,
Fig. 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 , where at . From the figure one can see that small particles are fully coupled to the gas, as for particles of , corresponding to mm-sized ice particles near the water ice line. Further, a slight change of the slope can be seen at . 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 as a condensation event occurs. As particle radius is connected to mass as , for any number of condensation events the radius increases as
The final radius is plotted as a function of as black crosses in Fig. 3. In this test 2000 ice superparticles were distributed evenly in 20 different initial size bins over the range , 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 . 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.
The modelled particle growth can be compared to the analytically expected . For a physical ice particle the total mass growth due to condensation during the time is
where is the material density of ice, is the initial and 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 , where 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,
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. As 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.
The right handside of Eq. 46 can be rewritten using
Here 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 handside Eq. 46 can be rewritten as
where and 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
with being the number of vapour superparticles added to the simulation. Eq. 46 therefore is equivalent to the very simple expression
which the code can be checked against, as 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 , equal for all particles, can be found. The final radius can be rewritten as , where 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. Eq. 46 is thereby rewritten as
The initial number of vapour and ice superparticles and are known, as they are input parameters to the simulation. The radius represented by a superparticle in the end of the simulation, , and when it is created, , are both given as output from the code. The radius growth for each particle can therefore be found. For and , a value of is found. In Fig. 3 the blue full line represents the expected growth of the 2000 initial ice superparticles, , as a function of initial radius . 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. averaging over many particles, the result can be considered to be correct (Zsom & Dullemond, 2008).
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 run 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 Section 6.4.
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 (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 . The pressure ice line in a vertically isothermal disc is shown in Fig. 4.
The height of the atmospheric ice line above the midplane is varied from to . 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 which 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.
Fig. 5 shows the time evolution of a simulation with the atmospheric ice line at . The distribution of red vapour and blue ice particles is shown together with the size distribution of ice particles, for , and . 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 decimeter-sized ice pebbles residing in a midplane layer.
Varying the atmospheric ice line position
Simulations were run with the atmospheric ice line placed at different heights above the midplane. Decreasing corresponds in a physical protoplanetary disc to placing the simulation box radially closer in towards the radial ice line, where a hypothetical 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 the particle growth for different is shown. The mean ice particle size in units of is shown as a function of time in units of . The curves show, from top to bottom, (black), (violet), (blue), (green), (red) and (yellow). 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. Fig. 2 shows how the particle scale height decreases with increasing particle radius . Small particles have a scale height comparable to the gas scale height , whereas for larger particles . Placing the ice line at a height larger than 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 is both slow, and to modest particle sizes. For ice line positions of the particle growth is on the order of or less in . For the larger extent of the disc, where , 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.
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 is shown as a function of the height of the ice line above the midplane , starting from an initial particle size of . The system is shown when it has approximately reached an equilibrium, at for , shown in blue, at for , shown in red, and at for , 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 as 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
Setting gives an analytical estimate of the maximum particle size achievable by condensation, shown in Fig. 7, as
We chose , but even in this case the analytical expression underestimates the resulting slightly. Nevertheless our simulations show that particles grow approximately to a size where their scale-height is 1/3 of the ice line height.
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 , 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 in the vertical direction. In the radial direction the box size is varied from to , 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 and , as 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 conditions such that ice sublimates prevail at all , whereas at 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, causing a too large particle growth.
Varying the atmospheric ice line position
Including the radial ice line greatly increases the growth efficiency. As shown in Sec. 5.1, models with only the atmospheric ice line gives negligible growth for ice line positions of . However, taking also the radial ice line into account leads to particle growth beyond , corresponding to decimeter-sized pebbles at the ice line, within .
Fig. 8 shows how the mean particle size as a function of time varies with atmospheric ice line position when also the radial ice line is included. The colours show different heights of the atmospheric ice line above the midplane, with black denoting , violet , blue , green , red and yellow . Growth is somewhat faster for close to the midplane and slower for further 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 position above , a similar growth as for is found.
Growth stops at for all . This is due to radial drift towards the star, which eventually removes all ice particles from the simulation. From Eq. 30 it can be seen that the radial drift inwards peaks for particles of . As 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 for a smaller () and a larger () number of particles in Fig. 8, in addition to particles. For particle sizes 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 decimeters, but may be unphysical in the regime of meter-sized boulders and beyond.
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. Fig. 9 shows the time evolution of a run where the radial extent of the simulation box has been increased, so that . The most significant growth takes place close to the radial ice line, but due 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, . The growth timescale 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 () and right () panels of Fig. 11. These low-turbulence simulations show that, although slowly, particles can grow via condensation also in dead zones.
Our results show that ice condensation is an efficient way to form centimeter-sized and decimeter-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 which we discuss here.
6.1 Other ice lines
This model has been developed primarily to investigate the water ice line at . Similar ice lines, or condensation fronts, exist also for other chemical species, and as this model is scalable to any 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 further out in the protoplanetary disc, and the numerous silicate condensation fronts further in towards the central star. A global disc model, including 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, as the ice line probably moves both in a systematic way over long time scales and due to shorter heating events (Armitage, 2011). However, as the time scale for growth by condensation found in this work is on the order of , which is very short compared to the life time of the disc , 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).
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 ice 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 the 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 nanometer-sized (Draine, 2003). The largest dust particles of relevance are the ones of the largest size that are strongly enough coupled to the gas to be present at the relevant vertical distance from the midplane. This is around micrometer-scale (see Fig. 3). The Kelvin curvature effect states that the equilibrium vapour density for a curved surface is higher than for a flat surface, meaning that vapour more easily condenses onto a flat surface than a curved one, an effect which is most important for the smallest particles, around nanometer-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. As 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 . 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 ), 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.
Further, 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 millimeters, 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 () 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 () (Wettlaufer, 2010).
Collision experiment for higher velocities have been used to derive criticial relative velocities above which fragmentation can occur, . Higa et al. (1996) found a critical velocity of , but other groups have found significantly larger values (Arakawa, 1999; Arakawa et al., 2002). Computer simulations suggest critical velocities of up to (Benz & Asphaug, 1999).
Collision velocities for millimeter-sized particles are set by the turbulent velocities and are expected to be on the order of a few , whereas collisions involving larger particles approach the radial drift velocity, (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 Section 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.
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 centimeters, 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 millimeter-sized dust grains to decimeter-sized pebbles is . Significant growth is obtained for a range of turbulent -values from to , 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. Similarly, the outer disc regions have 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 vapour can still exist.
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 oftemperature 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 decimeter-sized pebbles in protoplanetary discs, thereby providing the missing link to further growth into planetesimals and planets.
Acknowledgements.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 improved the quality of the paper. This research was partially funded by the European Research Council under ERC Starting Grant agreement 278675-PEBBLE2PLANET.
- email: firstname.lastname@example.org
- email: email@example.com
- Arakawa, M. 1999, Icarus, 142, 34
- Arakawa, M., Leliwa-Kopystynski, J., & Maeno, N. 2002, Icarus, 158, 516
- Armitage, P. J. 2011, ARA&A, 49, 195
- Barge, P., & Sommeria, J. 1995, A&A, 295, L1
- Blum, J., & Wurm, G. 2008, ARA&A, 46, 21
- Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859
- Benz, W., & Asphaug, E. 1999, Icarus, 142, 5
- Bridges, F. G., Supulver, K. D., Lin, D. N. C., Knight, R., & Zafra, M. 1996, Icarus, 123, 422
- Carballido, A., Fromang, S., & Papaloizou, J. 2006, MNRAS, 373, 1633
- Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368
- Chiang, E., & Youdin, A. N. 2010, Annual Review of Earth and Planetary Sciences, 38, 493
- Cuzzi, J. N., Dobrovolskis, A. R., & Champney, J. M. 1993, Icarus, 106, 102
- Cuzzi, J. N., & Zahnle, K. J. 2004, ApJ, 614, 490
- Draine, B. T. 2003, ARA&A, 41, 241
- Dullemond, C. P., & Dominik, C. 2004, A&A, 417, 159
- Fleming, T., & Stone, J. M. 2003, ApJ, 585, 908
- Fromang, S., & Nelson, R. P. 2009, A&A, 496, 597
- Goumans, T. P. M., Catlow, C. R. A., Brown, W. A., Kästner, J., & Sherwood, P. 2009, Physical Chemistry Chemical Physics (Incorporating Faraday Transactions), 11, 5431
- Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385
- Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35
- Higa, M., Arakawa, M., & Maeno, N. 1996, Planet. Space Sci., 44, 917
- Hogerheijde, M. R., Bergin, E. A., Brinch, C., et al. 2011, Science, 334, 338
- Hughes, A. L. H., & Armitage, P. J. 2010, ApJ, 719, 1633
- Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022
- Johansen, A., Brauer, F., Dullemond, C., Klahr, H., & Henning, T. 2008, A&A, 486, 597
- Johansen, A., Youdin, A., & Klahr, H. 2009, ApJ, 697, 1269
- Klahr, H. H., & Bodenheimer, P. 2003, ApJ, 582, 869
- Krivov, A. V., Müller, S., Löhne, T., & Mutschke, H. 2008, ApJ, 687, 608
- Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32
- Lecar, M., Podolak, M., Sasselov, D., & Chiang, E. 2006, ApJ, 640, 1115
- Lodders, K. 2003, ApJ, 591, 1220
- Meijerink, R., Pontoppidan, K. M., Blake, G. A., Poelman, D. R., & Dullemond, C. P. 2009, ApJ, 704, 1471
- Morbidelli, A., & Nesvorny, D. 2012, A&A, 546, A18
- Mulders, G. D., & Dominik, C. 2012, A&A, 539, A9
- Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375
- Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43
- Oishi, J. S., & Mac Low, M.-M. 2009, ApJ, 704, 1239
- Papoular, R. 2005, MNRAS, 362, 489
- Qi, C., D’Alessio, P., Öberg, K. I., et al. 2011, ApJ, 740, 84
- Safronov, V. S. 1969, 1969.,
- Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
- Sirono, S.-i. 2011, ApJ, 735, 131
- Stevenson, D. J., & Lunine, J. I. 1988, Icarus, 75, 146
- Supulver, K. D., & Lin, D. N. C. 2000, Icarus, 146, 525
- Vitense, C., Krivov, A. V., Kobayashi, H., Löhne, T. 2012, A&A, 540, A30
- Weidenschilling, S. J. 1977, MNRAS, 180, 57
- Wettlaufer, J. S. 2010, ApJ, 719, 540
- Wilner, D. J., D’Alessio, P., Calvet, N., Claussen, M. J., & Hartmann, L. 2005, ApJ, 626, L109
- Windmark, F., Birnstiel, T., Güttler, C., et al. 2012, A&A, 540, A73
- Wurm, G., Paraskov, G., & Krauss, O. 2005, Icarus, 178, 253
- Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459
- Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588
- Zsom, A., & Dullemond, C. P. 2008, A&A, 489, 931
- Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57