Introduction

Compact exoplanetary systems have multiple planets that are mostly solid by mass, with approximately co-planar and circular orbits at small radial distances of only about 10 to 102 times their host star’s radius (R*)1. Compact systems also display a similar ratio of the total mass of planets in each system compared to the stellar mass (blue box, Fig. 1). While there could be undetected systems with smaller mass ratios, the largest close-in planets have likely been detected, and so the maximum mass ratio of approximately 10−4 appears significant. This is not explained solely by dynamical stability: e.g., a Trappist-1-like star could accommodate three 5 Earth-mass planets totaling about 6 × 10−4 times the stellar mass, M*, but such large mass ratio systems seem rare.

Fig. 1: Estimated total mass of transiting compact systems, Mtot, scaled to the stellar mass, M*.
figure 1

a (Mtot/M*) for compact systems having ≥3 known planets that orbit a single star within a < 0.5 au (circles, blue box). Points are ordered left-to-right by ascending stellar mass. Data are from www.exoplanetarchive.ipac.caltech.edu. For cases without mass estimates, we use the observed planet radius, increase the estimated radius uncertainty by a factor of 2, and then apply a radius vs. mass relation46. Light [medium] blue circles are systems with all [some] planetary masses estimated from this relation, while dark blue circles are systems with measured planetary masses. We include only systems with resulting errors ΔMtot/Mtot < 1. b Distribution of compact system mass ratios. Over a wide range of stellar masses (M* ≈ 0.1 to 1.3 times the Sun's mass), compact multi-planet systems display a common mass ratio, with 90% of systems having 3 × 10−5 < (Mtot/M*) < 3 × 10−4 (gray shaded region). This mass ratio is more similar to that of the gas giant satellite systems (squares, yellow box) than to the inner or outer planets in our Solar System (triangles, red box). Source data are provided as a Source Data file.

Full size image

How compact systems formed and why they are so different than our Solar System is a fundamental question. A puzzle is how planets with orbital periods of only days-to-weeks survived within a circumstellar gas disk that persisted for about 106 yr. Gravitational interactions with a gas disk typically cause a planet’s orbit to spiral inward, which would lead to rapid planet loss in a compact disk with a stellar metallicity. For some disk conditions, orbital migration is outward2, but these must not have predominated in compact systems whose planets remained close to the star. Prior works have shown that planet loss may be prevented if the gas disk has an inner cavity that halts migration of planets near its outer edge3,4,5,6,7. Such cavities are common around T Tauri (Class II, post-infall) stars8.

However, there has been no explanation for the common mass ratio among compact systems seen in Fig. 1. Prior models presume a priori a starting mass of accreting disk solids equal to that needed to yield the observed compact system planet masses. Yet what established that initial mass—and why resulting compact systems would display a preferred mass ratio despite varied stellar masses and disk evolutions—has been a mystery. If the preferred ratio were simply due to some process having selected for a common disk-to-star mass ratio, then the mass of solids available for forming compact planets would be proportional to the stellar metallicity. Instead, the occurrence of compact systems shows a weak dependence on metallicity, in contrast to the strong metallicity dependence of hot/cold Jupiters9. These issues motivate additional consideration of the processes that establish the initial conditions for compact system formation.

As a molecular cloud core collapses, a circumstellar disk forms because infalling gas and entrained grains have too much angular momentum to fall directly onto the star. The timescale over which infall occurs is constrained by surveys of young stellar objects. Estimated median lifetimes (half-lives), t1/2, for Class 0 + I objects undergoing infall range from about 0.1 to 0.5 Myr, with an additional 0.1 to 0.4 Myr spent in a subsequent “flat SED" phase10,11, which may be associated with remnant infall10,12.

A nearly universal assumption is that planet accretion begins only after infall ends. Yet, increasing observational evidence points to an earlier start to accretion. The mass in sub-cm particles in 1 to 3 Myr old disks seems too low to explain observed compact systems, suggesting that most solids have accreted into larger sizes by this time13,14. The flux density at different wavelengths in systems undergoing infall is indicative of grain growth15, and there is evidence for mm-sized pebbles in the collapsing envelopes and disks of protostars that are only about 105 yr old16,17. Further, large-scale dust rings and gaps in disks undergoing infall hint that large bodies may exist or be forming by this time18,19.

In this work we develop a concept in which compact system planets form during the end stages of infall. Planets accrete close to the star with their masses regulated by a balance between mass accumulation from the infalling solids and the inward gas-driven migration that becomes faster as the planetary mass increases. Across a wide range of disk conditions, we find that this balance naturally leads to total planetary system masses of a few  × 10−5 to 10−4 times the stellar mass, in agreement with the observed properties of compact systems (Fig. 1).

Results

We adopt a standard model in which gas and tiny grains infall to a circumstellar disk interior to the centrifugal radius, rc (Supplemental Figure 1), with the gas then spreading viscously20. The detailed properties of the collapsing core would determine rc20,21 (see Infall description section in Methods). Observed disks undergoing infall have mean radii of about a few  × 10 au, with wide scatter down to observational limits at approximately 10 au22. Because some disk material viscously expands beyond rc during infall, observed disk radii would be larger than rc. We here consider that compact systems result from disks with small rc ≤ 1 au, much less than the 10 au often assumed for our Solar System.

Once in stellar orbit, grains delivered to the disk by the infall undergo mutual collisions, particle growth, and gas-driven radial drift23. How grains grow into km-sized planetesimals and larger is a complex question24. A premise of our model is that planetesimals form somewhere within the region of infall (i.e., at or interior to rc) during the final infall stage.

Planetesimal formation in inner disks is plausible given recent results (see also Dust-to-planetesimal accretion section in Methods). Experiments find a 2 to 3 orders-of-magnitude increase in rocky grain adhesive forces at high inner disk temperatures25, and observations indicate particle growth in young disks interior to the ice line26,27,28, suggesting that cm-scale and larger pebbles may form interior to rc. Planetesimal formation via streaming instability requires pebble concentration29,30, and the inner disk is a favorable region for this, due to inward radial drift of outer disk particles23,31,32, particle trapping within vortices that form near rc during infall33, and/or particle buildup in the outward gas velocity region near rc34 or within a pressure maximum near the silicate evaporation radius35.

Accretion within infall-supplied disks has been studied in the context of gas planet satellite origin, including models that predict a preference for satellite systems with 10−4 times the planet’s mass36,37 (Fig. 1, yellow box). A key difference between satellite systems and compact exoplanetary systems is that for the latter, the gas disk persists after the infall ends, which would tend to cause greater migration-driven planet loss. The approximate 10−4 maximum mass ratio seen in compact systems occurs for stellar masses that vary by an order-of-magnitude, suggesting that it does not reflect a maximum absolute planet mass, Mp, but instead a maximum (Mp/M*) value, as would be expected for systems regulated by gas-driven migration whose rate is proportional to (Mp/M*), per below.

Planet accretion vs. migration

We first estimate conditions that allow planets formed during late infall to survive until the gas disk dissipates, and the expected properties of surviving systems. We consider an infall flux that is radially uniform for r ≤ rc, has a gas-to-solids ratio f, and decreases exponentially with time, \({F}_{{{\rm{in}}}}={F}_{0}\exp (-t/{\tau }_{{{\rm{in}}}})\) (where F0 is the flux at the beginning of the final infall stage, and τin is related to an observationally determined half-life as τin = 1.4t1/2). We approximate the gas disk surface density as \({\sigma }_{g}(r,t)={\sigma }_{g}(r,0)\exp (-t/{\tau }_{{g}})\), where τg is a gas disk dispersal timescale and σg(r, 0) reflects an initial quasi-steady state38 between infall and spreading due to a viscosity ν = αcH, where α < 1 is a viscosity parameter, and c and H are the gas sound speed and vertical scale height, respectively. We define β ≡ (τg / τin) > 1, so that the gas disk persists after infall ends.

Planetesimals form interior to rc and subsequently grow into larger bodies via mutual collisions and sweep-up of small solids. The growth timescale of a planet, τacc, is regulated by the solid infall rate, with τacc ~(f/ε)Mp/(2πrΔrFin) (see Supplemental Discussion 2), where Δr is the planet’s feeding zone width and ε is the fraction of infalling solids incorporated into large bodies (either by directly participating in planetesimal formation or by being accreted by large bodies).

Growing planets undergo inward Type-I migration on a timescale39

$${\tau }_{{{\rm{I}}}} \sim \frac{1}{{C}_{a}\Omega (r)}\left(\frac{{M}_{*}}{{M}_{{{\rm{p}}}}}\right)\left(\frac{{M}_{*}}{{r}^{2}{\sigma }_{g}(r,t)}\right){\left(\frac{H}{r}\right)}^{2}$$
(1)

where Ca is a torque constant, \(\Omega (r)=\sqrt{G{M}_{*}/{r}^{3}}\) is orbital frequency, and \((H/r)\propto {(r/{R}_{*})}^{{\gamma }_{c}}\) with γc ≪ 1.

A planet accretes largely in-place until it approaches a critical mass, Mcrit, for which τacc =τI36. The planet then migrates inward before becoming substantially more massive. Setting τacc = τI and solving for Mp = Mcrit gives

$$\frac{M_{\rm crit}(t)}{M_*}\approx \ 4.1 \left(\frac{M_*}{\pi r_c^2F_0} \frac{f}{\varepsilon}\Omega(r_c)\right)^{1/9}\left(\frac{\pi}{C_a}\right)^{5/9} \left(\frac{H}{r}\right)^{26/9}\left(\frac{r}{r_c}\right)^{17/18} \left(\frac{\alpha\varepsilon}{f}\right)^{2/3}\\ \, \, e^{ -t/\tau_{\rm in}\left( \frac{5}{9}-\frac{2}{3\beta}\right)}\\ \approx \ 2.8\times10^{-5}\chi\left(\frac{1}{C_a}\right)^{5/9}\left(\frac{H/r}{0.05}\right)^{26/9} \left(\frac{r}{r_c}\right)^{17/18} \left(\frac{\alpha\varepsilon/f}{5\times10^{-5}}\right)^{2/3} \\ \, \, e^{-t/\tau_{\rm in}\left( \frac{5}{9}-\frac{2}{3\beta}\right)} $$
(2)

where \(\chi \equiv {\left(\frac{{M}_{*}/(\pi {r}_{c}^{2}{F}_{0})}{16 \;{{\rm{Myr}}}}\frac{f/\varepsilon }{100}\frac{100{{\rm{days}}}}{2\pi /\Omega ({r}_{c})}\right)}^{1/9}\) is a term of order unity that depends only weakly on disk and inflow parameters.

As a mass Mcrit planet migrates inward, a new mass Mcrit planet can form in its place on a comparable timescale due to ongoing solid infall. This regulates the total mass in planets, Mtot, to (Supplemental discussion 2)

$$\frac{M_{\rm tot}(t)}{M_*}= \int^{r_c}_{r_{\rm in}}\frac{M_{\rm crit}/{M_{*}}}{\Delta r} \rm{d}r \\ \approx \ 2.2\times10^{-4}\frac{1}{\chi}\left(\frac{1}{C_a}\right)^{4/9}\left(\frac{\alpha\varepsilon/f}{5\times10^{-5}}\right)^{1/3}\left(\frac{H/r}{0.05}\right)^{10/9}e^{-t/\tau_{\rm in}\left(\frac{4}{9}-\frac{1}{3\beta}\right)} $$
(3)

where rin is an inner loss boundary (we set rin/rc =0.1), and we neglect the weak dependence of (H/r) on r given that γc ≪ 1.

The predicted system mass ratio in eqn. (3) depends very weakly on the infall rate through χ. Of importance is the (αε/f)1/3 term. The fraction of the total infalling mass incorporated into solid planets is (ε/f), and if this fraction is increased for a fixed α, accretion becomes faster and planets grow more massive before migrating inward. For a fixed (ε/f), a higher α leads to faster disk-spreading and a lower gas surface density, which slows migration and allows for more massive planets. In addition (and differently from satellite origin models,36), there are time-dependent terms in eqns. (2) and (3) due to the mismatch between the infall timescale, τin, and the gas disk dispersal timescale, τg.

The balance between accretion and migration described by eqns. (2) and (3) continues until the migration time for a mass Mcrit planet becomes much longer than the gas disk dispersal timescale, i.e., until τI(Mcrit) ≫ τg. As this transition occurs, migration subsides and mass Mcrit planets can survive. We determine the time at which a mass Mcrit planet at rc has τI = 10τg, and evaluate eqn. (3) at this time to estimate the final planetary system mass ratio that survives against Type-I migration (Supplemental discussion 2):

$$\frac{{M}_{{{\rm{tot,final}}}}}{{M}_{*}}\approx {\left(\frac{{M}_{{{\rm{tot,0}}}}}{{M}_{*}}\right)}^{1+\Lambda }{\left(\frac{1}{2}\frac{{M}_{*}}{{M}_{{{\rm{in}}}}}\frac{f/\varepsilon }{10\beta }\right)}^{\Lambda }$$
(4)

where (Mtot,0/M*) is a maximum value from eqn. (3) for t ≪ τin, Λ ≡ (4β − 3)/(5β + 3) with 0.125 ≤ Λ ≤ 0.8 for 1 ≤ β < ∞, and \({M}_{{{\rm{in}}}}\equiv \pi {r}_{c}^{2}{F}_{0}{\tau }_{{{\rm{in}}}}\) is the total mass delivered during the final infall stage. As β increases, the surviving system mass ratio decreases due to the longer persistence of the gas disk.

Numerical simulations

We simulate this process using an N-body code40 that includes a time-dependent infall of gas and solids, and gas disk torques (see N-body planet accretion simulation section in Methods). There is strong observational evidence for inner magnetospheric cavities in the disks of Class II (post-infall) stars8,41, and such a cavity can halt migration of a planet near its edge (see Supplemental discussion 4). Whether disks of Class I stars undergoing infall have inner cavities is less clear, as they may instead undergo boundary layer accretion without a cavity42. Given this, we consider two limiting cases: a disk with an inner migration-stopping cavity, and a no-cavity disk with full inward migration.

Figure 2 shows results of simulations with full migration that consider different β for otherwise equal conditions. The system mass initially increases due to the solid infall until planets of mass about Mcrit form and migrate inward. As planets are lost by migration, new planets accrete in their place as the infall and gas disk properties evolve. The total planetary system mass (Fig. 2a, solid lines) oscillates around a value that decreases with time, as predicted by eqn. (3) (Fig. 2a, dotted lines), until the transition time at which τI(Mcrit) ≫ τg, after which the total mass stabilizes. As the relative disk dispersal timescale compared to the infall timescale (β) increases (varied colors in Fig. 2a), the final system mass ratio decreases, consistent with eqn. (4) (Fig. 2b).

Fig. 2: Total planetary system mass during infall and subsequent gas disk dissipation, in units of the stellar mass.
figure 2

a System mass ratio vs. time normalized by the gas disk dispersal timescale, τg. Simulations here assume (αε/f) = 5 × 10−5 (where α is the disk viscosity parameter, ε is the fraction of infalling solids incorporated into planets, and f is the infall gas-to-solids ratio), a common infall prescription, no inner cavity (i.e., inward migration of all planets), and varied β ≡ τg/τin values indicated by color (where τin is the infall decay timescale). Solid lines are simulation results; dotted curves show predictions from eqns. (3) and (4). b Final system mass ratio vs. β from 10 simulations compared to analytical prediction (eqn. (4), red dashed line, assuming f/ε = 102). Blue shaded region shows the range appropriate for compact systems with small rc, although planetary systems survive even for larger β values. Source data are provided as a Source Data file.

Full size image

To constrain the appropriate values for β for compact systems, we model the radial evolution of a gas disk subject to a time-dependent infall, viscous spreading, and photoevaporative loss. Disks with smaller rc have infalling material deposited at smaller radii where viscous timescales are shorter and photoevaporative losses can be greater, leading to shorter disk lifetimes20,43. We find that for rc < 1 au, α > few  × 10−4, and t1/2 > 2 × 105 yr (i.e., τin ≥ 3 × 105 y), the gas disk’s evolution may be approximated as having 1.3 ≤ β ≤ 2 (see Relative timescales of gas disk dispersal vs. infall section in Methods). For this range of β, eqn. (4) implies that surviving systems formed during infall would have about 20 − 40% of Mtot,0 (assuming f/ε = 102 and Min/M* = 0.03 as in Fig. 2-b). The higher β values associated with compact disks (1.3 ≤ β ≤ 2) compared with circumplanetary disks (thought38 to have β ≈1) would lead to a somewhat lower mass ratio for compact systems compared to satellite systems, which appears the case (Fig. 1). We note, however, that there is not a sharp cutoff as β increases, and compact systems survive even with β > 2 (Fig. 2-b).

Figure 3 shows results of simulations with varied (α/εf) and 1.3 ≤ β ≤ 2. For both disks with and without inner cavities, the final system mass ratio is between 10−5 and 2 × 10−4 across three orders-of-magnitude variation in (αε/f). For otherwise similar conditions, simulations with an inner cavity (triangles) yield systems with somewhat larger mass ratios than no-cavity cases (circles) and the no-cavity analytical estimate (eqn. (4), dashed lines). When a planet nears the cavity edge, its positive co-rotation torque may become comparable to its negative Type-I torque so that its gas disk interactions no longer cause it to migrate44. The inner planet may then hold back the migration of exterior planets via capture into mutual resonances5, allowing the resulting resonant chain of planets to grow somewhat more massive than in the no-cavity case. However, the barrier to planet infall is not absolute. As exterior planets continue to grow due to the infall, the inner planet’s co-rotation torque becomes insufficient to offset the outer planets’ negative Type-I torque, and the chain of planets collapses inward. As a result, the planet accretion-migration cycle, and its dependencies on disk and infall parameters, are broadly similar in the no-cavity and cavity cases, although a cavity allows a given mass ratio system to be achieved with somewhat less restrictive parameters, i.e., with lower (αε/f) and/or higher β (see Supplemental discussion 4).

Fig. 3: Results of planet accretion simulations with varied disk and infall properties.
figure 3

Final planetary system mass scaled to the stellar mass (with M*=1M⊙) as a function of (αε/f) (α is the viscosity parameter, ε is the fraction of infalling solids incorporated into planets, and f is the infall gas-to-solids ratio). The infall rate decays with timescale τin = 5 × 105 yr, while the gas disk disperses over a longer timescale, τg = 1.3 to 2τin (colors, legend). The simulations assume either an inner disk cavity with radius rcav = 20R* = 0.13 au (triangles) or no cavity (circles). Grey region shows range for 90% of observed compact systems shown in Fig. 1. Dashed lines show eqn. (4) predictions for the no-cavity case with β = 1.3 (light blue) and β = 2 (dark blue), with H/r = 0.05 and Min/M* = 0.03. Horizontal bars show plausible viscosity ranges based on observed stellar accretion rates77, models of magnetorotational instability (MRI)78 and infall-driven shocks79, assuming f/ε = 102. Source data are provided as a Source Data file.

Full size image

The appropriate values for ε and α are highly uncertain, and could vary across different systems. However, per eqn. (4), the final system mass is not strongly dependent on these parameters: for 1 < β ≤ 2, (Mtot,final/M*) ∝ α3/8 to α1/2, and (Mtot,final/M*) ∝ ε1/13 to ε1/4. The infall metallicity (1/f) would nominally be comparable to the stellar metallicity. All other factors being equal, a higher infall metallicity (smaller f) leads to more massive planets and a larger system mass ratio, although here too the dependence is weak, with (Mtot,final/M*) ∝ f5/9Λ−4/9, so that, e.g., a factor of 2 metallicity enhancement yields only a factor of 1.2 to 1.3 increase in mass ratio for 1.3 ≤ β ≤ 2. This appears consistent with the nearly equal occurrence of warm super-Earths around stars of wide-ranging metallicities9 that has been unexplained.

Figure 4 a shows resulting planet distributions at 5 to 10 Myr from a subset of simulations with full migration. Observed compact systems have similar planet sizes within each system, a “peas-in-a-pod” architecture. In our model, this similarity emerges naturally through the regulation of planet mass by Type-I migration (eqn. (2)). Size diversity can be quantified by \({\sigma }_{{{\mathcal{R}}}}^{2}={{\rm{Var}}}\left[{\log }_{10}\left({R}_{p}/{R}_{\oplus }\right)\right]\)1, where (Rp/R⊕) is planetary radius scaled to Earth’s radius, or with a mass partitioning coefficient45, \({{\mathcal{Q}}}=N/(N-1)\mathop{\sum }_{i=1}^{N}{({M}_{\rm {p}}/{M}_{{{\rm{tot}}}}-1/N)}^{2}\), where N is the number of planets in the system. Known multi-planet compact systems have a median \({\sigma }_{{{\mathcal{R}}}}=0.14\)1, with most systems having \({\sigma }_{{{\mathcal{R}}}}=0.08\) to 0.19, and \({{\mathcal{Q}}}\, < \,0.2\). Our simulated systems (both with and without an inner cavity; see Supplemental discussion 4 and 5) compare favorably to these, having \({{\mathcal{Q}}}\, < \,0.2\) (Fig. 4a) and a median value of \({\sigma }_{{{\mathcal{R}}}}\approx 0.16\) (assuming a mass-radius scaling law46). As is also the case for prior (post-infall) compact system formation models5,6, our systems soon after gas disk dispersal have, on average, a greater number of planets (typically 5 to 7 with Mp > 10−6M*) than observed systems (Fig. 4a).

Fig. 4: Predicted system architectures and period distributions.
figure 4

a Example simulation results for a 1M⊙ star, full inward migration, and β = 1.3 (see Supplemental discussion 4 and 5 for all results including inner cavity cases). Systems are arranged by increasing (αε/f) (α is the viscosity parameter, ε is the fraction of infalling solids incorporated into planets, and f is the infall gas-to-solids ratio). Symbol color and size scale with planet mass ratio, (Mp/M*); grey indicates planets estimated to be undetectable (Supplemental discussion 1); and planets orbiting within 3 mutual Hill radii have been merged. The outer distance at which planetesimals form (assumed here to be equal to rc, blue vertical line) establishes the radial scale of the planet system, although growing planets can be scattered to somewhat larger distances too. The mass partitioning parameter, \({{\mathcal{Q}}}\), for each system is similar to those of observed systems that have \({{\mathcal{Q}}}\, < \,0.2\). b Period ratio distribution of adjacent planets in simulated systems versus observed systems from Fig. 1, including locations of mean motion resonances (vertical dashed lines). The dark [light] blue curve shows the period distribution soon after gas disk dispersal for simulations with an inner cavity [without a cavity]. The red dashed curve shows this distribution for systems that underwent a dynamical instability during longer, 100 Myr orbital integrations. As in prior works5,6, we find that later instability leads to improved agreement between model results (red dashed curve) and observations (dark and light grey curves). Source data are provided as a Source Data file.

Full size image

Another key comparison is the planet period distribution. Planets in observed compact systems are not predominantly in mean motion resonances (Fig. 4b, dark and light grey lines)47. In contrast, planet formation models involving migration (including ours and prior, post-infall models5,6) predict more resonant pairs and generally more compact systems near the time of gas disk dispersal. Our cases without an inner cavity (Fig. 4b, light blue curve) show fewer resonant pairs than cases with a cavity (Fig. 4b, dark blue curve), because the cavity favors the formation of resonant chains, but both distributions are initially distinct from observations (Fig. 4b, grey curves).

A strong case has been made that on longer timescales, dynamical instabilities cause the period ratio distribution to evolve and become more consistent with the observed distribution5,6. To illustrate such an evolution, we performed 100 Myr integrations of our systems assuming the gas disk had completely dispersed. As in prior works5,6, we find that systems that undergo instability have a period ratio distribution that more closely matches that of observed systems (Fig. 4b, red curve vs. grey curves). The system mass ratios are not substantially affected by this later evolution. Often instability involves planet merger (which leaves the system mass ratio unchanged), and in cases where planets are ejected, their mass is typically a small fraction of the total system mass. Of the systems that underwent instability in our 100 Myr integrations (shown in the (Fig. 4b, red dashed curve), 90% had their mass ratio change by  <5%. As argued previously, continued interactions on longer timescales and/or other effects48 may drive continued evolution of the period distribution and further decrease the number of observable planets5,6.

Discussion

It has been proposed7 that compact system planets accrete within a narrow ring of planetesimals at the silicate condensation distance and then migrate inward until they are stopped by an inner cavity49. The ref. 7. model is similar to ours in that inward Type-I migration produces a uniformity in planet mass within each system. A primary difference is that ref. 7 adopts a smaller rc, so that the silicate condensation line and planetesimal formation occur exterior to rc where the disk is not directly supplied by infall. The final planet system mass is then set by the assumed mass of the outer planetesimal ring. Why this initial ring mass would have the value needed to produce a preferred system mass ratio across varied stellar masses is not addressed. In contrast, we consider that planetesimals and planets form interior to rc in the region of ongoing infall. It is this condition that leads to a common system mass ratio (eqn. (4)).

Predictions of our model should be increasingly testable by observations. Accretion during infall implies a common compact system mass ratio independent of stellar mass, and that systems with somewhat lower estimated mass ratios might host yet undetected planets. An unusually weak dependence of planet and planet system mass on stellar metallicity is predicted, with (Mtot/M*) proportional to metallicity to the 0.2–0.3 power, in contrast to other models that generally imply a stronger dependence. While we here considered planet accretion interior to rc, some small solids would be carried beyond rc by expanding gas and would be available for accretion in that region. Results here imply a difference between compact system planets formed interior to rc – that display a characteristic planet mass and system mass ratio—and those that might form beyond rc, whose masses would not be regulated in the same manner. Such a transition in planet mass ratio might be detectable.

While the masses and orbital properties of compact systems formed during infall prove to be insensitive to disk conditions (e.g., the specific value of rc, infall rate, and viscosity), resulting planet compositions will be more sensitive to such parameters. Our model argues for rocky planet formation interior to the ice line7,50. As infall slows, the disk cools38 and condensation fronts shift inward. While early-formed planets would tend to be refractory rich (and more likely to be lost), some conditions (e.g., small α, larger rc) might allow for ice stability during late infall. Coupled composition and dynamical models will be needed to assess this.

Overall, whether planets that form during infall survive, and what the characteristics of surviving planets would be, will depend on rc (the radial extent of disk infall) and β (the ratio of the gas disk lifetime to the infall timescale). Here we focus on small rc and small β cases, consistent with compact systems. However, rc depends sensitively on the angular momentum of the parent cloud core and interactions of infalling material with the disk and stellar magnetosphere, so that even stars with similar masses may have substantially different rc values. For larger rc, long accretion timescales in the outer infall region may promote the formation of giant planets after infall has ended. Larger rc systems will also have larger β, due to their longer gas disk lifetimes43, implying that surviving inner planets that accreted during infall would have a lower total mass compared to the star; e.g., per Fig. 2, a mass ratio of inner planets similar to our terrestrial planets might result for β ≈8. Most broadly, our results suggest that the long-standing assumption that planet accretion commences only after infall has ended may not be valid for all systems, and consideration of this early accretionary phase is warranted.

Methods

Relative timescales of gas disk dispersal vs. infall

We show that the properties of the final planetary system formed during infall depend on the ratio of the gas disk dispersal timescale τg, to the timescale over which disk infall ends, τin (eqn. (4) and Fig. 2). To estimate the range of β ≡ τg/τin appropriate for accretion disks with small rc < 1 au, we model the evolution of a gas disk subject to a time-dependent infall, viscous spreading, and photoevaporative loss51 using a semi-analytical code that tracks the gas surface density, σg(r, t) as the last 10% of the stellar mass infalls to the disk and continues through the gas disk dispersal phase.

We consider a uniform infall flux per disk area, Fin, that decays exponentially with timescale τin, with 105 < τin < 106 yr. We consider an α-viscosity model with ν = αcH = αc2/Ω (where c and H are the gas sound speed and vertical scale height, and Ω is the orbital frequency) and few  × 10−4 < α < 0.1. We adopt typical X-ray luminosity values (LX = 1029erg s−1 for a mass 0.1M⊙ star and LX = 1030erg s−1 for a 1M⊙ mass star52).

To set the starting gas surface density and radial temperature profile, we assume that the disk is initially in a quasi-steady-state between infall and viscous spreading, where the latter depends on temperature in the α-model38. At each disk radius r, we solve for the mid-plane disk temperature, T, that balances: 1) radiative cooling from the disk’s top and bottom surfaces; 2) heating due to the star’s luminosity with \({\dot{E}}_{*}=9/4{\sigma }_{{{\rm{SB}}}}{T}_{*}^{4}{({R}_{*}/r)}^{2}(c/r\Omega )\) (where σSB is the Stephan-Boltzmann constant and T* is the stellar temperature); 3) heating due to viscous dissipation, \({\dot{E}}_{\nu }=(9/4){\sigma }_{g}\nu {\Omega }^{2}\); and 4) dissipation due to the energy difference between infalling gas and that of a local keplerian orbit, \({\dot{E}}_{{{\rm{in}}}} \sim {F}_{{{\rm{in}}}}G{M}_{*}/r\)53:

