Recently seven Earthsized planets were discovered around the Mdwarf star TRAPPIST1. Thanks to transittiming variations the masses and therefore the bulk densities of the planets have been constrained suggesting that all TRAPPIST1 planets are consistent with water mass fractions on the order of 10. These water fractions as well as the similar planet masses within the system constitute strong constraints on the origins of the TRAPPIST1 system. In a previous work we outlined a pebbledriven formation scenario. In this paper we investigate this formation scenario in more detail. We used a Lagrangian smoothparticle method to model the growth and drift of pebbles and the conversion of pebbles to planetesimals through the streaming instability. We used the Nbody code MERCURY to follow the composition of planetesimals as they grow into protoplanets by merging and accreting pebbles. This code is adapted to account for pebble accretion typeI migration and gas drag. In this way we modelled the entire planet formation process pertaining to planet masses and compositions not dynamical configuration. We find that planetesimals form in a single early phase of streaming instability. The initially narrow annulus of planetesimals outside the snowline quickly broadens due to scattering. Our simulation results confirm that this formation pathway indeed leads to similarlysized planets and is highly efficient in turning pebbles into planets. Our results suggest that the innermost planets in the TRAPPIST1 system grew mostly by planetesimal accretion at an early time whereas the outermost planets were initially scattered outwards and grew mostly by pebble accretion. The water content of planets resulting from our simulations is on the order of 10 and our results predict a Vshaped trend in the planet water fraction with orbital distance from relatively high innermost planets to relatively low intermediate planets to relatively high outermost planets.
['accretion', 'accretion disks', 'turbulence', 'methods: numerical', 'planets and satellites: formation', 'protoplanetary disks', 'Astrophysics - Earth and Planetary Astrophysics']
Pebbledriven planet formation for TRAPPIST1 and other compact systems
{'Pebble-driven planet formation for TRAPPIST-1 and other compact systems': 'Djoeke Schoonenberg 1 , Beibei Liu 2 , Chris W. Ormel 1 , and Caroline Dorn 3  \n- 1 Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1090 GE Amsterdam, The Netherlands e-mail: d.schoonenberg@uva.nl\n- 2 Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University, Box 43, 221 00 Lund, Sweden\n- 3 University of Zurich, Institute of Computational Sciences, Winterthurerstrasse 190, CH-8057, Zurich, Switzerland  \nJune 18, 2019', 'ABSTRACT': "Recently, seven Earth-sized planets were discovered around the M-dwarf star TRAPPIST-1. Thanks to transit-timing variations, the masses and therefore the bulk densities of the planets have been constrained, suggesting that all TRAPPIST-1 planets are consistent with water mass fractions on the order of 10%. These water fractions, as well as the similar planet masses within the system, constitute strong constraints on the origins of the TRAPPIST-1 system. In a previous work, we outlined a pebble-driven formation scenario. In this paper we investigate this formation scenario in more detail. We used a Lagrangian smooth-particle method to model the growth and drift of pebbles and the conversion of pebbles to planetesimals through the streaming instability. We used the N-body code MERCURY to follow the composition of planetesimals as they grow into protoplanets by merging and accreting pebbles. This code is adapted to account for pebble accretion, type-I migration, and gas drag. In this way, we modelled the entire planet formation process (pertaining to planet masses and compositions, not dynamical configuration). We find that planetesimals form in a single, early phase of streaming instability. The initially narrow annulus of planetesimals outside the snowline quickly broadens due to scattering. Our simulation results confirm that this formation pathway indeed leads to similarly-sized planets and is highly e GLYPH<14> cient in turning pebbles into planets. Our results suggest that the innermost planets in the TRAPPIST-1 system grew mostly by planetesimal accretion at an early time, whereas the outermost planets were initially scattered outwards and grew mostly by pebble accretion. The water content of planets resulting from our simulations is on the order of 10%, and our results predict a 'V-shaped' trend in the planet water fraction with orbital distance: from relatively high (innermost planets) to relatively low (intermediate planets) to relatively high (outermost planets).  \nKey words. accretion, accretion disks - turbulence - methods: numerical - planets and satellites: formation - protoplanetary disks", '1. Introduction': "The cool dwarf star TRAPPIST-1 has been found to be the host of an ultra-compact system of seven Earth-sized planets (Gillon et al. 2016, 2017; Luger et al. 2017). All seven planets are similar in size (about an Earth radius), which is puzzling to explain with classical planet formation theories (Ormel et al. 2017). Moreover, measurements of transit-timing variations and dynamical modelling have constrained their bulk densities, and it has been shown that these densities are consistent with water fractions of a few to tens of mass percent (Grimm et al. 2018; Unterborn et al. 2018b; Dorn et al. 2018). This poses another challenge to planet formation theory. If the planets have formed interior to the water snowline (the distance from the star beyond which water condenses to solid ice), one would expect their water content to be much lower; on the other hand, if they had formed outside of the snowline, one would expect even higher water fractions.  \nOrmel et al. (2017) proposed a scenario for the formation of a compact planetary system like TRAPPIST-1, assuming an inward flux of icy pebbles from the outer disk. In this scenario, planetesimals form in a narrow ring outside the snowline, due to a local enhancement in solids that triggers the streaming instability. The solids enhancement materialises because of outward di GLYPH<11> usion of water vapour and condensation (Stevenson & Lunine 1988; Ciesla & Cuzzi 2006; Ros & Johansen 2013; Schoonenberg & Ormel 2017; Dra¸˙zkowska & Alibert 2017). The resulting  \nicy planetesimals merge and accrete pebbles until they are large enough to start migrating towards the star by type-I migration. Once the protoplanet has migrated across the snowline, the pebbles it accretes are water-poor. The total mass in icy material a protoplanet accretes outside the snowline versus the total mass in dry material it accretes inside the snowline determines its eventual bulk composition. In the case of the TRAPPIST-1 system, this process repeats itself until seven planets are formed, such that the snowline region acts as a 'planet formation factory'. In this formation model, the similar masses are a result of planet growth stalling at the pebble isolation mass, and the moderate water contents are a result of the combination of wet and dry accretion. Unterborn et al. (2018a) have also proposed a formation model for TRAPPIST-1 where inward migration is key to explaining the moderate water fractions, however in their model planet assembly processes are not included.  \nIn this paper we present a numerical follow-up study of Ormel et al. (2017). Dust evolution and planetesimal formation are treated with the Lagrangian smooth-particle method presented in Schoonenberg et al. (2018). The pebble flux from the outer disk is no longer a free parameter as in Ormel et al. (2017), but follows from the dust evolution code. This code is coupled to the N-body code MERCURY (Chambers 1999), which is adapted to account for pebble accretion, gas drag, and type-I migration (Liu et al. 2019). Our model does not follow the later dynamical evo-  \nlution - the re-arrangement of the planetary system architecture during disk dissipation (Liu et al. 2017), or the characteristics and long-term stability of the resonant chain (Papaloizou et al. 2018; Izidoro et al. 2017, 2019) - but it does follow the entire planet formation process from small dust particles to full-sized planets, whilst keeping track of the solid bodies' water content. Any model that aims to explain the origin of the TRAPPIST-1 planets must, of course, be able to explain their observed masses and compositions.  \nIn Sect. 2, we summarise the features of the codes that we employ in this work and describe how they are coupled. We present results of our fiducial model runs in Sect. 3.1 and describe how the results depend on parameter choices in Sect. 3.2. We discuss our results in Sect. 4 and summarise our main findings in Sect. 5.", '2.1. Disk model': 'The gas disk in which planet formation takes place is modelled as a one-dimensional axisymmetric disk, assuming a viscouslyrelaxed (steady-state) gas surface density profile GLYPH<6> gas:  \nGLYPH<6> gas = ˙ M gas 3 GLYPH<25>GLYPH<23> ; (1)  \nwhere ˙ M gas is the gas accretion rate, and the viscosity GLYPH<23> is related to the dimensionless turbulence parameter GLYPH<11> (Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974) as  \nGLYPH<23> = GLYPH<11> c 2 s GLYPH<10> GLYPH<0> 1 ; (2)  \nwith cs the sound speed and GLYPH<10> the Keplerian orbital frequency. The sound speed cs is given by  \ncs = s kBT ( r ) GLYPH<22> ; (3)  \nwhere kB is the Boltzmann constant and GLYPH<22> is the mean molecular weight of the protoplanetary disk gas, which we set to 2.34 times the proton mass. Concerning the temperature T as a function of radial distance to the star r , we consider viscous heating and stellar irradiation, and define T ( r ) as  \nT ( r ) = [ T 4 visc ( r ) + T 4 irr ( r )] 1 = 4 ; (4)  \nwhere T visc( r ) is the viscous temperature profile, which in our standard model is given by  \nT visc( r ) = 180 GLYPH<18> r 0 : 1 au GLYPH<19> GLYPH<0> 1 K ; (5)  \nand T irr( r ) is the irradiation temperature profile, given by  \nT irr( r ) = 150 GLYPH<18> r 0 : 1 au GLYPH<19> GLYPH<0> 1 = 2 K : (6)  \nWe assume a vertically isothermal disk, leading to a disk scale height H gas of  \nH gas = cs = GLYPH<10> : (7)  \nThe gas moves with a velocity v gas = GLYPH<0> 3 GLYPH<23>= 2 r .  \nArticle number, page 2 of 16', '2.2. Dust evolution and planetesimal formation': 'The evolution of dust and formation of planetesimals are treated with a Lagrangian smooth-particle method. Here, we only give a summary of the characteristics of this model. For more information we refer to Schoonenberg et al. (2018).', '2.2.1. Dust evolution': "Weassume that the dust surface density GLYPH<6> dust is initially a radially constant fraction of the gas surface density GLYPH<6> gas. This radially constant fraction is called the metallicity Z , which we set to 2%. Outside of the water snowline, grains consist of 50% water ice and 50% refractory material; interior to the water snowline, the dust composition is completely refractory. Dust grains all start out with the same size (0.1 GLYPH<22> m). At any time, the dust grain size distribution at a given distance from the star is mono-disperse: particles are described by a single size at a given time point and location in the disk. Dust grains grow collisionally assuming perfect sticking (Krijt et al. 2016), and fragment when their mutual impact velocities become larger than the fragmentation threshold velocity, which we set to 10 m / s for ice-coated particles and 3 m / s for refractory particles, motivated by laboratory experiments (Blum & Münch 1993; Sirono 1999; Gundlach & Blum 2015). Although it has recently been reported that the sticking properties of ice vary with temperature (Musiolik & Wurm 2019), for simplicity we keep the fragmentation thresholds constant with the semi-major axis.  \nInitially, when the dust grains are small, the dust is vertically distributed in the same way as the gas; the solids scale height H solids is initially equal to the gas scale height H gas. When dust grains have grown large enough to start decoupling aerodynamically from the surrounding gas, however, they settle to the disk midplane. Vertical settling results in a solids scale height that is smaller than the gas scale height (Youdin & Lithwick 2007):  \nH solids = r GLYPH<11> GLYPH<28> + GLYPH<11> H gas ; (8)  \nwhere GLYPH<28> is the dimensionless stopping time of the particles. In calculating GLYPH<28> , we take into account the composition (water fraction) of the particles, and treat both the Epstein and the Stokes drag regimes. Particles with a non-negligible dimensionless stopping time ( GLYPH<28> & 10 GLYPH<0> 3 ) are called pebbles 1 . Besides settling vertically, pebbles also move radially due to angular momentum loss as a result of experiencing a headwind from the gas. Taking into account the back-reaction of the solid particles onto the gas, the radial velocity of solid particles is given by (Nakagawa et al. 1986)  \nv p = GLYPH<0> 2 GLYPH<17>v K GLYPH<28> GLYPH<0> v gas(1 + GLYPH<24> ) (1 + GLYPH<24> ) 2 + GLYPH<28> 2 ; (9)  \nwhere GLYPH<24> is the midplane solids-to-gas volume density ratio, v K the Keplerian velocity, and GLYPH<17> the 'sub-Keplerianity' or headwind factor,  \nGLYPH<17>v K = GLYPH<0> 1 2 c 2 s v K @ log P @ log r ; (10)  \nwhere P is the gas pressure. The radial pebble flux ˙ M peb is calculated by  \n˙ M peb = GLYPH<6> dust2 GLYPH<25> r j v p j : (11)  \n1 The solids surface density GLYPH<6> dust covers all solid particles. Pebbles belong to the same class of particles as dust grains, which have much smaller stopping times.  \nFig. 1. Cartoon of the code construction used to simulate the planet formation process from dust grains to full-sized planets. The evolution of dust and pebbles and planetesimal formation are treated by the Lagrangian code ('L-code'), which covers the entire disk. Planetesimals form just outside the snowline, and their growth and migration are followed with an N-body code ('N-code'), which accounts for dynamics, planetesimal and pebble accretion, gas drag, and type-I migration. When protoplanets migrate across the inner domain radius of the N-code r domain, their final growth stage is calculated with a semi-analytic model ('A-code'). The way in which the codes are coupled in time is discussed in Sect. 2.4.  \n<!-- image -->  \nThe location of the water snowline depends on the temperature structure, as well as on the water vapour pressure - which in turn depends on the flux of icy pebbles from the outer disk that deliver water vapour to the inner disk (e.g. Ciesla & Cuzzi (2006); Piso et al. (2015); Schoonenberg & Ormel (2017)). We therefore evaluate ˙ M peb outside the water snowline to calculate the position of the water snowline (for more details, see Schoonenberg &Ormel (2017); Schoonenberg et al. (2018)).", '2.2.2. Planetesimal formation': 'The pebble mass flux outside the water snowline also regulates planetesimal formation outside the snowline. The combination of a strong gradient in the water vapour distribution across the water snowline together with some degree of turbulence leads to an outward di GLYPH<11> usive flux of water vapour (Stevenson & Lunine 1988; Ciesla & Cuzzi 2006; Ros & Johansen 2013; Schoonenberg & Ormel 2017; Dra¸˙zkowska & Alibert 2017). The vapour that has been transported outwards condenses onto inward-drifting icy pebbles, leading to a locally enhanced solidsto-gas ratio. This enhancement can be large enough to reach the conditions for streaming instability, such that planetesimals can form in an annulus outside the snowline. We have quantified the critical pebble flux ˙ M peb ; crit needed to trigger streaming instability as a function of GLYPH<28> , GLYPH<11>; and ˙ M gas (Schoonenberg & Ormel 2017; Schoonenberg et al. 2018). If the pebble flux measured outside the water snowline exceeds the critical value ˙ M peb ; crit for a given set of disk conditions, planetesimals form at a rate dM pltsml = dt = ˙ M peb GLYPH<0> ˙ M peb ; crit; the excess pebble flux is transformed to planetesimals. More generally, at any time  \ndM pltsml dt = max[ ˙ M peb GLYPH<0> ˙ M peb ; crit ; 0] : (12)  \nWe assume that all formed planetesimals have the same initial size. In our fiducial model, this initial size is 1200 km, corresponding to a mass of 1 : 8 GLYPH<2> 10 GLYPH<0> 3 Earth masses (for a planetesimal internal density of 1.5 g cm GLYPH<0> 3 , corresponding to a composition of 50% ice and 50% silicates (Schoonenberg et al. 2018)). Simulations show that the initial mass function of planetesimals formed by streaming instability can be described by a power law (possibly with an exponential cuto GLYPH<11> ), such that the typical size of formed planetesimals is a few hundred kilometers in size (Johansen et al. 2015; Simon et al. 2016; Schäfer et al. 2017; Abod et al. 2018), although there is considerable variation amongst these simulations, related to the choice of parameters such as simulation box size and disk metallicity. Our initial planetesimal size of 1200 km is a few factors larger. However, smaller initial planetesimal sizes lead to a larger number of planetesimals and therefore to longer computational timescales. For computational reasons, therefore, we have chosen this value 2 . We have verified that our results do not depend greatly on the initial planetesimal size (Sect. 3.2). Practically, an initial planetesimal size of 1200 km means that a single planetesimal forms every time enough pebble mass has been transformed to planetesimal mass to materialise a new 1200 km-sized planetesimal. The location where this planetesimal forms is assumed to be the location where the solids-to-gas ratio outside the snowline peaks, according to the results of Schoonenberg & Ormel (2017). Planetesimals therefore form in a narrow annulus just exterior to the snowline.', '2.3. Planetesimal growth and migration': 'The planetesimals that form in the Lagrangian dust evolution code described above are followed with an adapted version of the MERCURY N-body code (Chambers 1999; Liu et al. 2019). In this code, planetesimals can grow by planetesimal accretion as well as by pebble accretion. Moreover, gas drag and type-I migration are taken into account. More details can be found in Liu et al. (2019). For planetesimals with the initial size of 1200 km, the migration timescale is still very long. Only after planetesimals have grown larger (by accreting pebbles and other planetesimals) do they start to migrate inwards. Computational times become longer as planetesimals migrate inwards and their orbital periods decrease. To reduce the computational time we therefore take protoplanets out of the N-body code when they cross some (quite arbitrary) distance r domain. We take r domain = 0 : 7 r snow to ensure that the removed bodies are in the water-free growth regime (see Sect. 3.1.1 for further discussion). The last growth phase of protoplanets, which occurs after they have crossed r domain, is treated with a semi-analytic model and is discussed in Sect. 2.5.  \nWe do not take into account the e GLYPH<11> ects on the disk structure of a rapid opacity variation across the snowline - due to the abundance of small grains released by evaporating icy pebbles - leading to a pressure bump (e.g. Kretke & Lin (2007); Bitsch et al. (2014)). A pressure bump - or more generally a local flattening of the pressure profile - would lead to a pile-up of pebbles, thereby aiding planetesimal formation at the snowline. The formation of more icy planetesimals would result in higher planet water fractions, since planet growth by planetesimal accretion would be more e GLYPH<14> cient. The e GLYPH<11> ects of an opacity transition on the migration rate remains unclear, but one possibility is that the snowline acts as a migration trap (e.g. Bitsch et al. (2014); Cridland et al. (2019)). If that were the case, our planet formation model that is based on the inward migration of protoplanets across the snowline would not hold. However, it must be appreciated that these e GLYPH<11> ects heavily depend on the exact thermodynamic state of the disk, which carries major uncertainties (e.g. the e GLYPH<14> ciency of the magneto-rotational instability and the rate at which micron-sized grains coagulate to lower the opacity). For simplicity, therefore, we have sidelined these subtleties and assumed a simple disk structure.', '2.4. Coupling of the two codes': "Wecouple the Lagrangian dust evolution code ('L-code'), the Nbody code ('N-code'), and the semi-analytic 'last growth phase' code ('A-code') in a self-consistent way, enabling us to model the entire planet assembly process (we do not model the longterm dynamical evolution of the planets). The gas disk model (Sect. 2.1) is the same in all codes. Figure 1 provides a cartoon of the simulation setup. The L-code covers the entire disk and describes the evolution of dust and pebbles and planetesimal formation outside the water snowline, as well as the evolution of dry (silicate) pebbles interior to the snowline. The dynamics and growth of planetesimals (by planetesimal and pebble accretion) are followed by the N-code. Interior to r domain, the growth of protoplanets by dry pebble accretion is treated with a semi-analytic model ('A-code'), described in more detail in Sect. 2.5. In this section we describe how the L-code and N-code are coupled.  \nInformation about the pebbles (the pebble flux, as well as the size and composition of pebbles as a function of distance to the star) is provided by the L-code to the N-code. When streaming instability is taking place, the pebble flux that reaches interior to the snowline is reduced because part of the pebble flux out-  \nFig. 2. Sketch of the pebble mass flux as function of semi-major axis in the N-body code. The pebble mass flux decreases with decreasing semimajor axis because of pebble accretion. Interior to the snowline, the pebble mass flux is halved due to evaporation. When streaming instability (SI) is active, the pebble mass flux at the snowline is decreased even further due to planetesimal formation, as dictated by the Lagrangian dust code.  \n<!-- image -->  \nside the water snowline (the excess flux; ˙ M peb GLYPH<0> ˙ M peb ; crit) is converted to planetesimals. Therefore, during streaming instability, the pebble flux reaching inside of the water snowline is given by 0.5 ˙ M peb ; crit, where the factor 0.5 takes into account the evaporation of the pebbles' water content. This is sketched in Fig. 2. Pebbles drift in from the outer disk and are accreted by planetesimals just outside the snowline. The pebble flux therefore decreases with decreasing semi-major axis.  \nThe formation times and locations of planetesimals are provided for the N-code as well. The formed planetesimals are injected into the N-code at their formation times. The locations where planetesimals are injected are picked randomly within an annulus centred on the formation location dictated by the Lcode, but with a finite width of GLYPH<17> r ; which is the typical length scale of streaming instability filaments as found by simulations (Yang & Johansen 2016; Li et al. 2018). The eccentricities and inclinations of the planetesimals are initialised according to a Rayleigh distribution (Liu et al. 2019). When planetesimals (or proto-planets) migrate across the snowline, the characteristics of the pebbles they can accrete change (typically, in our simulations, icy pebbles outside the snowline have dimensionless stopping times GLYPH<28> GLYPH<24> 10 GLYPH<0> 1 , corresponding to a physical size of a few centimetres, and silicate grains interior to the snowline have GLYPH<28> GLYPH<24> 10 GLYPH<0> 3 , corresponding to physical sizes of less than a millimetre); this information is provided to the N-code by the L-code.  \nThe pebble accretion e GLYPH<14> ciency of the planetesimal population outside of the snowline that follows from the N-code is communicated back to the L-code. Here we assume that the fraction of pebbles that are accreted by the already existing planetesimals outside the snowline, as found by the N-code, are not available to form new planetesimals in the L-code. In other words, anywhere exterior to the snowline, pebble accretion happens before planetesimal formation.  \nThe mutual feedback between the two codes compels us to run them in a 'zig-zag' way. To avoid any biases in the code outputs, we switch from one code to the other after fixed time intervals. The length of these time intervals is chosen such that the output values (e.g. pebble mass flux, pebble accretion e GLYPH<14> -  \nFig. 3. Sketch of the code coupling. The code runs in blocks, starting at block 1. For illustration purposes the total number of blocks is 12 in this sketch; in reality, there are many more blocks in one simulation. After some planetesimals have formed during the first run of the Lagrangian dust evolution code (top row), we switch to the N-body code and inject the formed planetesimals (middle row). After a fixed time interval we switch back to the Lagrangian dust code and update the value for the pebble accretion e GLYPH<14> ciency as was output by the N-body code. This 'zig-zag' behaviour continues until all planetesimals have left the N-body code, either because they have been accreted by bigger objects, or because they have grown into bigger objects and have migrated across the inner radius of the N-body domain r domain (see text for more details). The final stage of our calculation concerns the growth of the protoplanets that have crossed r domain, and is discussed in Sect. 2.5.  \n<!-- image -->  \nciency) do not change significantly, but that a significant number of planetesimals are still formed within an interval. We choose a time interval of 5000 years, which meets both these requirements. The coupling of the two codes is illustrated in Fig. 3. The top row depicts the Lagrangian dust evolution and planetesimal formation code; the middle row depicts the N-body code. The code starts at the block indicated by the number 1, then moves to block 2, and so on. The bottom row corresponds to the semianalytic part of the code that deals with the last phase of planet growth, after planetesimals have crossed the inner domain radius of the N-body code. We will discuss this semi-analytic calculation in the next section. Both the Lagrangian code and the Nbody code cover the entire time during which planetesimals are present beyond r domain.", '2.5. Last growth phase': "As discussed in Sect. 2.3, bodies are removed from the N-body code when they migrate interior to r domain, to reduce computational time. After protoplanets have left the N-body code domain, information about their dynamics is lost. Therefore, we do not obtain the eventual planet period ratios, and we make the assumption that the planets leaving the N-body code domain keep their order and do not merge. However, we can still compute their final masses and water fractions, which is what we are ultimately after. In order to do this, we have constructed a semi-analytic model of the very last pebble accretion phase ('Acode'; Fig. 1). During this last growth phase, protoplanets are in the dry region. We assume that the planets are on perfectly circular and planar orbits because we do not have the dynamical information on the planets in the A-code. Assuming zero inclinations is justified because tidal damping timescales are very short (Teyssandier & Terquem 2014; Liu et al. 2019). However, due to resonant forcing, the planetary embryos are expected to have orbits with eccentricities ep on the order of the disk aspect ratio h (Goldreich & Schlichting 2014; Teyssandier & Terquem 2014). We check in Sect. 3.1.6 that our results do not change significantly if we assume ep = h instead of ep = 0 in the pebble accretion e GLYPH<14> ciency calculations for the very inner disk.  \nEach time a protoplanet has crossed r domain, we integrate the mass growth rate of each planet interior to r domain, dM planet ; i = dt , in time, until the next protoplanet crosses r domain and is added to the sequence of protoplanets in the semi-analytic model. The growth rate of protoplanet i is given by  \ndM planet ; i dt = GLYPH<15> PA ; i( t ) ˙ M peb ; rdomain ( t ) j = i GLYPH<0> 1 Y j = 1 1 GLYPH<0> GLYPH<15> PA ; j( t ) ; (13)  \nwhere i 2 [1 ; N ] with N the total number of protoplanets that crossed r domain (an incremental function of time), and where i = 1 and i = N correspond to the outermost and the innermost protoplanet in the A-code, respectively. The symbol ˙ M peb ; rdomain denotes the pebble mass flux at the inner boundary of the Nbody code domain r domain. The value of ˙ M peb ; rdomain follows (as a function of time) from the N-code and the L-code. To calculate the pebble accretion e GLYPH<14> ciency GLYPH<15> PA in the A-code we use the same expression as used in the N-code, which is based on Liu & Ormel (2018) and Ormel & Liu (2018).  \nIn calculating GLYPH<15> PA, we assume that the protoplanets stall at r domain ( GLYPH<24> 0.07 au). We check, however, that putting the planets at random (but ordered) locations between r domain and the current position of the innermost TRAPPIST-1 planet ( GLYPH<24> 0.01 au) during the integration of the growth rates does not matter much for the final outcome, because the pebble accretion e GLYPH<14> ciency does not vary much between 0.01 au and 0.07 au (because the disk aspect ratio h = H gas = r is nearly constant over this distance). Therefore, planetary migration in the final growth phase followed by the A-code is not so important for the final planet masses and compositions.  \nWhen a planet grows large enough to induce a pressure bump exterior to its orbit, the pebble flux is halted: no pebbles can reach the planet's orbit. The mass at which this occurs is called the pebble isolation mass (Morbidelli & Nesvorny 2012; Lambrechts et al. 2014; Bitsch et al. 2018; Ataiee et al. 2018). We do runs both including pebble isolation and without. When we include pebble isolation, the growth of a protoplanet stalls when that protoplanet reaches the pebble isolation mass, and no pebbles reach protoplanets interior to the pebble-isolated planet.  \nTable 1. Fiducial model parameters.  \nNotes. The disk outer radius r out follows from the choices for ˙ M gas, GLYPH<11> , and Mdisk.", '3.1. Fiducial simulation': "We have chosen the fiducial values of our model parameters equal to the values used in Ormel et al. (2017). Our aim is not to perfectly reproduce the TRAPPIST-1 system, for multiple reasons. First of all, the exact compositions of the TRAPPIST-1 planets are uncertain, and we do not account for the possibility of water loss. Due to water loss, simulated planetary systems that we would consider too wet to resemble the TRAPPIST-1 system could become drier systems post formation (see Sect. 4.2 for further discussion). Secondly, computational times are too long to run an optimisation model (e.g. a Markov chain Monte Carlo model), for which many realisations are necessary. Thirdly, it is not trivial to define a quantitative metric for the similarity between synthetic and observed planetary systems (but see Alibert (2019) for a pioneering work). Fourthly, the stochastic component in our model, the injection location of planetesimals in the N-body code (Sect. 2.4), leads to quite some variations between results of simulations with the same initial conditions, as we will demonstrate shortly. Therefore, with 'typical TRAPPIST1 system' we generally mean a compact system of several (515) planets, which have moderate water fractions on the order of 10%.  \nThe input parameters of our model are the gas accretion rate ˙ M gas, the total gas disk mass Mdisk, the disk outer radius r out, the metallicity Z , and the dimensionless turbulence strength GLYPH<11> , four of which are independent. The combination of ˙ M gas and GLYPH<11> defines the gas surface density profile (Eq. (1)). The metallicity of TRAPPIST-1 has been measured to be close to solar (Delrez et al. 2018) and we therefore take Z = 0 : 02. Observed gas accretion rates of low-mass M-dwarf stars are typically on the order of 10 GLYPH<0> 10 solar masses per year (Manara et al. 2017), which is our fiducial value. Observations of gas disk outer radii are notoriously di GLYPH<14> cult. Setting the total gas disk mass to Mdisk = 0 : 038 MTrappist GLYPH<0> 1 leads to an outer disk radius r out = 200 au. The total solids budget of our fiducial disk is Z Mdisk = 18 : 3 M GLYPH<8> , half of which is rocky material. Our fiducial model parameter values are listed in Table 1.", '3.1.1. First growth phase': "In this section we present the results of our fiducial model for the first, predominantly wet growth phase, with which we mean the growth and migration of planetesimals that we follow with the Lagrangian dust evolution code and the N-body code (represented by the two top rows in Fig. 3). The results of the semianalytic procedure for the final, dry growth phase after protoplanets have crossed the inner domain radius of the N-body code (Sect. 2.4) are discussed in Sect. 3.1.2.  \nFig. 4. Protoplanet trajectories (solid lines) for one realisation of the fiducial model. The water snowline (depicted by the dashed blue line) is located at around 0.1 au and moves outwards in time due to a decreasing icy pebble flux. Shaded regions depict background planetesimal densities (the protoplanets are not accounted for in these background planetesimal densities). The first planetesimals form after GLYPH<24> 10 3 years. The colour of the solid lines corresponds to the eccentricity of that protoplanet: the redder, the more eccentric its orbit; the blacker, the more circular. The line width is proportional to protoplanet mass.  \n<!-- image -->  \nIn Fig. 4 we show the evolution trajectories of growing planets for our fiducial simulation parameters. The shading corresponds to the density of planetesimals: the darker a region is shaded, the higher the density of planetesimals in that region. The snowline is depicted by the dashed blue line, and moves outwards in time because the flux of icy pebbles delivering water vapour to the inner disk decreases over time. As in Schoonenberg et al. (2018), we find that planetesimals form outside the snowline at an early time (in the first GLYPH<24> 10 5 years). We note that Fig. 4 gives the impression that planetesimals are present at t = 0. This is because the first planetesimals form after GLYPH<24> 10 3 years, which is close to t = 0 on the timescale plotted. The planetesimal belt that is initially narrow, quickly broadens due to scattering.  \nWe find that not only large protoplanets migrate across the N-body inner domain radius r domain; small planetesimals cross this radius as well due to scattering. The solid lines in Fig. 4 correspond to the largest bodies that amount to more than 99% of the total mass that crosses r domain. These are the bodies for which we follow their subsequent dry growth phase (Sect. 3.1.2). The smallest bodies - corresponding to less than 1% of the total mass that crossed r domain - are ignored in the semi-analytic model of the last growth phase (A-code). In Fig. 5, a histogram shows the masses of bodies that crossed r domain, from a total of ten simulations with the fiducial input parameters. The red line denotes the cumulative mass percentage of bodies that crossed r domain. The vertical dotted line corresponds to the threshold below which bodies correspond to 1% of the total mass that crossed r domain.  \nThe masses of protoplanets are accreted from pebbles as well as from planetesimals. The planetesimal fraction (the mass frac-  \nFig. 5. Histogram of the masses of bodies that crossed the inner domain radius of the N-body simulation r domain. The peak at the smallest mass bin is due to bodies of the initial planetesimal size. The red line denotes the cumulative mass percentage of crossed bodies. The vertical dotted line corresponds to a cumulative mass percentage of 1%. Bodies to the left of this threshold are neglected after they crossed r domain; only the bodies more massive than the threshold mass are followed. This plot stacks the results of ten realisations of the fiducial simulation.  \n<!-- image -->  \n⊕  \ntion contributed by planetesimal accretion) and the water fraction of the total amount of bodies that crossed r domain are plotted as a function of time in Fig. 6. Each line corresponds to one realisation of the fiducial simulation, and each 'crossing event' is denoted with a dot. We find that the general trend is that the planetesimal mass fraction of the mass that crossed r domain (red lines) goes down with time. This is because the later bodies cross the domain radius r domain, the more their growth has been dominated by pebble accretion rather than planetesimal accretion. This is also reflected in Fig. 4, where it is clear that the bodies that migrate inwards at a relatively late time, have been scattered outwards at early times. After having been scattered outwards, these protoplanets find themselves in a planetesimal-free region, and grow mostly by pebble accretion. In contrast, the protoplanets that migrate inwards at an early time have grown in a planetesimal-rich region, and have therefore grown mostly by planetesimal accretion rather than pebble accretion. This is also the reason why in Fig. 6 the water fraction of the total mass in bodies that crossed r domain (black lines) goes down with time. Whereas wet planetesimals were available for the 'early migrators' even interior to the snowline, 'late migrators' could only accrete dry pebbles after migrating past the snowline.  \nFig. 6. Water mass fraction (black) and planetesimal mass fraction (red) of the total mass in protoplanets (and planetesimals) that have crossed the inner domain radius of the N-body simulation r domain, as a function of time. Here we have not yet taken into account the final growth stage in the dry inner disk (A-code). Results of ten simulations with the fiducial parameters are plotted. Each crossing event corresponds to a red and a black dot. Crossing events belonging to the same simulation are connected by a line.  \n<!-- image -->", '3.1.2. Last growth phase': "After protoplanets migrate across the N-body inner domain radius r domain, we select the largest bodies (corresponding to > 99% of the total mass in bodies that crossed r domain, see Fig. 5) and calculate their final growth stage using a semi-analytic model (A-code; Sect. 2.5). In this last growth phase, protoplanets are well interior to the snowline and grow by accreting dry pebbles.  \nIn Fig. 7 we show the results of the semi-analytic computation of this last growth phase for one simulation with the fiducial model parameters. In this calculation we have not taken into account a pebble isolation mass, so for a given protoplanet, pebble accretion continues until the pebble mass flux at that protoplanet's position dries out. In the left panel, the masses of protoplanets at the time of their crossing the N-body inner domain radius r domain are plotted against their crossing time by the grey dots, and their final masses are plotted by the connected black dots. We note that by virtue of our model, the ordering in crossing times is the same as the final ordering in planets: the planet that migrated across r domain first (last) ends up as the innermost (outermost) planet in the system. We see that the first and the third planet are much less massive than the other eight planets, with masses of 0.03 and 0.06 Earth masses, respectively. Their masses at r domain just exceeded our threshold mass (Fig. 5). Due to their relatively low masses, they are less e GLYPH<14> cient at accreting pebbles than the more massive protoplanets, and therefore their masses stay low. We also note that even though the second planet  \n<!-- image -->  \nFig. 7. Properties of simulated planets from one simulation with the fiducial model parameters. Left: Protoplanet masses at the time of their migration across the inner domain radius of the N-body code r domain (grey) and their final masses (black), plotted against the time at which they migrated across r domain. Right: Water mass fractions of protoplanets at the time of their migration across the inner domain radius of the N-body code r domain (grey) and their final water fractions (black), plotted against the time they crossed r domain. In the right panel, the dot size is proportional to the final planet mass.  \n<!-- image -->  \nentered the last growth phase with a higher mass than the fourth planet, the fourth planet eventually becomes more massive than the second planet. This is because the pebble flux at the second planet's position dries out more quickly, due to exterior planets accreting pebbles, than the pebble flux at the fourth planet's position.  \nThe right panel of Fig. 7 shows the water fractions of the planets at the time at which they migrated across the inner domain radius of the N-body simulation r domain (in grey) as well as their final water fractions (in black), plotted against the time they crossed r domain. In this plot the dot sizes are proportional to the final masses of the planets. The water fractions go down during the last growth phase because planets only accrete dry pebbles. The first and the third planet start and end the last, dry growth phase being more water-rich than the other eight planets. This reflects the fact that these planets crossed r domain at relatively low masses, consisting predominantly of wet planetesimals.  \nIgnoring the two smallest planets for now, we find that the innermost and outermost planets generally end up more waterrich than the middle planets, leading to a 'V-shape' in the planet water fraction with orbital distance. With orbital distance (again ignoring the two smallest planets), the planet water fraction starts relatively high at GLYPH<24> 0.19 for planet c, gradually goes down until it reaches its lowest value of GLYPH<24> 0.07 at planet g, beyond which it goes up again to a value of GLYPH<24> 0.17 at planet k. Planet j is an exception to the V-shape, as we discuss below. In Fig. 8 we plot the constituents of each planet (wet planetesimals, wet pebbles, and dry pebbles) for each planet from this simulation. The 'Vshape' is also visible in this figure. The planets that formed earliest (the innermost planets) accreted the most (wet) planetesimals (Sect. 3.1.1; Fig. 6). The later a planet formed, the more its growth was dominated by pebble accretion. The middle planets are therefore generally the driest, because they accreted a lot of dry pebbles after crossing the snowline. The outermost planets also grew mostly by pebble accretion, but by the time they crossed the snowline and could accrete dry pebbles, the pebble flux had already dried out to a great extent. In this simulation, the second-most outer planet does not obey this trend of water fraction with planet order. Because the outermost planet formed  \nFig. 8. Stacked bar chart of planet constituents for each planet of the fiducial model run. Dark blue bars correspond to planetesimals (which are formed outside the snowline and therefore contain 50% water); light blue bars denote wet pebbles (also 50% water) accreted outside the snowline; brown bars depict dry pebbles accreted interior to the snowline. The 'V-shape' in the water fraction with planet order (right panel of Fig. 7) is visible in this figure as well.  \n<!-- image -->  \nrelatively late, the second-outermost planet had quite a large pebble flux to accrete from for quite a long time, without an exterior planet 'stealing' part of the pebbles.", '3.1.3. Comparison to the TRAPPIST-1 system': "In Fig. 9 we compare the simulated planets from the fiducial simulation with the TRAPPIST-1 planets on a mass-radius diagram. The solid, dotted, and dashed lines correspond to the mass-radius relations for a rocky composition with 0%, 10%, and 20% water mass fractions, respectively (Dorn et al. 2018). The TRAPPIST-1 planet data are taken from Dorn et al. (2018). The two smallest simulated 'planets', the innermost ('b') and the third innermost ('d') (see Fig. 7), fall outside the mass domain plotted here (it is conceivable that in reality these small bodies were swallowed by the other, bigger planets). The arrows starting from each simulated planet data point depict the trajectory that each simulated planet would follow on the massradius diagram when its water content evaporated completely (see Sect. 4.2). Earth and Venus are plotted by the black dots and fall on the mass-radius line corresponding to a water fraction of 0%. The water fractions of simulated planets are similar to those of the TRAPPIST-1 planets. The simulated planetary system contains four planets that are slightly more massive than the most massive TRAPPIST-1 planet, but the simulated system agrees with our rough definition of a 'typical TRAPPIST-1 system', which is quite remarkable given our lack of parameter optimisation.", '3.1.4. Multiple realisations': 'In Fig. A.1 the mass-radius diagrams for nine other realisations of the fiducial simulation are shown. Di GLYPH<11> erent realisations of the same simulation do not lead to identical results. The scatter is due to the stochastic component of our model, which is the exact injection location of planetesimals in the N-body code. However, the characteristics of the resulting planetary systems are similar. All realisations of the fiducial simulation lead to a number of planets between 9 and 13 with masses between 0.02 and 3.6 Earth mass. In Table 2 we list the planetary system characteristics, such as the number of planets and the mean planet water fraction, whic result from di GLYPH<11> erent model variations. For each quantity, the mean and the standard deviation are provided, based on the results of di GLYPH<11> erent realisations of each model. The mean and standard deviation of the individual planet mass and water fraction are weighted by the individual planet mass, in order to reduce the importance of very small planets to these quantities (which in reality may not have survived anyway). The results of our fiducial simulation are listed in the first row. The second row provides results for the fiducial model where we have accounted for the pebble isolation mass, as we discuss in Sect. 3.1.5. In Sect. 3.2 we discuss the results of the other model runs, in which we have varied several input parameters.', '3.1.5. Pebble isolation mass': "Our fiducial simulations have produced planets with masses of approximately Earth mass without the implementation of the pebble isolation mass (PIM). Taking the PIM into account imposes a maximum mass on the planets. The PIM has been parametrised as a function of the disk aspect ratio h = H gas = r , the turbulence strength GLYPH<11> , and the stopping time of pebbles GLYPH<28> by using hydrodynamical simulations (e.g. Lambrechts et al. (2014); Bitsch et al. (2018); Ataiee et al. (2018)). Filling in the expressions for the PIM provided by Bitsch et al. (2018) and Ataiee et al. (2018), we find PIM values of 5.0 and 1.8 Earth masses at 0.05 au for our fiducial disk model, respectively. Because these values di GLYPH<11> er by a factor of GLYPH<24> 2, and are larger than  \nthe TRAPPIST-1 planet masses, we here assume a pebble isolation mass of one Earth mass, independent of location 3 and time, with the goal of demonstrating the e GLYPH<11> ect of imposing a maximum planet mass. We take our fiducial simulation results and run the semi-analytical model (A-code), in which we now limit the planet growth to 1 M GLYPH<8> .  \nIn the left panel of Fig. 10, the water fractions of the planets at the time at which they migrated across r domain (in grey) as well as the resulting final water fractions (in black), are plotted against their crossing times. The dot sizes are proportional to the final planet masses. We find that compared to the fiducial simulation without a PIM (Fig. 7), the nine innermost planets are more water-rich. Five of these nine planets have reached the PIM, so that they have accreted fewer dry pebbles than in the case without a PIM, leading to higher final water fractions. The other four planets have been starved of pebbles once an exterior planet became pebble isolated, and have thus also accreted fewer dry pebbles than in the no-PIM case. The outermost planet has not reached PIM nor has it been starved of pebbles by an exterior pebble-isolated planet, and therefore it has the same final water fraction as in Fig. 7. In the right panel of Fig. 10 the mass-radius diagram is presented for the simulation including the PIM. The simulated planets that were too massive to match the TRAPPIST-1 planets in our fiducial model without a pebble isolation mass (Fig. 9) are now limited to one Earth mass. In Table 2 the general results for our simulations including the PIM are listed ('model 2').", '3.1.6. Non-zero eccentricity in very inner disk': 'Due to orbital resonances, the protoplanets in the very inner disk are expected to have non-zero eccentricities, on the order of the disk aspect ratio: ep GLYPH<24> h (Goldreich & Schlichting 2014; Teyssandier & Terquem 2014). We have already discussed that the exact orbital distance of a protoplanet in the very inner disk ( GLYPH<24> 0 : 01-0 : 07 au, covered by the A-code) does not matter much for the pebble accretion e GLYPH<14> ciency, as the disk aspect ratio h is nearly constant over this distance range. Therefore, as long as the ordering of the chain of planetary embryos remains intact, the exact locations of the embryos in the very inner disk are not important for the final planet masses and water fractions. We now check whether our results change when we take ep = h into the calculations of the pebble accretion e GLYPH<14> ciencies in the very inner disk, instead of assuming zero eccentricities in the very inner disk (which is what we have done up to now). In Fig. 11 we show that the results change only slightly if we take ep = h , thereby justifying our assumption that planets in the very inner disk move on circular orbits when calculating their pebble accretion e GLYPH<14> ciencies. The reason why the results with ep = h are very similar to those with ep = 0 is that the relative velocity due to eccentricity ( ep v K ) is not dominant over the Keplerian shear velocity ( GLYPH<24> 0.01 v K ) between the planetary embryos and the pebbles (Liu & Ormel 2018). The pebble accretion e GLYPH<14> ciencies for ep > 0 are always larger than for ep = 0 because we kept the planet inclinations equal to zero (Liu & Ormel 2018).  \nFig. 9. Mass-radius diagram of simulated planets (blue dots) and TRAPPIST-1 planets (red crosses), data taken from Dorn et al. (2018). The solid, dotted, and dashed lines correspond to mass-radius relations for a rocky composition and a water fraction of 0%, 10%, and 20%, respectively. The mass-radius relations are the same as the ones used in Dorn et al. (2018). The arrows starting from the simulated planets depict the trajectories on the mass-radius diagram that the simulated planets would cover if they were to lose their entire water content.  \n<!-- image -->  \n⊕', '3.2. Dependence on model parameters': "In this section we discuss the outcomes of simulations in which we vary one input parameter at a time. The results are listed in Table 2. In all models except 'model 9' and 'model 10', we keep the initial gas and solids disk masses fixed at the fiducial values (Table 1), by adjusting the disk size accordingly.  \nWe first check the dependence of our results on the initial planetesimal size, the fiducial value of which is 1200 km. Decreasing the initial planetesimal size to 1000 km ('model 4') does not change the results, which justifies our fiducial choice of 1200 km. Increasing the initial size to 1500 km ('model 3') leads to slightly fewer and larger planets, but the total mass in planets remains unchanged.  \nA lower value of the turbulence strength GLYPH<11> implies a smaller disk, because the gas surface density GLYPH<6> gas is inversely proportional to GLYPH<11> in our framework of a viscous gas disk. A lower GLYPH<11> -value also implies a larger pebble flux, because the dust surface density is initially a constant fraction of the gas surface density. Because we only vary GLYPH<11> and keep the gas accretion rate ˙ M gas fixed, the pebble-to-gas flux ratio is also larger for smaller GLYPH<11> . Planetesimal formation by streaming instability outside the snowline is strongly dependent on the pebble-to-gas mass flux  \nratio (Schoonenberg & Ormel 2017). It is much less dependent on the disk size, because solid material from the very outer parts of the disk does not take part in the streaming instability outside the snowline, which happens at an early stage of the evolution of the disk. Therefore, as was already noted in Schoonenberg et al. (2018), even though we keep the solids budget fixed, a smaller GLYPH<11> leads to the formation of more planetesimals (disregarding secondary e GLYPH<11> ects such as the subtle dependencies on GLYPH<11> of the pebble size and of the streaming instability threshold on the pebble flux). Indeed, decreasing the value of GLYPH<11> to 5 GLYPH<2> 10 GLYPH<0> 4 ('model 5') results in more mass in formed planetesimals MSI. We note that although we lowered GLYPH<11> by a factor of 2, the mass in formed planetesimals increased by a larger factor. This is because the planetesimal formation rate is proportional to the pebble flux ˙ M peb that is in excess of the critical pebble mass flux ˙ M peb ; crit (Eq. (12)). Increasing the pebble flux by a factor of f therefore results in the planetesimal formation rate being increased by a factor greater than f (Schoonenberg et al. 2018). In the case of GLYPH<11> = 5 GLYPH<2> 10 GLYPH<0> 4 , most of the icy pebble flux is converted to planetesimals. Due to computational limitations, we performed the simulations of this model variation with the larger initial planetesimal size of 1500 km, in order to curb the number of planetesimals in the N-body code. To get rid of the planetesimal-sized bodies  \n<!-- image -->  \n<!-- image -->  \n⊕  \nFig. 10. Properties of planets from the fiducial simulation including a pebble isolation mass of 1 M GLYPH<8> . Left: Same as the right panel of Fig. 7, but now including a pebble isolation mass of 1 M GLYPH<8> . Right: Same as Fig. 9, but now including a pebble isolation mass of 1 M GLYPH<8> .  \n<!-- image -->  \nFig. 11. Comparison of planet properties from the fiducial simulation and the fiducial simulation with non-zero protoplanet eccentricities. Left: Final planet masses against the time at which they migrated across r domain for the fiducial simulation (black), and for the fiducial simulation accounting for non-zero protoplanet eccentricities ep in the very inner disk (by setting ep equal to the disk aspect ratio h in the A-code) (red). Right: Same as left panel, but now for final planet water fractions instead of masses. Taking into account non-zero eccentricities in the very inner disk leads to only small variations in the final results.  \n<!-- image -->  \nthat were scattered interior to the N-body domain, before focusing on the protoplanets in the semi-analytical calculation of the last growth phase (A-code), we had to neglect the smallest 3% instead of 1%. On average, 45.8 planets form in these simulations, with an average total mass in planets of 14.8 M GLYPH<8> . This is a lot more than in the fiducial model runs, because in this case there are many more planetesimal seeds to grow planets from, and pebble accretion is more e GLYPH<14> cient for lower GLYPH<11> . The type-I migration timescale is inversely proportional to the gas surface density (Tanaka et al. 2002), and is therefore shorter for lower GLYPH<11> , leading to faster migration across the snowline. Additionally, planetesimals form closer to the snowline for lower GLYPH<11> (Schoonenberg & Ormel 2017). From these considerations one would expect lower final water fractions for lower GLYPH<11> , because protoplanets outside the snowline have to traverse a shorter distance at a larger migration velocity to reach the inner disk in which they can accrete dry pebbles. However, the e GLYPH<11> ect of the formation of many more wet planetesimals that accrete each other is  \nmore dominant, and we find that the planets resulting from the simulations with lower GLYPH<11> have much higher water fractions and are smaller than those produced in the fiducial model 4 .  \nA higher value of GLYPH<11> , on the other hand, implies a larger disk and a smaller pebble flux that is sustained for a longer time. Because planetesimal formation happens only in the early stages of the disk evolution, one would expect fewer planetesimals to form in this case, due to the smaller pebble flux. However, increasing the value of GLYPH<11> to 2 GLYPH<2> 10 GLYPH<0> 3 ('model 6') results in slightly more planetesimals compared to the fiducial model. This is because although the pebble flux is smaller, the streaming instability phase lasts longer, due to the subtle dependency of the critical pebble mass flux (needed to trigger the streaming instabil-  \nTable 2. Simulation results for di GLYPH<11> erent models.  \nNotes. The columns denote (from left to right): model description; total number of planets; total mass in planets; total mass in planetesimals formed by streaming instability; weighted mean of individual planet mass; mean planet water fraction (weighted with planet mass). Each column entry gives the mean and standard deviation of the corresponding quantity, calculated from multiple realisations. The model descriptions are elaborated on in the text.  \nity) on the balance between di GLYPH<11> usion, vertical settling, and radial drift, and therefore on the combination of GLYPH<11> and the dimensionless stopping time of pebbles GLYPH<28> (Schoonenberg & Ormel 2017; Schoonenberg et al. 2018). The migration timescale is longer for higher GLYPH<11> , and planetesimals are formed further away from the snowline. Therefore, protoplanets reside exterior to the snowline for a longer time compared to the fiducial model. These e GLYPH<11> ects result in slightly higher final water fractions compared with the fiducial model. The final planet masses are a bit lower than in the fiducial model, because pebble accretion is less e GLYPH<14> cient for higher GLYPH<11> . For GLYPH<11> values higher than 5 GLYPH<2> 10 GLYPH<0> 3 , the region outside the snowline does not reach the conditions for streaming instability and no planets are formed, but this could be resolved by changing other model parameters, such as the metallicity, as well.  \nVarying the gas accretion rate ˙ M gas has similar e GLYPH<11> ects as varying GLYPH<11> . Lower (higher) ˙ M gas means a larger (smaller) disk with lower (higher) surface density and pebble flux. However, when we change ˙ M gas, we also change the viscous temperature profile (Eq. (5)), such that the snowline is located further outwards (inwards) for higher (lower) ˙ M gas. The pebble-to-gas mass flux ratio is inversely related to the snowline location (Schoonenberg et al. 2018). The planetesimal formation rate in turn depends on the pebble-to-gas mass flux ratio, but not linearly: we remember that only the pebble flux in excess of the critical pebble flux is converted to planetesimals (Eq. (12)). With all these complications, we find that reducing the gas accretion rate to ˙ M gas = 5 GLYPH<2> 10 GLYPH<0> 11 M GLYPH<12> yr GLYPH<0> 1 results in approximately twice as much mass in planetesimals formed by streaming instability. As in the case of lower GLYPH<11> , this leads to a larger number of planets than in the fiducial case, because there are more planetesimal seeds from which to grow. The larger amount of wet planetesimals overshadows the e GLYPH<11> ect of lower migration velocity (which tends to increase the planet water fraction), and the final water fractions are a bit higher than in the fiducial model. A higher gas accretion rate of ˙ M gas = 2 GLYPH<2> 10 GLYPH<0> 10 M GLYPH<12> yr GLYPH<0> 1 in turn leads to an amount of formed planetesimals similar to the fiducial case. Because in this model the type-I migration timescale is shorter than in the fiducial model, the planets resulting from this model are slightly more numerous, less massive, and more water-rich than in the fiducial case.  \nIn our model, the parameter values are fixed for the duration of the simulation. If in reality protoplanetary disks feature de-  \ncreasing values of GLYPH<11> during their evolution, we expect that planetesimal formation outside the snowline starts when the conditions for streaming instability are marginally reached, such that only a small fraction of the pebble flux is converted to planetesimals. Therefore, our fiducial model may be more physical than the models with lower and higher GLYPH<11> values, in which larger fractions of the icy pebble flux were converted to planetesimals and planets formed with higher water fractions.  \nIn 'model 9' and 'model 10' we change the disk outer radius r out to 300 au and 100 au, respectively. The gas disk mass changes accordingly, in contrast to all models discussed up to now, in which we kept the gas disk mass fixed. Planetesimal formation takes place in the early stage of the disk evolution and is therefore not a GLYPH<11> ected by these changes in the disk outer radius: the pebbles from which planetesimals are formed originate at distances to the star closer than 100 au. However, pebble accretion can go on for a longer (shorter) time if the disk is larger (smaller). This leads to more massive (less massive) planets for larger (smaller) r out compared with the fiducial model. The water fractions for smaller r out are also higher than in the fiducial model, because of less dry pebble accretion in the inner disk. We note that the opposite is not true: the mean water fraction for larger r out is only slightly lower than that in the fiducial model. This is because a few planets 'steal' the long-lasting flux of dry pebbles and become relatively water-poor, leaving the other planets without many dry pebbles to accrete; these other planets therefore stay water-rich. We note that a model with a smaller disk size may be more realistic than our fiducial model, as observations seem to suggest that disks around very low-mass stars are quite small ( < 100 au) (Ricci et al. 2013; Hendler et al. 2017), however, primordial disk radii may have been larger (e.g. Rosotti et al. (2019)).  \nAlthough the number of planets and their final masses and water fractions depend on our model parameters, all planetary systems resulting from our di GLYPH<11> erent model variations feature the same trend in the planet water fraction with planet order: the inner- and outermost planets are generally more water-rich than the middle planets, for the reasons explained in Sect. 3.1.2. For a very large disk (as in 'model 9'), the uprise in water fraction in the outermost planets is less strong, because in that case dry pebble accretion is still important even for the outermost planets.", '4.1. Evolution of the snowline location': 'In our model, the water snowline location depends on the temperature and vapour pressure. The vapour pressure depends on the influx of icy pebbles delivering water vapour to the inner disk, and therefore varies with time. The pebble flux goes down with time, and the snowline therefore moves outwards with time (Fig. 4). We have assumed that the gas accretion rate ˙ M gas and therefore the viscous temperature profile are constant in time. However, taking into account that the gas accretion rate decreases with time leads to an inward movement of the snowline (Garaud & Lin 2007; Oka et al. 2011; Martin & Livio 2012). Schoonenberg et al. (2018) have shown that the planetesimal formation phase is not a GLYPH<11> ected by the migration of the snowline, because the viscous evolution timescale is much longer than the duration of the planetesimal formation phase. The general picture of planets acquiring moderate water fractions by originating and growing outside the snowline and then migrating inwards would also not be a GLYPH<11> ected by a moving snowline, as long as the early growth and migration of planets is faster. Starting with a single embryo of 0.01 Earth masses and a simplified planetary growth function, Bitsch et al. (2019) calculated planetary growth tracks and resulting planetary compositions in an evolving disk, in which the snowline moves inwards with time. When they only include inward migration, they find that embryos indeed grow and migrate inwards faster than the evolution of the snowline. However, as they mention, the eventual composition of a planet does depend on the growth rate (pebble flux). For example, one can imagine that late-stage pebble accretion might be a GLYPH<11> ected by snowline migration, if the pebble flux coming from the outer disk remains considerable on the viscous evolution timescale. The outermost planets - that migrate inwards at late times, when the snowline might have moved inwards considerably - could then be surrounded by icy pebbles rather than dry pebbles for a slightly longer time. This would lead to even more water-rich outer planets and therefore to a strengthening of the trend in the water fraction with planet ordering (inner- and outermost planets being more water-rich than middle planets; see Sect. 3.1.2).', '4.2. Water loss': "In our simulations we have not taken into account physical processes by which water may be lost from the planets during or after their formation. When protoplanets migrate across the water snowline, their surface escape velocity is much larger than the thermal motions of a water vapour atmosphere (this is not true only for bodies smaller than a few hundred kilometres). We therefore do not expect protoplanets to lose water when they migrate past the water snowline. Collisions may present a bigger threat to water retention. We assume each collision between protoplanets and / or planetesimals leads to perfect merging, not taking into account the possibility of more destructive outcomes. According to Burger et al. (2018), up to 75% of the water can be lost in a hit-and-run collision between Mars-sized bodies (in a gas-free environment), and according to Marcus et al. (2010), an entire icy mantle can be stripped o GLYPH<11> a rocky core by a violent giant impact. If violent collisions happen outside the snowline, a large fraction of the 'lost' water may still end up in the planets via condensation on pebbles, but collisions interior to the snowline would lead to lower planet water fractions.  \nBolmont et al. (2017) and more recently Bourrier et al. (2017) have estimated the amount of water loss from the  \nTRAPPIST-1 planets due to photodissociation as a result of high-energy stellar radiation after the disk disappeared. Extrapolating the current high-energy stellar irradiation back in time using an evolutionary model, Bourrier et al. (2017) report an upper limit of tens of Earth oceans on the water loss of TRAPPIST-1g since its formation. The other planets have lost lesser amounts. Tens of Earth oceans corresponds to a few percent of an Earth mass, which is still smaller than the water contents of our simulated planets and of the present-day TRAPPIST-1 planets (Grimm et al. 2018; Dorn et al. 2018).  \nIf a significant fraction of the TRAPPIST-1 planets' water contents were lost during or post formation, a disk model in which planets form more water-richly than the currently observed TRAPPIST-1 planets may be more realistic (e.g. a smaller, and therefore less massive disk than our fiducial model, such as 'model 10', see Table 2).", '4.3. Planetary systems other than TRAPPIST-1': "We have shown that compact systems of Earth-sized planets can form around M-dwarf stars under di GLYPH<11> erent disk conditions. The occurrence rate of small planets is higher around M-dwarfs than around FGK-stars (Mulders et al. 2015), but more systems like TRAPPIST-1, with as many as seven small planets, have not been discovered yet. However, forming a large enough number of planets has proven to be no problem in our model, which leaves us wondering why no other compact 'many-planet' systems have been discovered. One reason could be that a large fraction of compact multi-planet systems originally formed with more planets in a resonance chain, but became unstable (Izidoro et al. 2017; Lambrechts et al. 2019; Izidoro et al. 2019).  \nAcrucial feature of our model is that the icy pebble flux from the outer disk reaches the snowline such that the streaming instability outside the snowline can be triggered. If a massive planet outside the snowline were to open a gap and halt the pebble flux, planetesimal formation would shut down. We hypothesise that this is the main reason for the di GLYPH<11> erent architecture of compact multi-planet systems like TRAPPIST-1 (in which all planets have similar sizes) on the one hand, and a system like the solar system on the other. In the solar system, the early formation of Jupiter may have stalled icy planetesimal formation just exterior to the snowline, leaving only dry material available to form the rocky planets in the inner solar system (Morbidelli et al. 2015; Lambrechts et al. 2019).", '5. Conclusions': "We have performed a numerical follow-up study regarding the formation model for compact planetary systems around M-dwarf stars, such as the TRAPPIST-1 system, which was presented in Ormel et al. (2017). We have coupled a Lagrangian dust evolution and planetesimal formation code (Schoonenberg et al. 2018) to an N-body code adapted to account for gas drag, type-I migration, and pebble accretion (Liu et al. 2019). This coupling enabled us, for the first time, to self-consistently model the planet formation process from small dust grains to full-sized planets, whilst keeping track of their water content. The main conclusions from this work can be summarised as follows:  \n- 1. Our work has demonstrated that the planet formation scenario proposed in Ormel et al. (2017) indeed produces multiple Earth-sized planets (even without taking into account the pebble isolation mass), with water mass fractions on the order of 10%.  \n- 2. In contrast to what was hypothesised in Ormel et al. (2017), we find that all planetesimals form in one, early streaming instability phase.\n- 3. Because the N-body code contains a stochastic component (the exact injection location of planetesimals), we find scatter in the characteristics of planetary systems (such as the number of planets, the average planet mass, and water fraction) resulting from simulations with the same input parameters. The amount of scatter is quantified in Table 2.\n- 4. Our planet formation model shows a universal trend in the water fraction with planet order: inner- and outermost planets are generally more water-rich than middle planets, independent of the model parameters. This 'V-shape' in the variation of the planet water fraction with orbital distance emerges because the three di GLYPH<11> erent growth mechanisms (planetesimal accretion, wet pebble accretion, and dry pebble accretion) have di GLYPH<11> erent relative importance for the inner, middle, and outer planets. The importance of wet planetesimal accretion for planet growth diminishes with planet order, whereas pebble accretion becomes more important for planets with larger orbital distance. The e GLYPH<11> ect is a decreasing water fraction with planet orbital distance. However, the ratio of wet pebble accretion to dry pebble accretion goes down with planet orbital distance, such that the water fraction goes up again for the outermost planets.  \nWe reiterate that, due to computational limitations, in this work we have solely focused on the final masses and compositions of the planets, and have not followed the final dynamical configuration of planetary systems. In order to extend our formation model such that it also treats the long-term dynamics of the system, the N-body code should be extended into the very inner region of the disk (interior to r domain), or a dedicated dynamics study should be performed of this inner region, by making use of the results of the simulations presented in this paper.  \nAcknowledgements. D.S. and C.W.O are supported by the Netherlands Organization for Scientific Research (NWO; VIDI project 639.042.422). B.L. thanks the support of the European Research Council (ERC Consolidator Grant 724687PLANETESYS), and the Swedish Walter Gyllenberg Foundation. C.D. acknowledges the support of the Swiss National Foundation under grant PZ00P2\\_174028 and the Swiss National Center for Competence in Research PlanetS . We thank the anonymous referee for his or her comments that helped improve the manuscript.", 'Appendix A: Mass-radius diagrams for multiple realisations of the fiducial model.': 'Weattach as an appendix the mass-radius diagrams of the planets resulting from nine realisations of the fiducial simulation.', 'References': 'Abod, C. P., Simon, J. B., Li, R., et al. 2018, arXiv e-prints, arXiv:1810.10018  \nAlibert, Y. 2019, arXiv e-prints, arXiv:1901.09719 Ataiee, S., Baruteau, C., Alibert, Y., & Benz, W. 2018, A&A, 615, A110 Bitsch, B., Morbidelli, A., Johansen, A., et al. 2018, A&A, 612, A30 Bitsch, B., Morbidelli, A., Lega, E., Kretke, K., & Crida, A. 2014, A&A, 570, A75 Bitsch, B., Raymond, S. N., & Izidoro, A. 2019, A&A, 624, A109 Blum, J. & Münch, M. 1993, Icarus, 106, 151 Bolmont, E., Selsis, F., Owen, J. E., et al. 2017, MNRAS, 464, 3728 Bourrier, V., de Wit, J., Bolmont, E., et al. 2017, AJ, 154, 121 Burger, C., Maindl, T. I., & Schäfer, C. M. 2018, Celestial Mechanics and Dynamical Astronomy, 130, 2 Chambers, J. E. 1999, MNRAS, 304, 793 Ciesla, F. J. & Cuzzi, J. N. 2006, Icarus, 181, 178 Cridland, A. J., Pudritz, R. E., & Alessi, M. 2019, MNRAS, 484, 345 Delrez, L., Gillon, M., Triaud, A. H. M. J., et al. 2018, MNRAS, 475, 3577 Dorn, C., Mosegaard, K., Grimm, S. L., & Alibert, Y. 2018, ApJ, 865, 20 Dra¸˙zkowska, J. & Alibert, Y. 2017, A&A, 608, A92 Garaud, P. & Lin, D. N. C. 2007, ApJ, 654, 606 Gillon, M., Jehin, E., Lederer, S. M., et al. 2016, Nature, 533, 221 Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456 Goldreich, P. & Schlichting, H. E. 2014, AJ, 147, 32 Grimm, S. L., Demory, B.-O., Gillon, M., et al. 2018, A&A, 613, A68 Gundlach, B. & Blum, J. 2015, ApJ, 798, 34 Hendler, N. P., Mulders, G. D., Pascucci, I., et al. 2017, ApJ, 841, 116 Izidoro, A., Bitsch, B., Raymond, S. N., et al. 2019, arXiv e-prints, arXiv:1902.08772 Izidoro, A., Ogihara, M., Raymond, S. N., et al. 2017, MNRAS, 470, 1750 Johansen, A., Mac Low, M.-M., Lacerda, P., & Bizzarro, M. 2015, Science Advances, 1, 1500109 Kretke, K. A. & Lin, D. N. C. 2007, ApJ, 664, L55 Krijt, S., Ormel, C. W., Dominik, C., & Tielens, A. G. G. M. 2016, A&A, 586, A20 Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35 Lambrechts, M., Morbidelli, A., Jacobson, S. A., et al. 2019, arXiv e-prints, arXiv:1902.08694 Li, R., Youdin, A. N., & Simon, J. B. 2018, ApJ, 862, 14 Liu, B. & Ormel, C. W. 2018, A&A, 615, A138 Liu, B., Ormel, C. W., & Johansen, A. 2019, arXiv e-prints, arXiv:1902.10062 Liu, B., Ormel, C. W., & Lin, D. N. C. 2017, A&A, 601, A15 Luger, R., Sestovic, M., Kruse, E., et al. 2017, Nature Astronomy, 1, 0129 Lynden-Bell, D. & Pringle, J. E. 1974, MNRAS, 168, 603 Manara, C. F., Testi, L., Herczeg, G. J., et al. 2017, A&A, 604, A127 Marcus, R. A., Sasselov, D., Stewart, S. T., & Hernquist, L. 2010, ApJ, 719, L45 Martin, R. G. & Livio, M. 2012, MNRAS, 425, L6 Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icarus, 258, 418 Morbidelli, A. & Nesvorny, D. 2012, A&A, 546, A18 Mulders, G. D., Pascucci, I., & Apai, D. 2015, ApJ, 814, 130 Musiolik, G. & Wurm, G. 2019, ApJ, 873, 58 Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375 Oka, A., Nakamoto, T., & Ida, S. 2011, ApJ, 738, 141 Ormel, C. W. & Liu, B. 2018, A&A, 615, A178 Ormel, C. W., Liu, B., & Schoonenberg, D. 2017, A&A, 604, A1 Papaloizou, J. C. B., Szuszkiewicz, E., & Terquem, C. 2018, MNRAS, 476, 5032 Piso, A.-M. A., Öberg, K. I., Birnstiel, T., & Murray-Clay, R. A. 2015, ApJ, 815, 109 Ricci, L., Isella, A., Carpenter, J. M., & Testi, L. 2013, ApJ, 764, L27 Ros, K. & Johansen, A. 2013, A&A, 552, A137 Rosotti, G. P., Tazzari, M., Booth, R. A., et al. 2019, MNRAS, 486, 4829 Schäfer, U., Yang, C.-C., & Johansen, A. 2017, A&A, 597, A69 Schoonenberg, D. & Ormel, C. W. 2017, A&A, 602, A21 Schoonenberg, D., Ormel, C. W., & Krijt, S. 2018, A&A, 620, A134 Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337 Simon, J. B., Armitage, P. J., Li, R., & Youdin, A. N. 2016, ApJ, 822, 55 Sirono, S. 1999, A&A, 347, 720 Stevenson, D. J. & Lunine, J. I. 1988, Icarus, 75, 146 Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257 Teyssandier, J. & Terquem, C. 2014, MNRAS, 443, 568 Unterborn, C. T., Desch, S. J., Hinkel, N. R., & Lorenzo, A. 2018a, Nature Astronomy, 2, 297 Unterborn, C. T., Hinkel, N. R., & Desch, S. J. 2018b, Research Notes of the American Astronomical Society, 2, 116 Yang, C.-C. & Johansen, A. 2016, The Astrophysical Journal Supplement Series, 224, 39  \nYoudin, A. N. & Lithwick, Y. 2007, Icarus, 192, 588  \n<!-- image -->  \nArticle number, page 16 of 16'}