A semi-analytic model of the bouncing barrier for protoplanetary dust aggregates

Abstract: Collisional bouncing limits the growth of dust aggregates in protoplanetary disks, but its dependence on aggregate size, collision velocity, and filling factor remains poorly understood. Here we develop a semi-analytic model for the sticking probability of colliding dust aggregates. We divide each aggregate collision into two phases: a compression phase and a separation phase. The compression phase is described with an elastoplastic contact model, which determines the maximum contact radius and repulsive energy after compression. The separation phase is treated as fracture of a stochastic network of interparticle bonds, whose fracture energy is evaluated using weakest-link statistics. The model naturally predicts that larger aggregates bounce more readily because larger contact regions are more likely to contain weak bonds. Comparison with distinct element method simulations shows that the model reproduces the simulated sticking--bouncing boundary. Furthermore, applying the calibrated model to moderately porous aggregates inferred from ALMA observations of protoplanetary disks, we find that the predicted bouncing barrier passes through the observationally inferred size--velocity range. Thus, our semi-analytic model provides a useful framework for predicting the collisional evolution of protoplanetary dust aggregates.

• The paper develops a semi-analytic sticking model that combines elastoplastic contact mechanics during compression with weakest-link fracture statistics during separation.
• The model predicts that larger aggregates bounce more readily because heterogeneous bond networks make fracture energy grow more slowly than aggregate mass, and identifies α = 2 and β ≈ 1.3 as the best-fit parameters against DEM simulations.
• Applying the model to porous aggregates in the IM Lup disk predicts a 50% sticking radius near the observed 0.1–1 mm range at collision speeds below 1 m s⁻¹, while highlighting uncertainties from porosity, composition, oblique impacts, and erosion.

Overview and motivation

Collisional coagulation of submicron dust grains in protoplanetary disks is limited by the bouncing barrier, but the conditions separating sticking from bouncing have remained uncertain, particularly because DEM simulations typically use smaller aggregates than laboratory experiments. Arakawa et al. (2606.26491) develop a semi-analytic model for the sticking probability of colliding dust aggregates that treats the collision as two phases: a compression phase described by elastoplastic contact mechanics, and a separation phase treated as fracture of a stochastic network of interparticle bonds evaluated with weakest-link statistics. The central result is a natural explanation for the size dependence of the bouncing barrier: larger contact regions are statistically more likely to contain weak bonds, so fracture energy grows more slowly than aggregate mass.

Model structure

The model considers head-on collisions between equal-mass aggregates of radius $R_{\rm agg}$, filling factor $\phi$, and monomer radius $r_1$. In the compression phase, the aggregates are treated as effective elastoplastic spheres with reduced Young's modulus $E^*$ and yield stress $\sigma_y$, following Andrews' classical contact model [doi:10.1080/14786443008565033]. This yields the maximum compression length $\delta_{\rm max}$, maximum contact radius $a_{\rm max} = \sqrt{R^* \delta_{\rm max}}$ (with $R^* = R_{\rm agg}/2$), and repulsive energy $K_{\rm rep}$ available after maximum compression. Because $\delta_{\rm crit} \propto R^*$ and $\phi$0 at fixed velocity, the contact geometry scales self-similarly: $\phi$1, consistent with prior DEM results [2024arXiv240815573A, 2025ApJ...995..207T].

In the separation phase, the contact region is idealized as $\phi$2 chains of $\phi$3 bonds connected in parallel, with $\phi$4 and $\phi$5, where $\phi$6 is the mean coordination number. Each bond has a heterogeneous breaking energy drawn from a power-law cumulative distribution $\phi$7 for $\phi$8, where $\phi$9 combines the normal bond-breaking energy $r_1$0 and rolling dissipation energy $r_1$1. Weakest-link statistics give a mean chain fracture energy $r_1$2 — the key nonlinearity: chain strength does not scale linearly with length when bonds are heterogeneous. Summing over parallel chains and applying the central limit theorem, the sticking probability is

$r_1$3

so the 50% sticking condition corresponds to $r_1$4.

Scaling relations for size and velocity dependence

The model yields closed-form scaling predictions. At fixed velocity and filling factor, $r_1$5 while $r_1$6, so

$r_1$7

which directly implies that larger aggregates bounce more readily. In the elastic regime ($r_1$8), the ratio scales as $r_1$9, giving a 50% sticking radius $E^*$0; for $E^*$1 this is $E^*$2. In the plastic regime ($E^*$3), plastic dissipation suppresses $E^*$4 while the growing contact area raises $E^*$5; for $E^*$6 the ratio grows with velocity as $E^*$7, predicting that bouncing is favored at intermediate rather than highest velocities — consistent with the DEM findings of Oshiro et al. [2025ApJ...983...75O]. This non-monotonic velocity dependence is a distinctive prediction of the framework.