$${T}^{4}=\left(1+\frac{3{\kappa }_{R}{\sigma }_{g}}{8}+\frac{1}{2{\kappa }_{P}{\sigma }_{g}}\right)\frac{{\dot{E}}_{\nu }}{2{\sigma }_{{{\rm{SB}}}}}+\left(1+\frac{1}{2{\kappa }_{P}{\sigma }_{g}}\right)\frac{{\dot{E}}_{*}+{\dot{E}}_{{{\rm{in}}}}}{2{\sigma }_{{{\rm{SB}}}}}+{T}_{\rm {cd}}^{4}$$
(5)

where κR [κP] is the disk’s frequency-averaged Rosseland [Plank] opacity53 and Tcd = 15 K is the molecular cloud temperature. This expression assumes that heating due to stellar luminosity and infalling material is deposited at the disk surface, while viscous dissipation occurs primarily in the mid-plane53. Given T(r, 0), we calculate the initial viscosity and surface density (ν(r, 0) and σg(r, 0)) from eqn. (18) of ref. 38.

To calculate σg(r, t), we use standard methods for a viscous disk54 subject to photoevaporative loss51,52, assuming vertical thermal balance in the disk at each radius per eqn. (5). Figure 5 shows an example evolution.

Fig. 5: Gas disk surface density evolution.
figure 5

