We have previously narrated a formation scenario for the TRAPPIST-1 system, based on pebble accretion. We have now followed this up with a numerical model by combining three types of codes: (i) a Lagrangian code for dust evolution; (ii) an N-body code for planetesimal (collisional) dynamics; and (iii) an Analytical code for pebble accretion. In this way we model grain growth, pebble radial drift, planetesimal formation at the H2O iceline, planetesimal runaway growth, planet migration, pebble accretion, and follow the water fraction of the planets.

This model gives us in short a synthetic planetary system, for the pebble-driven formation scenario. In the second figure we compare the masses and radii of the simulated planet system (blue points) with the observations of TRAPPIST-1 (red error bars). Lines denote mass-radius relationships for a composition that amounts to H2O fractions of 0, 10, and 20% by mass. The blue arrows indicate the effects of the evaporation of water.

Our results reproduce by-and-large the physical properties of the TRAPPIST-1 planets (similar-sized planets and composition) and newly quantify the spread in these quantities. Note that in the figure, we have not accounted for pebble isolation, which could have limited growth of the larger planets. Post-formation evaporation would likewise change the planet properties.

The model setup.

Cloud formation is complex, but schematically the picture is clear. High in the atmosphere the conditions have become such that vapor species can condensate on tiny seed particles (aerosols or nuclei). These cloud particles growth, through further condensation or coagulation, whereafter they are transported down by gravity. In big exoplanets, without a solid surface, the cloud species then re-evaporate in the hotter lower regions and the vapor can be transported upwards to fuel new clouds, completing the cycle

The figure presents examples of cloud profiles, where we vary the nucleation rate (a parameter of the model) and switch coagulation on or off. Models accounting for coagulation produce larger grains. It can be seen that a higher nucleation rate results in small particles, which makes sense because there are more particles that share the vapor. However, accounting for coagulation ameliorates the differences somewhat.

grain size (x-axis) and intensity (color) of KCl condensates in GJ 1214b for several choices of the nucleation rate.

(Refractory) carbon is plentiful in the ISM and is the key building block of any organic compounds and life. Yet, there is almost none of it left on Earth. Where has all the carbon gone?

We therefore arrive at the significant conclusion that the low carbon abundance in the inner solar system provides a major constraints on its formation. Perhaps powerful FU Ori-like events flashed the carbon and/or the inner disk was dynamically isolated from the outer disk by an early formation of Jupiter.

Sketch of how carbon (black) gets destroyed

What is the likelihood of a drifting pebble to be accreted by a planet? This number — ε — depends on the ratio of two quantities: the pebble accretion rate by the planet and the drift speed of the pebble. These quantities in turn depend on many other parameters (planet mass, planet eccentricity and inclination, stellar mass, turbulence, disk headwind, etc.). In two works we systematically calculate &epsilon. In the first paper, we consider a global frame (star at the center), accurately accounting for the curvature of the disk. In Paper II, we follow this up with a stochastic equation of motion for the pebbles, with which we can mimic the forcing by a turbulent flow.

We finally formulate analytical expressions for ε — general enough to be used in N-body simulations. The expressions fit the result quite well as shown by the curves in the figure. Altogether, this study is in line with our previous works, namely that pebble accretion is generally not (super-)efficient and that a massive pebble reservoir (disk) is needed to grow planets.

pebble accretion efficiencies in 3D for various turbulence levels (colors) and particle stopping times (τ)

Because of settling, planets embedded in disks can accrete pebbles at a large cross section, under the condition that pebbles can actually reach the planet. Large planets may open a (shallow) gap in the disk, resulting in pebbles piling up at pressure maxima away from the planet. Small planets, likewise, affect the gas disk in their vicinity, resulting in streamlines that curve away from the planet, such as the well-known horseshoes. Because of their strong coupling to the gas, small particles are particularly affected.

accretion of pebbles in realistic gas flows

