{'fraction_non_alphanumeric': 0.03842162116366255, 'fraction_numerical': 0.02414589202126533, 'mean_word_length': 4.258269295021718, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Formation of multiple-planet systems in resonant chains around M dwarfs

Recent observations have revealed the existence of multiple-planet systems composed of Earth-mass planets around late M dwarfs. Most of their orbits are close to commensurabilities, which suggests that planets were commonly trapped in resonant chains in their formation around low-mass stars. We investigate the formation of multiple-planet systems in resonant chains around low-mass stars. A time-evolution model of the multiple-planet formation via pebble accretion in the early phase of the disk evolution is constructed based on the formation model for the TRAPPIST-1 system by Ormel et al. (2017). Our simulations show that knowing the protoplanet appearance timescale is important for determining the number of planets and their trapped resonances: as the protoplanet appearance timescale increases, fewer planets are formed, which are trapped in more widely separated resonances. We find that there is a range of the protoplanet appearance timescale for forming the stable multiple-planet systems in resonant chains. This range depends on the stellar mass and disk size. We suggest that the protoplanet appearance timescale is a key parameter for studying the formation of multiple-planet systems with planets in resonant chains around low-mass stars. The composition of the planets in our model is also discussed.


INTRODUCTION
Recent observations have revealed the existence of multiple-planet systems around low-mass stars, such as TRAPPIST-1 (Gillon et al. 2017;Luger et al. 2017), YZ Ceti (Astudillo-Defru et al. 2017;Stock et al. 2020), Teegarden's star (Zechmeister et al. 2019), and GJ 1061 (Dreizler et al. 2020). These planets have approximately one Earth mass and are orbiting around stars with masses of ∼ 0.1M , where M is the solar mass. Except for those of the planets orbiting around Teegarden's star, their period ratios are close to the period commensurabilities. Although more multiple-planet systems around low-mass stars are needed to perform statistical analyses, the fraction (3/4) of resonant or near resonant planetary systems around low-mass stars suggests that planets in or near resonant chains around low-mass stars would be more abundant than those around stars with masses of ∼ 1M (Fabrycky et al. 2014;Winn & Fabrycky 2015). Ormel et al. (2017) proposed a scenario for the formation of compact planetary systems around TRAPPIST-1, which has about 0.08M . This scenario includes the following stages: a planetesimal forms around the H 2 O iceline (hereafter iceline) owing to the streaming instability; the planetesimal grows through pebble accretion to a protoplanet with the isolation mass while migrating inward; the protoplanet stops migration around the magnetospheric cavity of the disk; the next few protoplanets are formed sequentially in the similar manner and are trapped in resonant chains; after the disk gas depletion, some mechanisms move planets out of resonance; finally, planets reside at the current orbits. This analytical model was further investigated by a numerical study (Schoonenberg et al. 2019). Furthermore, population synthesis simulations showed the mass distribution of protoplanets around low-mass stars based on the similar concept (Liu et al. 2019a(Liu et al. , 2020Dash & Miguel 2020).
The aforementioned studies provide us insights into the accretion of protoplanets around M dwarfs. However, the formation of multiple planets in the chains of resonant orbits in the pebble-driven core accretion scenario has not been fully understood. In particular, the configuration of planets in resonant chains would be affected by the planetesimal forming efficiency. Recent studies showed that the streaming instability can be either promoted or hampered by turbulent stirring (Chen & Lin 2020;Gole et al. 2020;Schäfer et al. 2020;Umurhan et al. 2020). The efficiency of the planetesimal formation through the streaming stability is also affected by the size distribution of dust (e.g., Bai & Stone 2010;Krapp et al. 2019). These uncertainties may lead to a wide range of the planetesimal formation time. The timescale of the subsequent growth of planetesimals to protoplanets depends on the uncertain size distribution of planetesimals formed via the streaming stability (e.g., Liu et al. 2019b). It is expected that a longer appearance interval of a protoplanet that grows to be a planet leads to systems composed of a smaller number of planets in more widely separated resonances, which would not cause an orbital instability (Matsumoto et al. 2012).
In this study, we investigate the formation of multipleprotoplanet systems in which the protoplanets are trapped in resonant chains. The time-sequential evolution of the protoplanets is considered based on the pebble-driven core accretion scenario Schoonenberg et al. 2019). Moreover, the condition of the resonant trapping of protoplanets growing through pebble accretion is presented as a function of the protoplanet appearance time, stellar mass, and disk size. The protoplanets trapped in resonant chains are referred to as "planets" in this paper. Our simulation results also show in which resonant chain the planets are trapped.
The remainder of this paper is organized as follows. The model is described in Section 2. We present the results of our analytical simulations in Section 3. Section 4 discusses the influences of our assumptions on the results. Finally, the conclusion is summarized in Section 5.

Outline
We basically follow the pebble-driven core accretion scenario proposed in Ormel et al. (2017). We consider a protoplanetary disk composed of gas and dust. First, icy dust grains grow through coagulation. As their masses increase, the grains become gradually decoupled from the gas according to their increasing Stokes number τ p , which is the dimensionless stopping time. We assume that the grains begin to drift inward when τ p = 0.05 (e.g., Birnstiel et al. 2012;Okuzumi et al. 2012). These grains are called "pebbles" in this study. Their growth timescale t peb is a function of the dust-to-gas mass ratio (Birnstiel et al. 2012;Okuzumi et al. 2012); t peb = ξ(Z 0 Ω Kep ) −1 is the timescale over which the pebbles begin to drift inward, where the global dust-to-gas mass ratio Z 0 is 0.02, Ω Kep is the Kepler frequency, and ξ = 10 . The radial drift of pebbles ensues from the inside to outside of the disk. The first icy pebbles form at 500/(2π)T Kep,ice 80T Kep,ice and the final icy pebbles form at t peb,out 80T Kep,disk , where T Kep is the Kepler time; T Kep,ice is T Kep at the iceline and T Kep,disk is T Kep at the disk outermost radius (R disk ).
The pebbles drift inward and create a pileup structure near the iceline (e.g., , which is simply described by a Gaussian function of the pebble flux in our model. The pebble flux in the region exterior to the iceline is derived from the growth and drift of pebbles. When the dust density at the pileup peak exceeds the gas density, planetesimals are formed via the streaming instability. An embryo is formed through the pebble accretion of scattered planetesimals and/or via the runaway growth of planetsimals (Liu et al. 2019b;Schoonenberg et al. 2019). We focus on the growth of the embryo, which is the precursor of a (proto)planet. We use t pl to express the appearance timescale of a 100 km sized embryo through the streaming instability and subsequent growth. Hereafter, we call embryos as protoplanets for simplicity. The protoplanet accretes the surrounding pebbles (e.g., Ormel & Klahr 2010;Lambrechts & Johansen 2012). After the protoplanet has gained sufficient mass, it migrates inward owing to the protoplanet-disk tidal interactions (e.g., Ward 1986;Lin & Papaloizou 1993). The mass growth of the protoplanet through pebble accretion lasts until its mass reaches the pebble isolation mass (Lambrechts et al. 2014). Afterward, the migration of the first protoplanet is stopped at the magnetospheric cavity of the disk. The following protoplanets are sequentially trapped in mean motion resonances based on the competition between the resonant libration timescale and the migration timescale through the resonant width (Ogihara & Kobayashi 2013).
We perform simulations of the protoplanets in this growth scenario by changing their appearance time (t pl ), stellar mass (M ), and disk size (R disk ). Our simulations provide the number of planets trapped in resonant chains and the combinations of various resonant chains in a planetary system, which affect their subsequent orbital evolutions (i.e., whether they cause orbital instabilities or not (Matsumoto et al. 2012;Matsumoto & Ogihara 2020)). The details of the proposed model are presented in the following subsections.

Disk
The disks are divided into the two regions because the pebble formation timescales ( 80T Kep ) are shorter than the viscous timescales . The inner disk, which is the region interior to the iceline, is viscously relaxed. In this region, the gas surface density is given by Σ g,in =Ṁ g /3πν (Lynden-Bell & Pringle 1974), where the viscosity is ν = α(hr) 2 Ω Kep , α = 10 −3 in our model, h is the disk aspect ratio, and r is the orbital radius (Shakura & Sunyaev 1973). The gas accretion rate depends on the stellar mass (Manara et al. 2015):Ṁ (1) Moreover, the disk aspect ratio is assumed constant; namely, h = h 0 = 0.03 ).
In the disk outside of the iceline, the gas surface density is described by the power-law function: where the total disk mass is M disk = 0.04M . The disk aspect ratio is modeled by the following power-law function: This function provides the continuous disk aspect ratios between the inner and outer disks and can be derived from the disk temperature T ∝ r −1/2 in the outer disk. The gas temperature is expressed as follows: for both the inner and outer disk parts. The iceline is set at the disk location where the temperature is 170 K (Hayashi 1981). In our model, the iceline location (r ice ) is proportional to M . Because the inner disk is considered viscously relaxed, our simulations start from the viscous timescale at the iceline: Thus, our model focuses on the early history of the formation of planetary systems. The inner boundary of the disk coincides with the magnetospheric cavity radius: where B is the magnetic field strength of the stellar surface, R is the stellar radius, and ζ = (1/8) 1/7 0.74 is a dimensionless factor ranging from 0.5 to 1 for an aligned dipole (e.g., Frank et al. 2002;Chang et al. 2010). We adopt B = 180 G, and R = 0.5R for 0.08M stars, where R is the Solar radius (Reiners et al. 2009). We assume that the stellar radius increases linearly with increasing stellar mass. This dependence is motivated by the estimation in Hayashi (1966) and the observed empirical mass-radius relationship of M dwarfs (Bayless & Orosz 2006). However, recent studies demonstrated that the stellar size is more complicated and the stellar mass dependence changes in the pre-main-sequence evolution when the inefficient injection of accretion heat is considered (e.g., Kunitomo et al. 2017). Ignoring this complication, we substitute the aforementioned values and Equation (1) into Equation (6), and we obtain which is approximately proportional to the stellar mass.

Pebble Flux
In our model, there is a pebble forming front at a certain orbital radius (Section 2.1). The orbital radius of the pebble front for t = t peb is given by After the pebbles have formed at r g , they drift inward with the velocity (Adachi et al. 1976;Weidenschilling 1977) v r (τ p ) = −2η where the dimensionless pressure gradient is η = 5h 2 /4 and v Kep is the Kepler velocity in the outer disk. With Equation (3), v r becomes where v Kep,ice is the Kepler velocity at the iceline. Because v r does not depend on r in the above equation, the pebble travel time from r g to r ice can be simply expressed as follows: When t = t peb , pebbles form and begin to drift inward at r g ; their accretion rates areṀ p (r g , t peb ) = 2πr gṙg Z 0 Σ g . These pebbles reach the iceline when t = t peb + t drift . The pebble mass flux at the iceline is given bẏ where dt peb /dt is the correction due to t drift . Owing to t drift ∝ r g ∝ t 2/3 peb , the following expression holds: Because t drift /t peb ∝ r −1/2 g , t drift /t peb decreases over time. The pebble-to-gas mass flux ratio is as follows: In addition, the pebble drift time is added as a correction term for the final pebble reach time: The contribution of t drift (R disk ) is 0.58 times of that of t peb (R disk ) when R disk = 100 au and M = 0.08M , and this ratio increases as the stellar mass increases. When t = t end , we assume that the pebble flux is equal to zero at the iceline. We note that our F p/g and t end are different from those in Ormel et al. (2017) since we consider the effect of the pebble drift. This makes F p/g smaller and t end longer compared to those in Ormel et al. (2017). It is worth noting that t drift /t peb ∝ r −1/2 g , which indicates that t drift /t peb is larger than 0.58 at an early time.
The pebble flux inside the iceline is calculated in a similar way. However, τ s = 10 −3 is used for silicate grains. The Stokes number of the pebbles affects the drift timescale of a pebble going from r g to r when r < r ice ; i.e., Thus, the pebble-to-gas mass flux ratio inside the iceline is computed based on this drift time.

Pebble-to-gas density ratio
The pebble density at the midplane is important in the formation of planetesimals through the streaming instability (e.g., Youdin & Goodman 2005;Johansen et al. 2007) and the growth of planetesimals and protoplanets through pebble accretion (e.g., Lambrechts & Johansen 2012). Making use of the pebble surface density Σ p =Ṁ p /(2πrv r ) and scale height H p = α/τ p H g (Dubrulle et al. 1995;Youdin & Lithwick 2007), the pebble-to-gas density ratio at the midplane is given by where f ice represents the pileup profile around the iceline. We adopt a Gaussian function for f ice as follows: where δ ice = 0.05 is the typical width of the iceline and f peak = 19. Our choice of f peak is similar to f Σ,peak in . Although f peak changes according to parameters such as α, F p/g andṀ g (e.g., Hyodo et al. 2019), in this study, a constant f peak value is used in the simulations for simplicity. We choose f peak = 19 to produce a similar peak value of the density ratio to that in Ormel et al. (2017) (Fig. 3) in the α = 10 −3 case. Our f peak value is, however, larger than Ormel et al. (2017) since our F p/g is smaller (Equation (14)). A larger f peak value would be explained back-reaction of silicate grains (Hyodo et al. 2019). Figure 1 shows the time evolution of ρ p /ρ g at the iceline for three cases. The pebble-to-gas density ratio at the iceline decreases with time because the pebble mass flux decreases with time (Equation (14)). When the disk size is small, ρ p /ρ g is high but decreases quickly to zero. In small disks, the pebble mass flux is high because the total disk gas mass (M disk ) and metallicity (Z 0 ) are fixed and t peb is short. Around the massive stars, the pebble-to-gas mass flux ratio F p/g is low due to its inverse dependence onṀ g (Equation (14)) and a small value of t peb .

Growth of protoplanets
We assume that a protoplanet forms through the streaming instability and subsequent growth at the iceline. Given its radius is ∼ 100 km, the protoplanet is assumed to have the initial mass of 4.2 × 10 21 g. The pebbles around the iceline pileup where r > r ice (1 − δ ice ) would be a mixture of ice and silicates. We assume that only icy pebbles exist in this region. Thus, the Stokes number of a pebble is τ p . In this region, the scale height of the icy pebbles is small. The pebbles are concentrated around the disk midplane, and the pebble accretion takes place in a 2D manner (i.e., 2D mode). In the region where r < r ice (1 − δ ice ), there are silicate grains without the icy mantle, which are relatively small in size and thus are well coupled with gas. Consequently, the scale height of silicate grains is almost equal to the gas scale height. The silicate pebbles are widely distributed in the vertical direction and only those around the midplane are able to accrete onto protoplanets in a 3D manner (i.e., 3D mode). In these two pebble accretion modes, the pebble accretion efficiencies are given by the fraction of the pebble accretion rate to the pebble flux (Guillot et al. 2014;Ida et al. 2016;Ormel 2017): where q pl = M pl /M is the protoplanet mass-to-central star mass fraction. The pebble accretion efficiencies are less than 20% even for a planet with the pebble isolation mass, and this is the reason why the effect of the pebble loss on the pebble density evolution caused by accretion is not considered in this study (e.g., Lambrechts & Johansen 2014; Guillot et al. 2014). The pebble accretion rate is given byṀ The protoplanets stop growing when the pebble isolation mass is reached (Lambrechts et al. 2014):

Planetary migration
As a result of planet-disk tidal interactions, protoplanets migrate inward based on either type I or type II migration depending on the gap opening in the gaseous disk. These migration timescales are where the factor γ I depends on the disk temperature and surface density structures (Tanaka et al. 2002;Kley & Nelson 2012), and the factor 1 + 0.04K is related to the gap around a protoplanet (Kanagawa et al. 2018), When the mass of a protoplanet approximately reaches M iso , a partial gap is opened, and the migration becomes slower than the type-I migration. Applying q pl = h 3 to the above equation, the maximum value of K is about 30, and thus the type-II migration timescale, t II , of the protoplanet becomes about twice as slow as t I . It would be worth noting that the above migration timescales are derived in the absence of dust feedback, which would otherwise affect the migration timescales (Kanagawa 2019;Hsieh & Lin 2020).
By considering δ ice t I ≈ t grow = M pl /Ṁ pl , where a protoplanet crosses the iceline inner edge when its mass reaches The fraction of the crossing mass (M cross ) to the isolation mass is This fraction reflects the composition of the protoplanets: the early formed protoplanets are almost entirely composed of iceline materials, and the later protoplanets have less iceline materials. Considering ρ p /ρ g ≥ 1 at the iceline, F p/g ≥ 0.62 (Equation (17)), and M cross /M iso 0.4. The fraction of iceline materials in our estimation is higher than that in Ormel et al. (2017) due to the larger f ice . After the protoplanets cross the iceline inner edge, a protoplanet is added within the iceline at every t pl . Hence, the true time span between the appearance of protoplanets is given by the timescale on which a protoplanet leaves from the iceline region plus the appearance timescale of the next protoplanet, which includes the formation and growth of an embryo: where t cross is the timescale over which a protoplanet crosses the iceline and is given by t cross ≈ δ ice t I where q pl ≈ M cross /M . We note that the contribution of t grow is included in M cross . Because the migration timescales do not depend on r (Equations (23) and (24)), the protoplanets migrate inward exponentially over time. After the first planet reaches r c , the other planets are trapped in mean motion resonances. The resonance in which planets become trapped is determined by the comparison between the migration timescale and the critical migration timescale (Ogihara & Kobayashi 2013). Planets around low-mass stars tend to be trapped in the closest first-order resonances of the resonances that exist in the inner orbits than in the orbits of planets before the resonant trapping. This is because the migration timescale increases and the critical migration decreases at r c with decreasing M ; this allows planets to migrate more slowly and to be trapped in any first-order resonances. As the orbital resonances narrower than 6:5 could be observationally rare (Fabrycky et al. 2014), we restrict ourselves to the resonances from 2:1 to 6:5 for the resonant trapping in this study.
We consider that the innermost planets are at r c , no matter how many planets are trapped in resonances, for simplicity. The location of the innermost planet would be determined by the torques on it. The innermost planet does not only feel the inward migration (negative) torques of outer planets but also feel opposite direction (positive) torques due to the disk discontinuity (Ogihara et al. 2010;Liu et al. 2017). The innermost planet stays at r c when the positive torques are balanced with the migration torques. Furthermore, planets in resonances do not cause orbital instabilities since their eccentricities are quickly damped by the tidal interaction with the disk gas (Iwasaki et al. 2001;Tanaka & Ward 2004).

Termination of Planetesimal Formation
The formation of planetesimals is stopped when one of the following three conditions is satisfied.
2. Pebbles from the outermost disk radius reach the iceline, t = t end .

The planets are trapped in resonances around r ice
When one of the above conditions is satisfied, the simulations are stopped after all formed protoplanets stop growing and migrating. The time at which Condition 1 is satisfied is estimated based on Equations (14) and (17): where t ρ,drift is t drift at r g (t ρ,peb ). Condition 2 is represented by Equation (15), and Condition 3 is satisfied when the separation between the outermost planet and the iceline inner edge is less than the separation of the 6:5 resonance. When this condition is satisfied, it is expected that the pileup structure near the iceline is affected by the outermost planet. This condition can be approximately expressed as The maximum number of planets that satisfies this relationship is N max : Condition 3 can be recast to the time condition related to t N : where t v,ice is the time at which the initial protoplanet is added (Section 2.2) and is shorter than t ρ and t end .
Considering that N max ∼ 6, we have t N 5(t cross + t pl ) + t v,ice . When t end or t ρ is shorter than t N , the number of planets is estimated to be where t ter = min (t end , t ρ ) is the termination timescale of the planetesimal formation. The timescale of t ter − t v,ice represents the duration for all the planetesimals to be able to form progressively at the iceline. We refer to this timescale as the planetesimal forming duration. These conditions are shown in Figures 1 and 2. When M = 0.08M and R disk = 100 au, t end = 4.4 × 10 5 yr is shorter than t ρ . The stellar mass affects t end and t ρ through t peb and F p/g . Around a massive star, the planetesimal formation duration is short because both t end and t ρ are short. The dependences of t end and t ρ on the disk size have opposite outcomes: t end increases with increasing disk size because the pebble forming time (t peb ) and the drift time (t drift ) increase (Equation (15)); by contrast, t ρ decreases with increasing disk size (Equation (30)). This dependence arises from F p/g ∝ Σ g ∝ R −1 disk ; thus, the pebble-to-gas mass flux ratio decreases with increasing disk size because the gas surface density decreases in our model setting. When M = 0.08M and Dependences of tρ (Equation (30)) and t end (Equation (15) R disk = 200 au, t ρ is estimated to be 2.2 × 10 5 yr based on Equation (30) and it is 2.8 × 10 5 yr according to the calculation (Figure 1). Thus, t ρ is shorter than t end = 1.1 × 10 6 yr.

Time Evolution
In this section, we present the time evolution of protoplanets. Figure 3 shows the time evolution of the protoplanets when M = 0.08M , R disk = 100 au, and t pl = 4.12×10 4 yr. In this case, eight planets are formed. The protoplanets grow rapidly around the iceline region and reach the pebble isolation mass at ∼ 0.1 au. More specifically, the first four protoplanets reach the pebble isolation mass before they leave from the iceline region due to the high ρ p /ρ g value (Figure 1). The later formed protoplanets have smaller crossing masses because the pebble-to-gas mass flux ratio decreases over time (Equations (14) and (28)). The final protoplanet has the mass ratio M cross /M iso = 0.84. Namely, the inner four planets have the same composition as the iceline region, and the outer four planets are mainly composed of 84 % and more of the iceline materials.
After a protoplanet leaves from the iceline region and t pl is passed, the next protoplanet begins to grow at the iceline. The typical time interval between the protoplanet appearance is about 5.5 × 10 4 yr, where t pl is the main component of the time interval due to the quick growth of protoplanets (Equations (26) and (29)). In this case, the protoplanet reaches the isolation mass in the interval, i.e., before the growth of the next protoplanet begins.
The protoplanets migrate inward and become trapped in resonances after the innermost planet reaches r c . The inner five pairs are trapped in 3:2 resonances and the outer two pairs are trapped in 2:1 resonances. The outer planet pairs are trapped in larger period ratio resonances due to the longer growth timescale. Because the resonant pair of the planets is given by the closest inner resonance around low-mass stars (Section 2.4.2), the evolution of their period ratio determines in which resonance planets are trapped. The period ratios of adjacent planets become the largest value when the outer planet reaches the isolation mass because their migration timescales become identical (Equations (23) and (24)). This maximum period ratio is estimated by (t cross + t pl + t grow )/t I ; i.e., the timescale for the next formed protoplanet reaching the isolation mass is divided by the migration timescale, which is then converted into the typical orbital separation. The orbital separation given by the time interval between the protoplanet appearance (t cross + t pl ) is almost the same between planets in individual systems. The separation between the planets becomes larger as F p/g decreases, which increases the growth timescale of the next formed protoplanet (t grow , Equation (26)). Accordingly, the inner planets are trapped in less separated resonances  (35)), and the dash-dotted line represents tcross (Equation (29)). and the outer planets are locked in more separated resonances.
The formation of planetesimals is stopped when the eighth planet is formed. In this case, Condition 3 is satisfied (Section 2.5); namely, the separation between the outermost planets and the innermost radius of the iceline (r ice (1 − δ ice )) is smaller than the separation of the 6:5 resonance.

Resonant Chains of Planets
This section presents the simulations with varied t pl . Figure 4 shows the results for the case with M = 0.08M and R disk = 100 au. Evidently, the planets are trapped in more widely separated resonances (i.e., small p resonances) and the number of planets decreases with increasing t pl . Another interesting feature is that the planets in individual systems are trapped in one or two kinds of resonances. This can be interpreted by the values of (t cross + t pl + t grow )/t I : While the outer planets have longer t grow , planets have similar values in individual systems, which means that the maximum separations of planets are similar. It is worth noting that the outermost planet pairs are often trapped in a closely separated resonance than the inner ones. This is because the separation between the second outermost planet and the iceline inner edge is small. In addition, the resonant trapping is considered after the protoplanets leave from the iceline region.
In Figure 4, t cross , t ρ,pl , and t end,pl are also plotted, where t ρ,pl = t ρ − t v,ice 1.5 log 2 (r/r c ) , t end,pl = t end − t v,ice 1.5 log 2 (r/r c ) . (35) These timescales are derived from Equations (34) and (31) with p = 1. Using these timescales, the termination Conditions 1 or 2 can be written as (t pl or t cross ) (t ρ,pl or t end,pl ), and we can present these conditions in the figure. For M = 0.08M and R disk = 100 au, the curve corresponding to t end,pl agrees with the orbital distribution of planets when t pl 10 5 yr. When t pl < 10 5 yr, the planets fill orbital radii between r c and r ice . The number of planets is well represented by Equation (32). When t pl < 2 × 10 4 yr, more than 10 planets are trapped in 4:3 or closer resonances. These planets could cause orbital instability and experience giant impacts after the disk gas removal (Matsumoto et al. 2012;Matsumoto & Ogihara 2020). The resulting giant impacts continue to grow the planet mass up to the so-called ejection mass and even eject planets, leading to non-resonant planetary systems and freefloating planets due to the large total mass of planets . This suggests that observed planets in resonant chains would be mainly composed of 3:2 or 2:1 resonances. In the following sections, we adopt N crit = 10 as the criterion for the stable planet systems in resonant chains and discuss the final configuration of planets.
The abundance of the iceline materials in the formed planets is almost as high as that in the case that t pl = 4.12 × 10 4 yr. This is because planets are formed in the high F p/g environment when M = 0.08M and R disk = 100 au. In this stellar mass and disk size, the termination Condition 1, ρ p /ρ g < 1, is not satisfied (Figure 4). All planets are formed in the environment of the high pebble-to-gas density ratio at the iceline (Figure 1), i.e., high F p/g (Equation (17)). The formed planets are mainly composed of the iceline material since planets grow quickly in the high F p/g cases (Equation (28)).

Disk Size Dependence
The disk size affects the pebble-to-gas mass flux ratio (F p/g ∝ R −1 disk , Equation (14)) and the termination condition of the protoplanet formation. Figure 5 shows the number of planets and the resonance in which the innermost planet pairs are trapped. These indicate the configurations of the final planets because the planets are trapped in resonant chains composed of one or two resonances (Section 3.2). In Figure 5, we show the results for the R disk = 50 au, 100 au, and 200 au cases around 0.08M stars. Interestingly, the relationships between the innermost resonance and t pl in these results are similar. The planets tend to be trapped in more widely separated resonances (i.e., smaller p resonances) as the disk size increases. This is because F p/g decreases and t grow increases with increasing disk size.
The number of planets in the R disk = 100 au case is the largest for any t pl . This is because the termination timescale of the planetesimal formation (t ter = min (t ρ , t end )) is the longest in the R disk = 100 au case for these three cases (Figure 2). The termination timescales are 1.8 × 10 5 yr in the R disk = 50 au case, 4.4 × 10 5 yr in the R disk = 100 au case, and 2.2 × 10 5 yr in the R disk = 200 au case. While the planets fill orbital radii in the R disk = 100 au case when t pl < 10 5 yr, the numbers of planets are regulated by t end in the R disk = 50 au case and by t ρ in the R disk = 200 au case. As a result, the numbers of planets of the R disk = 50 au and R disk = 200 au cases are similar and are smaller than in the R disk = 100 au case.
The maximum disk size for the formation of planets in resonant chains in this scenario is determined based on t ter (= t ρ ) < t cross . When this condition is satisfied, only one planet is formed and the next planetesimal cannot be formed. The maximum disk size is R disk 320 au around 0.08M stars. Furthermore, no planets are formed when R disk ≥ 360 au because ρ p /ρ g is less than 1 even at t = t v,ice .
The disk size can affect the resonant chains. In the R disk = 100 au case, after the gaseous disk is dissipated, the planets can escape from resonant orbits due to orbital instabilities when they are trapped in 4:3 or more closely separated resonances because the number of planets exceeds N crit . However, in the R disk = 50 au and R disk = 200 au cases, planets trapped in 4:3 or more closely separated resonances do not always cause orbital instabilities because the number of planets is less than N crit . If any exoplanetary systems composed of Earth-mass planets trapped in 4:3 resonances are observed around 0.08M stars, they would have formed in disks with sizes of R disk ∼ 50 au or ∼ 200 au.
The compositions of the formed planets depend on the disk size since F p/g ∝ R −1 disk . This indicates that while the number of planets and their resonances are similar between the R disk = 50 au and 200 au cases, their compositions are different. In the R disk = 50 au cases, planets are composed of the iceline materials. In contrast, in the R disk = 200 au cases, the mass fractions of the iceline materials in planets are between 0.77 and 0.49.

Stellar Mass Dependence
The stellar mass dependences on the resonance of the innermost pairs and the number of planets are shown in Figure 6. As the stellar mass increases, the innermost planets become trapped in wider separation resonances (small p). The stellar mass dependence on the period ratio can be estimated using (t cross + t pl + t grow )/t I . We note that the stellar mass dependence is included in the pebble-to-gas mass flux ratio, F p/g ∝ M −7/15 in addition to the insignificant contribution from t drift /t peb . The stellar mass dependence on (t cross + t pl + t grow )/t I is positive; i.e., the period ratio increases with increasing stellar mass. Fewer planets are formed around more massive stars. This is because the planetesimal forming duration becomes shorter. Around massive stars, the termination timescale of the planetesimal formation becomes shorter ( Figure 2) and the viscous timescale at the iceline becomes longer (Equation (5)). Because the planetesimal formation takes place from t v,ice to t ter , the number of planets decreases with increasing stellar mass. In our results, no planets are formed when M ≥ 0.4M because ρ p /ρ g is less than 1 at t v,ice . This is one possible explanation of why the fraction of resonant or near resonant planetary systems around low-mass stars seems to be high.
Next, we consider the condition in which N ≤ N crit planets are formed. The protoplanet appearance time that satisfies N ≤ N crit is longer around less massive stars. This condition is satisfied for the following cases: 1.6 × 10 5 ≤ t pl /T Kep,ice ≤ 1.6 × 10 6 in the M = 0.08M case, 3.8 × 10 4 ≤ t pl /T Kep,ice ≤ 6.2 × 10 5 in the M = 0.16M case, and 8.1 × 10 2 ≤ t pl /T Kep,ice ≤ 1.0 × 10 5 in the M = 0.32M case. The protoplanet appearance time to form the stable planet systems in resonant chains varies with the stellar mass. We suggest that the tyical protoplanet appearance time can be estimated based on the fraction of planets near resonant chains as the function of the stellar mass.
The stellar mass affects the compositions of the formed planets due to the stellar mass dependence of F p/g . When we consider the planet formation around a more massive star, F p/g becomes smaller, and planets are more silicate-rich. In the case that M = 0.32M , the mass fractions of the iceline materials in planets are between 0.56 and 0.49.

DISCUSSION
In this paper, we present a time evolution model for multiple-planet formations based on the pebble-driven planet formation scenario . Although the relationship between t pl and the resonant chains of the planets can be determined, our results are based on a number of assumptions. In the following subsections, the assumptions and the formation of the TRAPPIST-1 system are discussed.

Disk Structure
In the results of this study, the number of planets depends on the termination time of the planetesimal formation (t ter ), which is the minimum of the final pebble reach time (t end ) and the time when ρ p /ρ g < 1 (t ρ ). We assume that the pebbles are formed through the coagulation of icy grains that originally exist in the outer disk.
Namely, we neglected the dust supply from the remnant star-forming core and surrounding interstellar medium, while these would be important to explain the discrepancy between the observed dust mass and observed exoplanet mass (Manara et al. 2018). These processes affect our results when the following conditions are satisfied: they provide sufficient pebble mass fluxes for keeping ρ p /ρ g > 1; the durations of these processes are longer than the termination timescale of the planetesimal formation (t ter < 10 6 yr). Such a high and long-term infall onto the star-disk system was suggested by the numerical simulation of turbulent molecular clouds (Padoan et al. 2014). While it is not clear whether this mass discrepancy is due to the dust supply or dust scattering (Ueda et al. 2020), a high and long-term infall would help to grow more planets.
We considered the evolution of protoplanets in the early phase of the disk evolution (t ter < 10 6 yr). This short timescale is the reason why the evolution of the gas disk is not taken into account. The lifetime of gas disks around low-mass stars is typically 3 Myr (Luhman 2012). The gas accretion does not change significantly before < 1 Myr. Consequently, the location of the iceline does not move since the viscous heating, which arose from the gas accretion, is almost unchanged (Oka et al. 2011). Moreover, Takahashi & Muto (2018) showed that the disk evolution can be influenced by the MHD disk wind (e.g., Suzuki & Inutsuka 2009;Bai 2014), which enables the creation of a ring-hole structures within 10 6 yr. Such disk evolution can also affect the inward migration of planets (Ogihara et al. 2015). The efficiency of the disk evolution due to the MHD disk wind around M dwarfs is a key parameter that should be considered in the final planet configuration, i.e., to determine whether the planets are in resonant chains.
We adopt f ice = 19 in all simulations. If we consider smaller f ice , our results change as follows: termination Condition 1 is quickly satisfied since t ρ,peb ∝ f 3 ice (Equation (30)); planets are formed in narrower parameter spaces of M and R disk due to short t ρ,peb (Figure 2); the number of the planets in resonances becomes smaller due to short t ter and long t grow (Equation (26)); planets are more silicate-rich (Equation (28)).

Isolation Mass
We adopt a simple formula as the pebble isolation mass (Equation (22)). According to recent hydrodynamical simulations, which provided more accurate expressions (e.g., Bitsch et al. 2018;Ataiee et al. 2018), the pebble isolation mass is about 59% of our expression (here, we use the expression of Ataiee et al. 2018). This less massive isolation mass affects our results, slightly. When the isolation mass is less massive than our expression, the growth timescale becomes shorter. Since the resonances in which planets are trapped depend on the timescale over which the next protoplanet reaches the isolation mass (Section 3.1), the relation between the resonances and t pl is slightly changed. For example, in the parameter set of Figure 3, eight planets are formed and all are trapped 3:2 resonances when we consider the isolation mass given by Ataiee et al. (2018). This also affects compositions of planets: all planets reach the isolation mass around the iceline region in this parameter set.

Filtering
We do not consider the pebble filtering effect, while we consider the formation of multiple-planet systems. In our model, protoplanets quickly grow up around the iceline region (Equations (26) and (28)). Earlier formed protoplanets obtain the isolation mass before the next protoplanets begin to grow since they satisfy the condition that the growth time to the isolation mass is shorter than t cross +t pl . This condition is not satisfied in the case of the later formed protoplanets in large R disk disks or around massive stars. These protoplanets accrete filtered pebbles. Filtering by a single protoplanet is inefficient (Equation (19), see also Lambrechts & Johansen (2014)). The efficiency of filtering would depend on the size distribution of planetesimals around the iceline (Guillot et al. 2014).

Appearance of Protoplanets
We consider that protoplanet precursors of 100 km in size successively appear, separated by a time t pl after the preceding protoplanet has crossed the iceline inner edge. In this section, we disucss t pl , although we need the hydrodynamical simulations that include planetesimals and pebbles around the iceline to consider a more realistic picture. The appearance timescale of the protoplanet would be given by the planetesimal formation timescale and its growth timescale. The planetesimal formation timescale via the streaming instability depends on the initial local pebbles-to-gas ratio, the Stokes numbers of pebbles, and the strength of turbulence (e.g., Youdin & Johansen 2007;Chen & Lin 2020;Umurhan et al. 2020). Some simulations show that the quick growth of clumps or filaments operate on timescales of 10 3 T Kep (e.g., Johansen et al. 2009;Yang et al. 2017). However, the turbulent diffusion prolongs the growth timescale of the streaming instability.
The initial growth of planetesimals that are precursors of protoplanets around the iceline is pebble accretion of scattered planetesimals and/or runaway growth of planetesimals (Liu et al. 2019b;Schoonenberg et al. 2019). These scattering and runaway growth timescales can be estimated by the dynamics of planetesimals (e.g., Kokubo & Ida 2000;Liu et al. 2019b). The timescale of the viscous stirring, which provides the timescale for planetesimals to scatter, is ∼ 10 5 T Kep . The timescale of the runaway growth is ∼ 10 4 T Kep -10 5 T Kep . These timescales agree with the protoplanet appearance time in Schoonenberg et al. (2019). These suggest that t pl would be ∼ 10 4 T Kep,ice -10 5 T Kep,ice . It would be worth noting that these timescales, especially for the runaway growth, depend on the surface density of planetesimals. The above timescales would change if the surface density of planetesimals changes significantly. The planetesimal forming efficiency is important to consider the subsequent growth of formed planetesimals.

Application to TRAPPIST-1
TRAPPIST-1 has a mass of ≈ 0.08M and hosts seven transiting Earth-mass planets with period ratios of 8:5, 5:3, 3:2, 3:2, 4:3, and 3:2 starting from the innermost pair (Gillon et al. 2017;Luger et al. 2017). It was suggested that the inner three planets are in each 3:2 resonant chain and that they experience the expansions of orbital separations induced by the stellar tide (Papaloizou et al. 2018) or by the magnetospheric rebound mechanism Liu et al. 2017). Based on these studies, we consider the planets in the initial resonant chain composed of 3:2 and 4:3 resonances. In our results, seven planets are trapped in 3:2 resonant chains in the R disk = 200 au case ( Figure 5). In this case, the particular protoplanet appearance time at which the planets are trapped in the 3:2 resonances and 4:3 resonances is about 10 5 T Kep,ice . When t pl of the TRAPPIST-1 g is slightly shorter than 10 5 T Kep,ice and the others are slightly longer, our results reproduce the planets trapped in the 3:2, 3:2, 3:2, 3:2, 4:3 and 3:2 resonant chains.
Next, we discuss the compositions of TRAPPIST-1 planets. The studies of interior modeling suggested that the water mass fractions of TRAPPIST-1 planets have the following two features: their water mass fractions are 25 wt%; these fractions are uniform or increasing with orbital periods (e.g., Dorn et al. 2018;Unterborn et al. 2018;Agol et al. 2020). Our model predicts that planets are mainly composed of iceline material. The mass fractions of the iceline material contents in the planets are between 0.77 (the innermost) and 0.50 (the outermost). While our results suggest seemingly higher in the water mass fractions, the water mass fractions of TRAPPIST-1 planets can be explained if we consider the accretion of silicate grains released by the sublimation of icy peb-bles around the iceline (Hyodo et al. 2019). It would be worth noting that Agol et al. (2020) showed that the water mass fractions of TRAPPIST-1 planets are about equal or less than 5 wt% if they have Earth-like interior structures. Such small fractions of waters are explained only when almost all iceline materials are silicate. In our results, the outer planets have less iceline materials. To explain the dependence of the water mass fractions of TRAPPIST-1 planets on orbital periods, it is needed that the iceline materials become more water-rich as icy pebbles formed at the outer region reach.