Example evolution of the gas disk surface density around a 0.1M⊙ star at different times with α = 10−3, rc = 0.06 au, and τin = 0.5 Myr. By 4 Myr, a large gap has formed due to photoevaporative losses. Source data are provided as a Source Data file.

Full size image

For low viscosity cases with α = 10−4, the initial disk mass estimated in this manner is  >0.1M*. There is then an inconsistency, because such a massive disk would have a large, gravitationally-driven viscosity with α ~ O(10−1) until its mass had decreased substantially55. For such cases (indicated by the triangle markers in Fig. 6), we set α = 0.1 when calculating T(r, 0) and σg(r, 0), and then proceed with α = 10−4 to calculate the post-gravitational instability evolution.

Fig. 6: Gas disk lifetime compared to infall timescale for compact systems.
figure 6

Predicted ratio (β) between the gas disk lifetime, τg and the infall timescale, τin, as a function of the disk viscosity parameter (α) and τin for compact systems for a (a) solar and (b) sub-solar mass star. Resulting values of β ≡ τg/τin are shown by color and numerical labels. Circles [triangles] show cases with initial gas disk masses  < 0.1M* [ > 0.1M*]. Simulations depicted by triangles assumed an initial viscosity parameter of α0 = 0.1, replicating the expected high-viscosity for gravitationally unstable, high-mass disks. Simulations with very low viscosities and fast infall, depicted by the gray-shaded regions in the lower-left corners of the plots, exhibit a disk mass time-evolution that is not well approximated by a single exponential. For cases outside this region, the evolution of the inner gas disk mass can be reasonably well approximated as proportional to \({e}^{-(t/{\tau }_{{g}})}\) (see Fig. 7). Source data are provided as a Source Data file.