Sub-Neptune planets are special. While most of their mass resides in the "rocky" core, their radius is determined primarily by the state of the gaseous envelope. The standard assumption is that the radius is set by a combination of the composition and mass of the envelope. But in this Letter we investigated another possibility: the (remnant) heat of the core that is slow to "escape". The heat of formation (gravitational binding energy) represents a huge amount of energy, which will indirectly heat the envelope, causing its contraction to stall. Of key importance is here the time it takes for the heat to "leak out". If this occurs quickly, the effect on the radius is large, but the effect soon withers away (blue curves). Similarly, when the core releases its heat on a very long timescale such that the leakage is presently insignificant, there is no effect. Only when the release timescale is similar to the age of the planets (~Gyr) is there a noticeable effect on the radius of the planet (see inset).

So, what is the true timescale on which the core releases its heat? That is a very hard question, which requires us to know the thermodynamic state of the core. Unless the assumption is made that these cores somehow behave identical to Earth, it is anyone's guess...

Evolution of the radius with time for several internal heating models

Pebbles don't just get accreted when they enter a pre-planetary atmosphere — they also get heated! Since the pre-planetary atmosphere (the concentration of gas around a gravitating body in a gas-rich disks) can become very hot, the small pebble may even not reach the core of the planet! Instead, it would end up in vaporized form in the atmosphere. We conducted a 1D calculation to find the mass where this occurs of a pebble of silicate composition and found complete vaporization for a planet a fraction of Earth's. Clearly, growth of a planets accreting pebbles proceeds very differently than accretion of planetesimals!

See this astrobite article by Michael Hammer for a great summary!

The fraction of the pebble that gets ablated and absorbed by the atmosphere (blue). Note that the y-scale is log and the x-scale is linear.

FU Ori outburst are events where the stars becomes very luminous and the disk very hot, presumably for a short amount of time. This makes these systems ideal to study phenomena related to the H2O iceline. One such star is V883 and here we interpret recent ALMA observations, constraining key properties of the pebbles interior and exterior to the H2O iceline.

How do super-Earth and mini-Neptunes "survive" in gas rich disks? The problem is that an early formation of a super-Earth is preferred (because they often harbor a massive hydrogen/helium-rich atmosphere), yet not too much gas can be accreted, as this would have converted them into hot-Jupiters. Since super-Earth planets are a factor 100 more frequent than hot-Jupiters there should be a robust process to prevent rapid gas accretion. We have previously suggested that the exchange of material between the super-Earth' atmosphere with the gas of the circumstellar disk is the key to solve the problem. In this work, led by (then) undergraduate student Nicolas Cimerman, we specifically tested this idea by running hydrodynamical simulations of the atmosphere of super-Earth planets, located at a distance of 0.1 au from their star. The animation on the right illustrates the effect of recycling by plotting the entropy of the atmosphere gas. This entropy decreases with time, because of radiative cooling of the atmosphere. However, the decrease is much stronger in the 1D bound (and isolated) atmosphere models (solid curve) compared to our 3D hydrodynamical models (dashed curve). This clearly shows that atmospheric recycling interferes with radiating cooling and Kelvin-Helmholtz contraction of these envelopes, and that super-Earths can survive "exposure" of a gas-rich disk.

evolution of the entropy -- for a 1D bound, isolated atmosphere and a 3D "open" atmosphere

stages in the formation of the TRAPPIST-1 system

I wrote a review chapter about the pebble accretion mechanism, which will become part of a book on "Formation, Evolution, and Dynamics of Young Solar Systems" (Springer Verlag; edited by Oliver Gressel and Martin Pessah). It covers the physics of the process, assessing the cons- and pros, and reviews some recent applications. Enjoy!

Collisional cross section compared to the geometrical as function of the accretor (Y-axis) and the pebble (x-axis). White indicates geometric (slow) growth

How do planetesimals form from pebbles? According to the so-called Streaming Instability model this happens once the pebble-to-gas density is around unity. The problem with pebbles is, however, their mobility: they tend to drift towards the star in an inside-out-fashion, meaning that the inner disk is emptied first. Unfortunately, this behavior tends to decrease the pebble spatial density, instead of the desired increase. Something special is needed.

Compared to previous works, we have investigated the problem in a dynamical framework, where both solids and gas accrete onto the star. Because pebbles do not stop it is harder to get the required pileups, but we nonetheless show that this is feasible with our more realistic approach.

Formation of a high solids-to-gas density ratio at the iceline