CONCLUSION
Recent observations have revealed the existence of multiple-planet systems around low-mass stars, which comprises planets with masses of ∼ 1M ⊕ planets. Interestingly, planets are near resonant orbits in most of these systems. In this study, we considered the formation of multiple-planet systems around low-mass stars. We construct a time evolution model of protoplanets around low-mass stars based on the scenario proposed by Ormel et al. (2017). The model considers the formation of planetesimals at the iceline and their growth due to pebble accretion.
We find that the protoplanet appearance timescale (t pl ) plays an important role in the configuration of resonant trapping and in the number of planets. When t pl is short, many planets are formed, and they are trapped in more closely separated resonances. These planets cause orbital instabilities after the disk gas removal, and they are expected to end up being planets in non-resonant orbits and even free-floating planets. As t pl increases, the number of planets decreases, and the planets become trapped in more widely separated resonances. This indicates that there exists a range of t pl for forming planets in resonant chains. The formed planets are trapped in the resonant chains that are composed of one or two kinds of resonances. We predict that inner planets have richer iceline materials. It is worth noting that our simulation results are based on simple analytical estimations. Nevertheless, these estimations provide a way for expanding the formation of the single-planet system to the formation of the multiple-planet system.
The disk size and the stellar mass are used as parameters in simulations. We find that the number of planets varies with the disk size because the planetesimal forming duration, which represents the duration determined by the condition that the pebble density is larger than the gas density at the iceline, changes. The disk size also affects the compositions of planets. We can predict the disk size and t pl of observed planets from their features such as the resonances, number of planets, and their compositions. Our model can well reproduce the TRAPPIST-1 system when the disk size is 200 au and t pl is about 10 5 orbital periods around the iceline. We also find that the range of t pl for forming planets in resonant chains depends on the stellar mass. Resonant planets are not formed when the stellar mass is larger than 0.4M . This is a possible explanation of why there are not many planets near resonant orbits around stars with masses of ∼ 1M .