Full size image

For each disk evolution with a given α and τin, we calculate an effective τg as the time needed for the gas mass interior to rc to decrease by \(1/\exp\). Note that because τg is an e-folding time, not the disk lifetime, significant gas is still present at t = τg. In general, τg decreases [increases] for larger [smaller] viscosity α-values.

We use this model to estimate β ≡ τg/τin as a function of α and τin. Across a broad parameter range, specifically for τin ≥ 3 × 105 yr and α > few  × 10−4, we find that β is between 1.3 and 2 for compact systems with rc  < 1 au (Fig. 6a, b, non-shaded regions). Figure 7 shows the mass of gas interior to rc vs. time for example cases in this parameter range. The gas decrease is reasonably well approximated by an exponential decay with timescale τg ≈1.3 τin to 2 τin. While β is small for these compact disks, resulting gas disk lifetimes (defined as when the optical depth in the near-infrared region, T = 300 K, is smaller than unity) are  > 3 Myr (for τin ≥ 3 × 105 yr and α > few × 10−4), in agreement with previously studies43.

Fig. 7: Gas disk mass evolution over time for different infall timescales.
figure 7

The gas disk mass interior to rc scaled to the stellar mass as a function of time normalized by the infall timescale, for (a) τin = 3 × 105 yr; (b) τin = 5 × 105 yr; and (c) τin = 10 × 105 yr. The colored curved lines correspond to different viscosity values, per the legend, while the shaded area surrounding each curve shows the inner disk mass evolution predicted by an exponential decay with β between 1 and 2 for the α = 0.001 cases, and with β between 1.1 and 1.5 for the other viscosity values. For disks produced by infall with small rc  < 1 au, the effective β values are generally quite low. Source data are provided as a Source Data file.