The Kepler spacecraft has discovered an enormous number of exoplanets of the super-Earth type: bigger than Earth, up to Neptune in sizes, orbiting very close to their host stars. Mostly, these planets are not alone: they are often part of a compact, multi-planetary system. Intriguingly, some of these planets are seen in resonance (a resonance is an integer commensurability in orbital periods, like 2:1 or 3:2); but most of them are not.

Magnetospheric rebound.

Broadly, there are two distinct philosophies for solving the dust coagulation problem in disks: (i) make the model as precise as you can; (ii) include the key physics, but optimize it for speed. Because of its toy model nature, the latter approach is often preferred as it better lends itself to explore a vast parameter space. Here, we have developed such a panoptic model for planet formation, where we model the first stages of dust coagulation in the disk. The figure shows an example of trajectories of coagulating dust particles. As is well known, dust first grows (vertical lines) and then drifts (horizontal lines). For each of these trajectories we also follow the particle size, their porosity (not shown) and the density (background shading). Lines terminate when particles have grown large enough to form planetesimals. Here we see that the largest radius where this occurs is ~5 AU, which we predict is the radius where the first planets can form (because they will sweepup the remainder pebbles that drift in).

'Lifelines' of particles that grow and drift

It has been shown that pebble accretion can be very efficient for optimal conditions of the pebble (aerodynamical) and protoplanet properties. But when does pebble accretion start? In particular, when gravity is weak small particles are carried with the flow streamlines around planetesimals, hence avoiding accretion. In this work, we have conducted a systematic study on how efficient gas drag-mediated interactions are and determined the collisional cross section, see figure. Values >1 imply gravitational focusing (fast growth), values <1 aerodynamic deflection (grow stalls). The two dotted lines denote the transition from the geometric (collisional=geometric cross section) regime to aerodynamic deflection and Safronov regimes, respectively. Also, the curved dashed line separates ballistic from settling encounters. Pebble accretion start to the right of this line. For typical disk, we find that planetesimals of around 100 km grow the slowest, indicating a bottleneck for sweepup growth at this size.

The collisional efficiency factor: >1 implies gravitational focusing, <1 aerodyanmic deflection.

Small planets dominate the universe and even large planets were once small. As planets form in gas-rich disks a key question is how the planet's gravity shapes the gas around it, and in particular what the properties of the planet's primordial atmosphere are. In this study hydrodynamical simulations were conducted to solve for the flow pattern of the gas. Due to the gravity, the gas compresses close to the planet. This process is illustrated in the video to the right, where the white color indicates the density of the background gas (that of the circumstellar disk) and the red much higher densities. Naturally, the gas compresses near the surface of the planet, because the gravitational pull of the planet is strongest.

Solid curves — streamlines — give the trajectories of the flow. Streamlines near the planet center are closed: the gas is bound to the planet. Of particular interests is the blue — critical — streamline, which gives the boundary between the region where gas circles around the planet and gas that orbits the star or moves on so called horseshoe orbits (which make a U-turn). They thus give the size of the atmosphere.

In this work a correlation was found between the size of the critical streamline, that is, of the atmosphere of the planet, and the amount of gas mass that the atmosphere contains. For more massive planets, or for planet atmospheres that are highly compressible (meaning: cold), rotational motions are stronger. This result is a restatement of Kelvin's circulation theorem. A consequence of the increased rotation is that it sets an upper limit on the total amount of gas that these atmospheres can have. This upper limit is reached when the gas starts to rotate at Keplerian velocities, signifying the transition from pressure support to rotational support.

These findings are, however, based on two dimensional (2D) hydrodynamical simulations, meaning that the vertical dimension is omitted. In reality the streamline geometry for these embedded planets is 3D and very complex as shown in the figure below. Here, the projections of the 3D streamlines on two planes are shown. One such streamline (solid purple) originates from considerable height, but then spirals in to end up orbiting the planet. However, as gas does not pile up, gas is also expelled from the region near the planet and returns to the disk (black solid).