Calibration against DEM simulations

Comparing against the DEM simulations of Oshiro et al. [2025ApJ...983...75O] for icy aggregates with $E^*$8 nm at $E^*$9 and $\sigma_y$0, the authors find that the sticking–bouncing boundary depends strongly on both $\sigma_y$1 and $\sigma_y$2. For $\sigma_y$3, $\sigma_y$4 systematically underestimates $\sigma_y$5 (predicted $\sigma_y$6 too small), while $\sigma_y$7 overestimates it; $\sigma_y$8 gives broad agreement. With $\sigma_y$9 fixed, $\delta_{\rm max}$0 and $\delta_{\rm max}$1 fail, whereas $\delta_{\rm max}$2 reproduces the simulated boundary for both filling factors. Notably, for $\delta_{\rm max}$3 or $\delta_{\rm max}$4, no single value of $\delta_{\rm max}$5 simultaneously fits both $\delta_{\rm max}$6 and $\delta_{\rm max}$7, so the combination $\delta_{\rm max}$8 is uniquely selected among the explored parameter sets. Physically, $\delta_{\rm max}$9 implies that grains undergo more than one equivalent 90° rolling displacement before bonds break during separation.

The physical interpretation of $a_{\rm max} = \sqrt{R^* \delta_{\rm max}}$0 is plausible given that energy dissipation during separation occurs mainly through grain rearrangement via rolling motion, but its status as an effective fitting parameter should be kept in mind: it is calibrated to one set of DEM simulations rather than derived from first principles.

Application to observed protoplanetary disks

Applying the best-fit parameters to moderately porous aggregates ($a_{\rm max} = \sqrt{R^* \delta_{\rm max}}$1–$a_{\rm max} = \sqrt{R^* \delta_{\rm max}}$2) inferred from ALMA radiative-transfer modeling of the IM Lup disk [2024NatAs...8.1148U] — which suggests $a_{\rm max} = \sqrt{R^* \delta_{\rm max}}$3–1 mm, $a_{\rm max} = \sqrt{R^* \delta_{\rm max}}$4–$a_{\rm max} = \sqrt{R^* \delta_{\rm max}}$5, and $a_{\rm max} = \sqrt{R^* \delta_{\rm max}}$6 — the model predicts that the 50% sticking radius falls within the observationally inferred size range at velocities below 1 m s⁻¹. This agreement indicates that the observed (sub)millimeter-sized aggregates may be sized by the bouncing barrier itself, providing an observational consistency check that purely numerical studies could not reach. A significant advantage of the semi-analytic formulation is that it extends to millimeter-sized or larger aggregates that are computationally inaccessible to direct DEM simulation.

Limitations and open questions

The paper states several caveats explicitly. First, the parameters $a_{\rm max} = \sqrt{R^* \delta_{\rm max}}$7 and $a_{\rm max} = \sqrt{R^* \delta_{\rm max}}$8 were constrained using relatively compact aggregates ($a_{\rm max} = \sqrt{R^* \delta_{\rm max}}$9–$R^* = R_{\rm agg}/2$0); their transferability to the more porous $R^* = R_{\rm agg}/2$1–$R^* = R_{\rm agg}/2$2 regime relevant to IM Lup is untested and requires new DEM simulations. Second, the model assumes equal-mass, head-on collisions of non-sintered icy aggregates; impact parameter, mass ratio, aggregate shape, composition, and sintering are all neglected. Third, the model addresses only the sticking–bouncing outcome and does not treat erosive mass loss, which may shape size distributions over many collisions. Finally, the elastoplastic contact model assumes small deformations, so predictions for $R^* = R_{\rm agg}/2$3 (where $R^* = R_{\rm agg}/2$4) carry uncertain validity; the fracture-energy estimate also fails formally when $R^* = R_{\rm agg}/2$5. An important open question left by the paper is whether silicate dust aggregates — directly relevant to terrestrial planet formation — share the same $R^* = R_{\rm agg}/2$6 values, which the authors propose to test against the extensive existing laboratory data.

Conclusion

This work provides a compact, physically motivated analytic description of the bouncing barrier in which the size dependence arises from weakest-link statistics of heterogeneous interparticle bonds. Calibrated with only two parameters against DEM simulations, the model reproduces the simulated sticking–bouncing boundary and, when extrapolated to porous aggregates, yields a barrier consistent with ALMA-inferred aggregate properties in the IM Lup disk. Its main limitations — untested parameter transferability across porosities and compositions, and neglect of erosion and oblique impacts — define the immediate agenda for extending the framework toward predictive dust-growth models.