Full size image

For cases with low viscosities and fast infall rates, the gas disk evolution is not well approximated as a simple exponential decay; this regime is indicated by the gray-shaded region in Fig. 6. For these cases, the late disk mass decreases more slowly than implied by a single exponential, and the resulting disk masses are large, comparable to the mass of the star. More sophisticated models are required to accurately estimate disk evolution in this regime55.

N-body planet accretion simulation

To simulate planet accretion during and after infall, we use a modified version of SyMBA40, a symplectic integrator that allows for close encounters and collisions (treated as inelastic mergers). The background gas disk is treated analytically and its effects are added as small perturbations to the planetary orbits36. We continue the simulations for 5 to 10 times the gas disk dispersal timescale, by which time migration-driven effects have ended. Collisions among planets may occur on longer timescales, which would be expected to alter somewhat the final architecture of systems shown in Fig. 4a. As in other N-body models of compact system origin3,5,7,56, we do not include gas accretion by the growing planets, given that the planets in these systems are mostly solid by mass.

N-Body Gas Disk Description

Our accretion simulations adopt a simplified gas disk treatment based on results of the 1-D disk model described above. We assume σg(r, t) varies as a power-law in orbital radius with a constant power-law exponent (γg), and that the gas dissipates on an exponential timescale, τg, so that \({\sigma }_{g}(r,t)={\sigma }_{g}(0,{R}_{*}){({R}_{*}/r)}^{{\gamma }_{g}}\exp (-t/{\tau }_{g})\). We set γg = 0.65 and σg(0, R*) ∝ Fin/ν to approximate an initial quasi-steady state with a disk opacity of order unity38. The disk vertical aspect ratio is set as \((H/r)={(H/r)}_{{R}_{*}}{(r/{R}_{*})}^{{\gamma }_{c}}\) with γc = 0.04 (implying a temperature profile \(T(r)\propto {r}^{2{\gamma }_{c}-1}\)), and we make the simplifying assumption that (H/r) is constant with time. For simulations of accretion around a solar mass star in the main text, we set γc = 0.04, and \({(H/R)}_{{R}_{*}}=0.04\); additional simulations for a more compact, Trappist-like system with M* = 0.1M⊙ set \({(H/R)}_{{R}_{*}}=0.06\) (see Supplemental discussion 5).