Thus, in these 3D simulations the atmosphere is `open' in the sense that it is continuously being replenished by gas from the protoplanetary disk. For these planets the replenishment time is fast, especially for planets close to the star. This is unsurprising because the background gas is denser and the planet rotates faster closer to the star; gas proceeds through the planet's atmosphere at a faster pace. The fast replenishment will be a key determinant in the evolution of these atmospheres. One implication is that such atmospheres have little time to cool, contract, and to become more massive. Thus, although the 2D and 3D simulations differ qualitatively, both indicate that there are limits to the atmosphere mass that these embedded planets can support. This is an encouraging finding in the light of the ubiquity of super-Earths and (mini-)Neptune size planets: planets that are quite massive but nonetheless ended up with only a tiny atmospheres compared to their rocky cores.

Emergence of a steady 2D flow and a bound atmosphere around a planet embedded in the primordial, gas-rich disk

Example of a calculation with the new structure equation for grain growth. The left panel shows the normalized temperature and density while the right panel plots the grain properties: their size, abundance, and opacity. Results of two simulations are presented: (i) grains only enter at the top of the atmosphere (solid lines); (ii) grains are additionally produced by ablating planetesimals within the atmosphere (dashed lines). In the latter grains coagulate quickly, but the bigger particles suppress the grain opacity. As a result, the temperature and density profiles are very similar.

Comparison of turbulent motions (black dashed), strength curves (red solid and dashed) and the escape velocity (blue)

Pantha rhei ("Everything flows" — Heraclitus, ~500 BC); but how do gravitating bodies like planets affect the flow pattern of the gas that attempts to stream past? When protoplanetary embryos made from accumulation of solid particles exceed a certain threshold mass (corresponding to ~1000 km in size), they can start to bind the gas from the nebula disk. The atmospheres of these young, low mass planets are (presumably) hot; generally, the gas does not collapse onto the planet and the atmosphere is in pressure-equilibrium with the disk. In fact, the boundary between planet and nebular disks has to be determined from the velocity of the flow.

We have calculated the steady-state flow pattern which emerges in the vicinity of the planet (see figure). Left and right one notes the background shearing flow, which is almost unperturbed. Towards the top and the bottom, one sees, very clearly, the horseshoe region where streamlines make a "U-turn". At the very center the flow curls around the planet in the prograde (counterclockwise) direction. This is the atmosphere of the protoplanet.

Implications concern the migration behavior of protoplanets (for which the width of the horseshoe region is a key parameter), the thermodynamical structure of the protoplanet atmosphere, circumplanetary disks formation, and the accretion behavior of small particles (for which the gas drag is important).

A typical flow pattern past a low mass planet, where it is in steady state. Arrows indicate the velocity and gray contours are streamlines of the flow. Red/Pink circles are isodensity contours.

If you ever shot a gun, you probably noticed that the gun recoiled back on you. This is due to the conservation of momentum. In disks, a migrating planet likewise reflects that its (angular) momentum is changing. The gravitational interaction with the gaseous disks is a well-known effect (Type-I migration).

In a disk solid bodies (planetesimals) act as the bullets. The gravitational interaction with the much more massive planet will slingshot them to different orbits (scatterings). The planet feels the recoil, which causes it to migrate. Overall, this planetesimal-driven migration is analogous to the (more well-known) Type-I migration; but both can be understood as a consequence of dynamical friction.

A complication is that, to first order, interactions with the interior disk cancel interactions with the exterior disk. Therefore, higher order effects must be included. This means that the migration depends on the gradient in the surface density and eccentricity. We investigate the effects of an eccentricity gradient and find a strong dependence. In addition, we find a regime where the migration is self-sustained.

Mapping the three migration regimes

It is already difficult to make one planet — let alone two. Therefore, we have investigated the idea for triggered planet formation. Applied to the solar system this means: form Jupiter first, then Saturn.

To make life a bit easier, we have assumed that Jupiter did already form (without specifying how) and carved a wide gap in the primordial gas-rich protoplanetary disk. This gap causes a pressure maxima, whose location could coincide with Saturn's for plausible parameters. This is important because debris ('small stuff') will pile up at these pressure maxima. The debris originates from collisions among planetesimals from the outer disks.

The figure on the right tells the story: without a pressure maxima (dashed lines) Saturn will not grow big — it doesn't even get to the Earth! With a pressure maxima created by Jupiter (solid line), it easily jumps over the Earth in terms of mass and become a gas giant next to Jupiter. Note that this mechanism works better when the outer disk contains smalle bodies, because these are weaker and collide more frequent.

Evolution of protoplanet mass for several models

A flowchart of the toy model

Coagulation will affect the dust size distribution in dense molecular clouds. This, in turn, affects the dust opacity — a critical quantity required for any interpretation of observational data sets. Following previous work (below) where we computed the dust size distribution as function of time, we now present the corresponding mass-weighted opacities for infra-red wavelengths. The figure shows that the opacities change on timescales of ~Myr (or less if the cloud's density is higher than the assumed n=10⁵ cm^-3). At visible and near-IR wavelengths the opacitiy decreases, but at longer wavelengths it will increase with time. We have quantified this evolutionary trend in terms of the strength in the 9.7μm silicate feature vs near-IR color excess and in terms of the sub-mm slope β.

Opacity changes with time/coagulation state

Examples of orbits for particles experiencing various amounts of gas drag and gravitational forces

It is well known that a collisional cascade produces a powerlaw size distribution. The most notably example being the dust size distribution in the interstellar medium, the so called MRN distribution. This situation represents a steady-state: collisions among particles of a certain size deplete their number at the same rate as the replenishment rate by collisions among larger particles. A collisional cascade is applicable to situations where collisions results in fragmentation. However, in a protoplanetary disks the gas damps the motions of the smaller particles. As a result, the small particles coagulate whereas only collisions with higher-mass particles are energetic enough to produce fragments. In this situation a steady-state emerges but the size distribution is now determined by the characteristics of both the fragmentation and the coagulation.

We have considered a coagulation-fragmentation steady state for the dust in protoplanetary disks. Two key parameters are: i) the velocity field (Brownian motion or turbulence); and ii) the properties of the fragmentation event (the size distribution of particles within a collision). In case of a power law dependence on mass of these quantities we have derived the ensuing exponent for the size distribution.

Size distribution for different fragmentation states

When bodies are large enough to attract each other gravitationally the effective cross section for collisions becomes larger than the geometrical cross section. This enhancement factor is called the gravitational focusing factor (GFF) and is approximately given by the square of the ratio of the escape velocity of the protoplanet and the relative velocity at wich the bodies approach. The nature of this phenomenon is such that the largest bodies grow faster than other bodies, a situation referred as runaway growth (RG). RG indicates a positive feedback effect: due to the growth, the gravitational focusing increases. RG or, generally, a large GFF is required in order to grow protoplanets in a sufficiently short time span.

However, the positive feedback of RG is counteracted by viscous stirring, the (long-range) deflection of bodies' trajectories that causes the mean relative motion to increase. The growth of the protoplanet then slows down, switching to the much slower oligarch growth stage.

In this study we have investigated the conditions at which runaway growth turns into oligarchic growth. In the figure the GFF is plotted as function of the (evolving) radius of the protoplanet, which indicates time. First GFF increase (the RG-stage) but decrease after the transition size R₁=R_tr (the oligarchic growth stage). We provide a criterion for the start of the oligarchic growth phase (i.e. for R_tr) in terms of environmental conditions (radius and surface density of planetsimals, semi-major axis, etc.)

Evolution of the gravitational focusing during runaway growth

When bodies reach km-sizes (planetesimals), their collisional and dynamical evolution becomes dominated by gravity. Ideally, one would model the collisional evolution by an N-body methods; however the shear number of planetesimals limit these attempts in practise. However, to model the (runaway and oligarchic) growth correctly, its discrete nature must be taken into account. For this reason we have extended our Monte Carlo /superparticle code to treat dynamical interactions.

The figure shows the size distribution at three times. Protoplanets can be seen to separate from the population of (leftover) planetesimals. The color shows the amount of dynamical excitation with respect to the largest body in the simulation: blue colors indicate a dynamical cold system, whereas red colors indicate a dynamical hot system. Click the image for an mpeg movie!

Planetesimal accretion—click on image for animation

Very small, μm-size, dust particles easily stick, setting the first steps in a coagulation process that will eventually form planets. The sticking assuming becomes less obvious for larger particles, however, with laboratory experiments indicating that equally-sized dust particles often bounce off. But in other experiments, particles can still stick, especially if the size ratio of the particles involved is large. Yet in other experiments, fragmentation is observed at larger impact velocities. The porosity of the aggregates is also important to determine the outcome of a collision.

To investigate the implications of these diverse results on the coagulation process, we have calculated the collisional evolution with a Monte Carlo code. Each particle is characterized by two properties —its mass and porosity— and represents a certain share of the total dust's mass budget. The Monte Carlo code calculates the probability of a collision among any combination of particles with the outcome of this collisions being given by (or interpolated from) the laboratory results.

Initially we find (click image to start movie) that the growth is fractal: the porosity (=enlargement factor) of the bodies increases. However, at some stage the fractal growth stops, as collisions become more energetic. Thereafter, there is a sudden (almost runaway) growth stage that is triggered by a wide size distribution. However, this growth is subsequently negated by fragmentary or mass-transfer collisions. In the end, bouncing predominates and little evolution is present. We find that the bouncing result is a near-universal outcome, i.e., it seems very difficult to circumvent.

In a follow-up study we have extended the model to include a vertical dimension to investigate the sedimentation/diffusion behavior of the dust. One of the underlying questions is whether or not the disk atmosphere can be kept dusty (meaning that small particles are present) on long timescales. Dusty disk atmospheres are observationally favored to explain the presence of, among others, the 10μm silicate feature. We found that, provided turbulence is strong, particle fragmentation can indeed replenish small particles. However, the drawback is that particles at the midplane —where planets presumably form— also do not grow large!

Combined evolution of particle mass and porosity. Click image for animation

Cores in molecular clouds are dense enough for grains to collide — a process which has implications for the interpretation of their observations. Not all the collisions are necessarily sticking, however; the relative velocities among dust particles is relatively large and the material properties will likewise influence the growth process. In this study we have investigated the process of dust coagulation and fragmentation by combining an state of the art molecular dynamics model for the outcome of collisions among individual dust aggregates and a Monte Carlo collisional evolution code.

We consider two type of coagulation modes: i) among slicate-like grains; and ii) among ice-coated grains. It is found that the former state leads (quickly) to steady-state dust distributions, where destructive collisions among high-mass aggregates replenish the small grains. However, if the grains are ice-coated their growth can become significant (see figure). The material properties of ices are very conducive for grain growth, perhaps up to a factor 10³ in size.

Whether or not this strong growth materializes also depends on the lifetime of the molecular cloud (or core). If the condensation is simply a transient phenomenon, e.g., a fluctuation in a turbulent environment, it lifetime — essentially the free-fall time — is probably only ~10⁵ yr. In that case, the imprints of grain growth are probably very minor. However, magnetic fields may retard the collapse, perhaps by a factor of 10 or so, promoting growth.

Observationally, grain growth manifests itself by the reduction of the opacity. We have calculated the effects of grain growth on the (geometrical) opacity and found that the initial decline can be understood in terms of the initial collision timescale between dust grains. Thus, a higher density compensates a shorter lifetime to give a similar observational signature. At large times, fragmentation and the replenishment of smaller grains stabilizes the opacity. More sophisticated opacity calculations are upcoming.

Evolution of the grain size distribution

Monte Carlo methods simulating the coagulation process suffer from one fundamental shortcoming: their limited dynamic range. Coagulation removes particles from the distribution and the number of collisions that can be followed is necessarily less than the initial number of particles N. The latter is of course very large for astrophysical purposes, but is in practice small due to computational reasons.

In this study we have implemented a new algorithm, introducing the superparticle concept to Monte Carlo simulations; The number of superparticles N_g is then limited, but N can be virtually infinite. Collisions are then between particle groups; rather than between individual particles. The algorithm leaves the user the freedom to choose the relation between superparticles and physical bodies.

In our favorite implementation, we assign the superparticles equally over logarithmic mass space, such that the high-mass, exponentially declining tail of the distribution is well resolved too. This means that the algorithm is ideally suited to study runaway growth systems. An aplication is presented in the figure, where the timescale for runaway growth (a.k.a. gelation) is plotted vs. initial partice number N for two runaway kernels. The larger the box size the sooner the gelation proceeds (although the particle density is constant in all cases).

Runaway (gelation) growth time as function of system size