Interactions with the gas disk cause orbiting bodies to undergo inward Type-I migration (with timescale from eqn. (1) and Ca ≈O(1)), density wave damping of eccentricity and inclination, and aerodynamic gas drag, as in previous simulations36. Because the gap-opening mass57 exceeds the maximum planet mass (eqn. (2)), planets are not expected to open gaps and transition to Type-II migration.

A subset of our simulations consider an inner cavity in the gas disk of radius rcav = 20R* (in Fig. 3) or rcav = 15R* (in Supplemental discussion 4), such as may result from magnetospheric truncation8 or perhaps from a larger inner disk viscosity35. We do not calculate the detailed changes in the Lindblad and co-rotation torques that occur near the cavity outer edge; these are strongly affected by edge steepness and require modeling of the cavity-producing process58,59. We assume that the cavity properties are such that the co-rotation and Lindblad torques cancel each other once a planet is interior to rcav and deactivate its Type-I migration, similar to previous studies7,60. The co-rotation torque responsible for the planet trap58 arises from interaction with gas in the local region surrounding the planet’s orbit. The dominant gas interactions responsible for eccentricity and inclination damping arise in this same region61. We thus assume that planets within the cavity still undergo damping of eccentricity and inclination, as in prior work7. Supplementary Movie 1 and 2 show two example evolutions with and without an inner disk cavity.

Our N—body simulations do not account for a gap in the gas disk that can form due to photoevaporative loss (see Fig. 5), which could decrease the migration rate relative to that assumed in the code. Because this gap tends to form late when the total planetary system mass has neared its final value, its presence would not be expected to substantially change the mass ratio values predicted here.

Infall description

In a classic collapse model, the centrifugal radius is a sensitive function of stellar mass, cloud rotation rate (ωcd) and cloud temperature (Tcd), with \({r}_{c}\propto {M}_{*}^{3}{\omega }_{{{\rm{cd}}}}^{2}{T}_{{{\rm{cd}}}}^{4}\)20. Observations62 find 7 ≤ Tcd ≤ 40 K and few × 10−15 ≤ ωcd ≤ 10−13s−1, implying a wide range of rc values, including those comparable to the radial scale of compact systems. For example, for a 1M⊙ [0.08M⊙] star, Tcd = 15 K, and ωc = 5 × 10−15s−1 [ωcd = 5 × 10−14s−1], rc =1 au [rc = 0.1 au] is predicted. Reference 7 argues for rc ≈ 10−1 au due to strong magnetic breaking. N-body simulations in the main text adopt rc = 0.6 au for a solar mass star, while those in Supplemental discussion 5 use rc = 0.07 au for an 0.1M⊙ star, with these values motivated by the Trappist-1 and Kepler-11 systems. Similar results to eqns. (2) through (4) are expected for varied rc so long as rc lies beyond the condensation line for rocky materials and planetesimals form interior to rc.

Infall is a complex process, and we utilize a simplified, tractable model for the concept development work here. We assume that during infall of the last few percent of the stellar mass to the disk, rc is constant, the infall flux Fin is uniform with radius interior to rc, and that the infall rate decreases with time as \({F}_{{{\rm{in}}}}(r,t)={F}_{0}\exp (-t/{\tau }_{{{\rm{in}}}})\). The total mass delivered during the simulated infall is \({M}_{{{\rm{in}}}}=\pi {r}_{c}^{2}{F}_{0}{\tau }_{{{\rm{in}}}}\), and our accretion simulations consider Min = 0.03 to 0.05M* and τin = 5 × 105 yr. We find that results depend primarily on the ratio τg/τin (β), rather than on the specific value of τin, so long as τin ≥ few × 105 yr (see above).

The infall of solid material to the disk interior to rc is mimicked by the addition of new bodies of mass of about a few × 10−8M* to the N-body simulation with a rate and position set by F0(ε/f), where f is the infall gas-to-solids ratio and ε is the fraction of infalling solids that are incorporated into planets36. This is consistent with our assumption that planetesimals form interior to rc, discussed next.

Dust-to-planetesimal accretion

Formation of  >km-sized planetesimals is thought to involve: 1) growth of macroscopic particles (pebbles) via grain-grain collisions, 2) mechanism(s) to spatially concentrate pebbles, and 3) gravitational collapse to yield planetesimals. While we do not model the dust-to-planetesimal stages of growth here, in this section we make the case that planetesimal formation interior to rc is plausible given current understanding.

Collisions between orbiting dust grains initially lead to rapid growth23. A particle’s Stokes number, St, is the product of its orbital frequency and its stopping time due to gas drag, and St increases with particle size with a functional dependence that varies with drag regime.

As particles grow, their collision velocities increase due to size-dependent interactions with the gas disk. The Stokes number of the largest particles produced via grain collisions, \({{{\rm{St}}}}_{\max }\), can be estimated by equating the collision velocity to a critical velocity, vf, above which collisions lead to fragmentation rather than to growth. With relative velocities dominated by turbulent motions of the gas (and assuming that the α viscosity parameter characterizes the diffusion of both gas and solids), this gives e.g.,24

$${{{\rm{St}}}}_{\max }\approx {({v}_{f}/{c})}^{2}/(3\alpha )$$
(6)

with gas sound speed c and viscosity parameter α. \({{{\rm{St}}}}_{\max }\) is strongly dependent on vf, which is an uncertain parameter. For water ice, multiple works consider an enhancement of vf either near the condensation line of water63, or due to the addition of a thin layer of water molecules64. For rocky grains, multiple works adopt vf = 1 m/s [e.g.,23,49,65], but there is a wide range of reported values. Silicate dust collision experiments at room temperatures66, granular mechanics simulations67, and observations of a young massive disk too hot for ices68 imply vf values in the 10 to 30 m/s range. Further, new experiments suggest that vf is increased as temperatures approach the rock melting point (which perhaps is not surprising given increased “stickiness” of ice-rich grains known to occur near the ice condensation line). A key quantity in estimating the energy needed to break grain contacts is the surface energy69, γe, with \({v}_{f}\propto {\gamma }_{e}^{5/6}\). Experiments show that as temperature increases from 1000 K to 1400 K25, the effective surface energy of silicate grain aggregates increases by 2 to 3 orders-of-magnitude, implying vf ≈ 50 to 300m/s in the T ≳ 1000 K disks modeled here. Thus hot inner disks may be a favored region for grain growth, and plausible combinations of vf and α can yield \(1{0}^{-2} \, < \, {{{\rm{St}}}}_{\max } \, < \, 1{0}^{-1}\), with

$${{{\rm{St}}}}_{\max }\approx 0.05{\left(\frac{{v}_{f}}{30{{\rm{m}}}/{{\rm{s}}}}\right)}^{2}\left(\frac{1200\,{{\rm{K}}}}{T}\right)\left(\frac{0.001}{\alpha }\right)$$
(7)

Assuming the Epstein drag regime, the corresponding maximum particle size would be

$${a}_{\max }\approx 10\,{{\rm{cm}}}\left(\frac{{\sigma }_{g}}{1{0}^{3}\,{{\rm{g}}}\,{{{\rm{cm}}}}^{-2}}\right){\left(\frac{{v}_{f}}{30{{\rm{m}}}{{{\rm{s}}}}^{-1}}\right)}^{2}\left(\frac{1200\,{{\rm{K}}}}{T}\right)\left(\frac{0.001}{\alpha }\right).$$
(8)

For comparison, prior works often adopt \({a}_{\max } \sim 10\) cm for icy grain growth in the outer disk, due to the increased stickiness of warm ice near the snowline49. For rocky grains interior to the snowline, works that adopt vf = 1 m/s find \({a}_{\max } \sim\) mm49, while those that adopt vf = 10 m/s predict \({a}_{\max }\, \gtrsim\) cm70 or \({{{\rm{St}}}}_{\max }\, \gtrsim \,O(1{0}^{-1})\)71 in inner disk regions. The latter work notes that silicate dust may enter into a plastic regime at temperatures  > 1000 K that renders every collision more dissipative and increases vf, similar to the effect observed in recent experiments of ref. 25.

Planetesimal formation by steaming instability (SI) requires that the volume density ratio of solids-to-gas exceeds unity in the mid-plane; this condition can be alternatively expressed as a minimum solid-to-gas surface density ratio, Zcrit. Determining Zcrit, which depends on St (and also on α if the disk is turbulent), requires state-of-the-art hydrodynamical models with sufficient resolution and run time. Recent high-resolution simulations of laminar disks29 find 0.005 < Zcrit < 0.006 for St = 10−1, and 0.01 ≤ Zcrit < 0.02 for St = 10−2, implying that for St ≥ 10−2, collapse may occur with stellar metallicities or only modest solid concentration. Gas turbulence and particle self-gravity increase Zcrit for a given St. Recent simulations including these effects30 for St = 10−1 find Zcrit ≈ 0.02 [≈0.06] for α = 10−4 [α = 10−3]; with St = 10−2, the simulations find Zcrit ≈ 0.06 [≈0.2] for α = 10−4 [α = 10−3]. These Zcrit values should be viewed as upper limits, because SI and collapse tend to become more likely as numerical resolution is increased30.

Several mechanisms could enhance Z in the inner disk. Inward radial drift of particles formed at larger orbital radii increases Z in the inner disk due to its smaller cross-sectional area31; e.g., dust evolution models23, their Fig. 8 that include infall find this produces Z ≈  0.02 in the disk interior to 1 au. Larger solid enhancements are possible if particles are concentrated by additional effects. Hydrodynamical simulations show that during infall, vortices due to Rossby wave instability form near rc, trapping particles and enhancing the local dust-to-gas ratio by a factor of approximately 5 [10] for 1-cm [10-cm] particles33.

Another proposal7,72 is that planetesimal formation will be favored in regions where the gas radial velocity is outward and advection of particles by the gas can counter their inner radial decay, leading to a force balance for certain particle sizes that may then become concentrated. In an infall-supplied disk, the gas velocity becomes outward at a location interior to rc that is a function of the infall angular momentum73; for the uniform infall flux adopted here, this occurs for r ≳ 0.85 rc for rd/rc ≳ 10, where rd is the disk’s outer radius, so that vg > 0 may apply to the outer  ≈20% of the disk area interior to rc.

Additionally, several works propose that solids would be preferentially concentrated near the rock evaporation radius, revap, that falls within rc in our model74,75. Interior to revap, the disk’s effective α value may greatly increase due to magnetically-driven turbulence, creating a local pressure maximum that could trap particles and strongly increase Z35. Alternatively, refs. 7,49 argue that if the gas velocity is outward across revap, the vapor produced by evaporation of inwardly drifting grains would be carried back outward, allowing for recondensation and a preferred region for solid build-up. This effect would apply to infall-supplied disks if revap falls within the outer infall region described above.

Thus while tremendous progress in understanding planetesimal formation has been made, substantial uncertainties remain. Multiple concepts identify disk locales at which solid concentration and planetesimal formation may be preferred (see, e.g., review70), but models that explicitly treat grain growth up to planetesimal formation are lacking, and whether conditions for collapse are actually achieved is unclear, particularly for turbulent disks. The appropriate value of vf in different regimes and/or the importance of an additional bouncing barrier76 remain an important topic of debate.

However, observations of young disks increasingly provide evidence of early accretion. Recent developments highlighted above suggest 1) a higher vf in inner disks than previously assumed; 2) several mechanisms capable of substantially concentrating solids in inner disks; and 3) a widening of the parameter space that allows for clumping by SI even for St ≪ 1, due to higher-resolution, longer time hydrodynamical simulations. These results suggest planetesimal formation interior to rc during final infall is plausible.