A Pathway for Collisional Planetesimal Growth in the Ice-dominant Regions of Protoplanetary Disks

We present a semi-analytic model for the growth, drift, desorption, and fragmentation of millimeter- to meter-sized particles in protoplanetary disks. Fragmentation occurs where particle collision velocities exceed critical fragmentation velocities. Using this criterion, we produce fragmentation regions in disk orbital radius–particle size phase space for particles with a range of material properties, structures, and compositions (including SiO2, Mg2SiO4, H2O, CO2, and CO). For reasonable disk conditions, compact aggregate H2O, CO2, and CO ice particles do not reach destructive relative velocities and are thus not likely to undergo collisional fragmentation. Uncoated silicate particles are more susceptible to collisional destruction and are expected to fragment in the inner disk, consistent with previous work. We then calculate the growth, drift, and sublimation of small particles, initially located in the outer disk. We find that ice-coated particles can avoid fragmentation as they grow and drift inward under a substantial range of disk conditions, as long as the particles are aggregates composed of 0.1 μm-sized monomers. Such particles may undergo runaway growth in disk regions abundant in H2O or CO2 ice, depending on the assumed disk temperature structure. These results indicate that icy collisional growth to planetesimally relevant sizes may happen efficiently throughout a disk’s lifetime, and is particularly robust at early times when the disk’s dust-to-gas ratio is comparable to that of the interstellar medium.


INTRODUCTION
In the classical picture of particle evolution during the initial stages of planet formation, collisional growth of solids to the planetesimal size scale is thought to be prevented by several limitations that are jointly known as the "meter-size barrier".While other processes may be relevant to this problem (discussed in Section 8.2), the classical meter-size barrier is the barrier to continued particle growth beyond millimeter to sub-kilometer scales due to either collisional fragmentation or particle drift (see Chiang & Youdin (2010); Birnstiel et al. (2016) for reviews).Fragmentation occurs when particles collide at relative velocities that are energetic enough to exceed the surface energy keeping the particles intact.The particle-particle relative velocities are maximized in the meter-size regime at separations of ∼1AU from the host star, resulting in potentially destructive collisions (Blum & Wurm 2000).These particles also become large enough to begin to decouple from the dynamics of the background gas, such that significant gas drag causes particles to rapidly drift inwards towards the system's host star.The effects of both collisional fragmentation and inward drift prevent particles in this regime from collisionally growing to large sizes.
Given the challenges of continued collisional growth past the meter-size barrier, alternative models for the formation of ∼1-100km planetesimals have been proposed.Resonant Drag Instabilities (RDIs), such as the streaming instability, have been shown to create particle overdensities which can directly gravitationally collapse into planetesimals (e.g., Goldreich & Ward 1973;Youdin & Shu 2002;Youdin & Goodman 2005;Johansen et al. 2006b;Chiang & Youdin 2010;Simon et al. 2016;Squire & Hopkins 2018;Gerbig et al. 2020).While models of unstable RDIs successfully produce planetesimals, they rely on strict initial conditions that may not represent commonly occurring disk conditions.
Enough uncertainty remains, however, in the behavior of particle collisions in disks that collisional growth may still be a viable planetesimal formation mecha-nism.While some of this uncertainty comes from the protoplanetary disk conditions during formation, crucially, there are significant uncertainties in the material properties of solids.These include, but are not limited to, the ice-coated particles composition, structure (i.e., non-aggregate vs. aggregate with small or large monomers, sintered vs. unsintered monomer connections, crystalline vs. amorphous, etc.), physical surface processes between silica-ice and ice-ice boundaries (Fogarty et al. 2010;Nietiadi et al. 2020), and compaction behavior of porous particles upon impact (Paszun & Dominik 2009;Krijt et al. 2015).
Particles beyond the water ice line are expected to be coated in ices with material properties that differ from silicate grains, thus altering the particles' strength in withstanding fragmentation.For example, several laboratory studies have shown that the critical velocity for H 2 O ice, and grains which are coated in ice may be higher than that for silicates, indicating that ice-coated particles may be more robust against collisional fragmentation (Poppe et al. 2000;Blum & Wurm 2008;Gundlach et al. 2011;Gundlach & Blum 2015;Musiolik et al. 2016a,b; see Section 2 for further discussion).This is thought to be largely due to the increased surface energy of H 2 O ice, particularly at temperatures near the ice line (Gundlach et al. 2018;Musiolik & Wurm 2019).Musiolik et al. (2016a) investigate CO 2 ice particles and find that mixtures between CO 2 and H 2 O can increase the overall sticking during collisions.In terms of particle structure, Kataoka et al. (2013) find that micronsized icy dust aggregates may stick together forming fluffy planetesimals that are able to surpass the drift barrier without fully fragmenting, and numerical models by Wada et al. (2007Wada et al. ( , 2009) ) reveal that icy dust aggregates are able to retain material during collisions of varying impact parameters and velocities.Collision-induced heating of the ice particles may also be an important factor that changes the surface layer physics from dry to wet (Nietiadi et al. 2020).Ultimately these studies demonstrate that particle composition (both in grain and ice mantles) and the corresponding particle material properties are of first order importance when modeling particle evolution in disks.
In this work, we investigate the implications of a species-dependent fragmentation velocity on the growth and evolution of particles in disks.We build upon the critical fragmentation velocity frameworks from Wada et al. (2007Wada et al. ( , 2009) ) and Stewart & Leinhardt (2009) for a variety of particle compositions, porosities, and material properties.This model calculates particle relative velocities based on the relevant drag regimes present throughout the disk (i.e., Epstein, Stokes, and ram pres-sure drag regimes), and finds regions of fragmentation as a function of particle size and orbital radius by comparing the relative velocities with the respective critical fragmentation velocities.
We produce the expected fragmentation regions for SiO 2 , Mg 2 SiO 4 , H 2 O, CO 2 , and CO solid particles in TW Hya and in the Minimum Mass Solar Nebula (MMSN).We investigate the consequences of assumptions regarding particle material properties and gas drag regimes on the likelihood of particle fragmentation and growth, and discuss whether particles are strongly affected by the classical meter-size barrier.We find that compact aggregate H 2 O, CO 2 , and CO ice particles in the outer regions of protoplanetary disks may not undergo collisional fragmentation.Particle drift and desorption of ice are incorporated in the model such that the radial particle evolution and growth is compared to the fragmentation regions.We describe favorable locations for continued particle growth in the disk typically ∼1-10 AU, either near the CO 2 ice line or exterior to the H 2 O ice line depending on the assumed disk temperature profile.These regions may allow for efficient runaway solid growth that is not limited by fragmentation, drift, or desorption, where the classical meter-size barrier is unlikely to operate.
In Section 2 we discuss the relevant material properties that govern particle collisions and in Section 3 derive the corresponding critical velocities used in this study.In Section 4 we outline the surface density profile and temperature dependencies in protoplanetary disks used in the model.Particle relative velocities are then derived in Section 5.The resulting fragmentation region phase space is described in Section 6, with favorable regions for planetesimal formation detailed in Section 7. Implications of these results are discussed in Section 8, including other potential barriers to particle growth.A summary and conclusion of this work can be found in Section 9.

PARTICLE PROPERTIES AND GROWTH
We model the fragmentation of key volatiles in protoplanetary disks including H 2 O, CO 2 , and CO, as well as two rock species that are only volatile in the very innermost high-temperature regions of the disk.While previous studies have frequently focused on the fragmentation of silicates, such as SiO 2 , we note that in comets the majority of non-volatile solid material has a composition similar to that of Mg 2 SiO 4 (Wooden et al. 2007).We thus consider the fragmentation properties of both of these species.The following properties determine the strength of an ice-aggregate, thus determining the likelihood that a particle will fragment upon collision: surface energy (γ), Young's modulus (E), Poisson's ratio (ν), and material density (ρ i ); material properties for each ice species are given in Table 1.Young's modulus quantifies how deformable and stretchy a material is-the higher the value the less deformable the material (Heindl & Mong 1936).Poisson's ratio quantifies how much a material will stretch in the perpendicular direction of the impacting load, typically ranging from 0 to 0.5 (Greaves et al. 2011).The size of the monomers composing the ice-aggregate also determines the overall strength and is discussed in context of the material properties below.We further include the binding energy of H 2 O, CO 2 , and CO in Table 2 which, in this context, is the energy required to remove ice from the surface of a grain.In this work the binding energy determines the ice desorption rate and thus sets the ice line location for each species (see Section 4).A species' surface energy is the amount of energy per unit area required to make (or destroy) a surface, and is a crucial component in calculating the critical fragmentation velocity in our model.
References-a Sandford & Allamandola (1988) b Aikawa et al. (1996) c Pontoppidan (2006) Both the binding energy and surface energy can be calculated from either laboratory experiments or numerical models and a range of values exist for these quantities in the literature.In this work, we have conservatively adopted frequently chosen literature values though we note that other values for these parameters are sometimes available and can be different for crystalline versus amorphous molecular structures.In particular, we note that recent laboratory experiments have found that the surface energy of H 2 O, particularly at temperatures cooler than ∼175 K, is significantly lower than previous estimates based on studies performed at higher temperatures (Gundlach et al. 2018;Haack et al. 2020;Musiolik 2021).However, including these lower surface energy estimates in models of disks, particularly when there is significant turbulence, can inhibit the formation of ∼millimeter-sized particles needed to reproduce observed millimeter disk fluxes (Pinilla et al. 2020).This could potentially indicate that the added heat from the kinetic energy of collisions can partially melt the H 2 O ice mantle, causing an increase in ice stickiness, even at lower temperatures (Nietiadi et al. 2020).Given these uncertainties, future laboratory H 2 O ice collision experiments at relative velocities and low temperatures relevant for protoplanetary disks, would provide ideal constraints.
CO 2 ices also have a range of estimated surface energies from numerical and laboratory work (noting that the value we use in this study, γ CO2 = 60 erg cm −2 , has a large uncertainty of ±22 erg cm −2 ).For example, Musiolik et al. (2016a) find that pure CO 2 ice has a critical fragmentation energy comparable to that derived from the tensile strengths of the silicate-containing mixtures basalt and palagonite.However, Musiolik et al. (2016b) find that ice mixtures can have significantly higher surface energies than the value for CO 2 used in this work, which can change the collisional behaviour from a homogeneous ice mantle (Musiolik et al. 2016a).Fewer Figure 1.Visual representation of the three ice coating growth scenarios (see Section 2 for a detailed description).In all cases, particles coagulate until they reach a size where ice can form a stable coating on the surface.In the first case (a), the volatile material thinly coats aggregate particles that then undergo growth via coagulation to form larger ice-coated particles with material properties representative of an ice-coated aggregate with small monomer sizes.In the second case (b), the volatile gas is more abundant and able to quickly form particles with a large coating of ice.These particles then coagulate and have fragmentation properties representative of an ice-coated aggregate with large monomer sizes.In the third case (c), a limited number of large particles rapidly form significant ice coatings and grow via condensation.These particles have fragmentation properties representative of a non-aggregate silicate or icy particle .recent constraints exist for the material properties of CO ice, and as such we match and assume the Young's modulus and Poisson's Ratio to that of water as it is in a reasonable range for ice.In general, surface energies in the conditions relevant in protoplanetary disks are not well-constrained as these values are typically based on laboratory experimental setup or idealized numerical computations.
Furthermore, the structure of ice particles in disks is largely unconstrained.We thus detail below three physically-motivated cases that may occur in protoplanetary disks and lead to significantly different predictions for the size at which particles fragment (see Figure 1).These cases correspond to particles that have fragmentation properties that resemble those of: (a) an ice-coated aggregate comprised of small monomers, (b) an ice-coated aggregate comprised of large monomers, and (c) an ice-coated compact particle.We further discuss a corollary potential case of an aggregate particle behaving like an ice-free aggregate composed of small monomers.There are several different potential outcomes of ice formation depending on the relative solid to volatile gas abundance, which changes throughout the lifetime of the disk.We assume that either the particles are coated in ice beyond ice lines in the disk, or ice surfaces formed on particles within the molecular cloud from which the star and disk formed and are retained beyond ice lines (Furuya et al. 2016;Oberg & Bergin 2021).Uncertainty remains in how reprocessed these inherited particles from the molecular cloud are, and as such we make no assumption on the composition and structure of particles prior to the disk phase.
In the first potential case (Figure 1a), a small coating of ice covers the entirety of the aggregate particle, including the surface of the grain and in the cavities of the aggregate particle.This scenario is likely to occur when there is a sufficient number of particles that can quickly deplete the available supply of volatile gas (or if the gas is only marginally supersaturated) such that each particle is coated with a relatively thin layer of ice.Furthermore, experiments of ice formation on Earth indicate that ice preferentially forms in grain cavities and is thus likely to coat the particle as pictured in Figure 1a (e.g., Campbell et al. 2017;Campbell & Christenson 2018;Holden et al. 2021).Once ice formation depletes the supply of volatile material, the abundant large, ice-coated aggregates will undergo further growth via coagulation.This growth is likely to occur with icecoated aggregates with similar properties (e.g., Powell et al. 2019).The resultant large particles that are produced via this growth scenario will then resemble an aggregate particle comprised of small monomers where the connections between the grains are dominated by the material properties of the dominant ice species coating the surface.This case may resemble the structure of a sintered aggregate, however the connection strength is of different compositions (rocky or icy) and the collision outcome between the sintered and unsintered cases will also differ (further discussed in Sections 3 and 8.2).
In the second potential case (Figure 1b), a large coating of ice will cover the surface of the silicate grains.This scenario is likely to occur if there is an abundant supply of volatile gas than in the first case but there still exists a large number of sufficiently large aggregates that have become ice-coated.In this case, once the particles have depleted the available supply of volatile gas, further particle growth will be dominated by collisions with particles of similar sizes.Due to the large coating of ice on each particle, these collisions will likely result in fragmentation as they are structurally similar to ice-coated aggregates comprised of large monomers (discussed in detail in Section 8.1).
In the third potential case (Figure 1c), ice will form at a very efficient rate and will quickly dominate the particle's properties.This scenario is likely to occur if volatile gas is abundant and there is a limited supply of large particles.In this case, the limited number of large particles will abundantly form ice and are unlikely to grow further via coagulation due to their low relative number densities.The fragmentation properties of these particles are likely to resemble that of a compact icy grain and not those of an aggregate particle.
In the ice free case, the connections between the monomers are dominated by the material properties of the ice-free monomer cores which are likely composed of grains that are often silicate in composition.This case is likely for small particles and may also occur if ice formation is inefficient or the volatile gas supply is depleted.
These different structural cases will lead to varying critical fragmentation velocities, ultimately depending on not only the particle composition, but also the particle structure and in particular the size of monomers composing the aggregates (see Figures 2 and 11 for specifics).For this study, we will evaluate the particle fragmentation for varying compositions, while taking the nominal particle structure to be compact ice-coated unsintered aggregates composed of small 0.1 µm sized monomers.This monomer size corresponds to the average particle size expected in the interstellar medium (Blum et al. 2006;Gundlach et al. 2018).We take the first case to be our nominal case because large coatings of volatile ice on particles require significant supersaturations of condensible gas which likely only exists in disks at early times or near a species' iceline (Powell et al. 2022).A comparison of how particle structure and monomer size is demonstrated in reference to the critical fragmentation velocities in the following Section (Section 3) with Figures 2 and 3.

CRITICAL VELOCITY FOR COLLISIONAL DESTRUCTION
Particle fragmentation, in the context of the fundamental material properties and structures discussed in Section 2, occurs in qualitatively different ways for solid rocks versus porous aggregates.For non-aggregate solid rocks with minimal porosity (as in Figure 1c), cracks typically develop at weak points where atoms are less effectively bound to their neighbors.A sufficiently substantial crack can weaken a solid rock enough for it to fragment into pieces.For porous rocks or aggregates (as in Figure 1a,b), the connections between distinct, solid monomers are individually weaker than the bonds between component atoms within a monomer, thus fragmentation typically results from breaking connections between monomers.Perhaps counter intuitively, this behavior causes aggregates composed of very small monomers to be more resilient to destructive collisions than solid compact particles since collisions between porous grains tend to compact the aggregate rather than fully breaking it apart along cracks (see Seizinger et al. 2013, in the case of erosion).To break an aggregate, the majority of the connections making up the full aggregate must be broken.Sintered aggregates may be stronger to collisions since the connections become fused into a neck as compared to a contact point, but whether the collision will result in growth versus bouncing depends on how compact the aggregate is and where in the disk it is (e.g., Maeno & Ebinuma 1983;Blackford 2007;Sirono 1999;Sirono & Ueno 2017;Sirono & Kudo 2021, see Section 8 for a discussion of the sintered case).In protoplanetary disks, particles may begin their growth as fluffy aggregates (e.g., Smirnov 1990;Dominik & Tielens 1997;Blum & Wurm 2000), but the size at which growing bodies transition to low-porosity solids remains unclear (Kataoka et al. 2013).We therefore consider both regimes of particle porosity and structure.
For fragmentation of solid non-aggregate particles, we expand on the critical disruption criterion presented in Stewart & Leinhardt (2009).The kinetic energy per mass of a projectile required to destroy a target planetesimal in either the strength or gravity dominated regimes is The term on the left, , represents the strength regime (small particles bonded together), while the term on right, q g R 3µ C1 , represents the gravity regime (rubble piles gravitationally held together).Here q s , q g , µ, and ϕ define the material properties in cgs units; values for strong and weak rock are in Table 3 (see Housen & Holsapple 1990, 1999;Stewart & Leinhardt 2009, for details).The sum of the spherical radii of the target and projectile is R C1 , which we approximate as the size of the larger colliding particle.The disruption criterion for our particles is not expected to be dominated by the gravity regime, however it is still included in the model.The disruption criterion can be described as the binding energy keeping together the compact particle (or in the gravity regime, a porous rubble pile consisting of a gravitational aggregate of rigid spheres) (Stewart & Leinhardt 2009).Comparing this energy with the collisional kinetic energy of a projectile provides the critical velocity needed to make cracks in and ultimately break apart the target particle.Here the kinetic energy of the projectile is KE = (1/2)m p v 2 rel , while the binding energy of the target particle is BE = m t Q * RD , where m p and m t are the masses of the projectile and target respectively.Converting the masses in terms of the target radius, R C1 , and solving for the relative velocity v rel , we find that the critical relative velocity required to break apart a solid compact particle is where f is the size ratio between the two colliding particles.It should be mentioned that for solid compact particles, relative velocities that are nearly at critical fragmentation velocities may also result in bouncing rather than sticking, which is further discussed in Section 8.For fragmentation of small aggregate bodies (not gravitational rubble piles), we follow the work of Wada et al. (2007Wada et al. ( , 2009) ) to find the critical velocity for collisional destruction.Their numerical studies model the collisions between aggregates in the ballistic cluster-cluster aggregation (BCCA) and ballistic particle-cluster aggregation (BPCA) regimes, or more generally fluffy and compact aggregates.This kind of fragmentation focuses on pulling apart each monomer, in contrast to the solid fragmentation framework of Stewart & Leinhardt (2009).The energy required to break apart two contact monomers is where F c = 3πγR is the maximum force required to separate the contact and δ c = (9/16) 1/3 a 2 0 /(3R) is the critical separation or compression distance between the two monomers (denoted with numerical subscripts).These expressions depend on the contact circle radius a 0 = (9πγR/E * ) 2/3 (Wada et al. 2007) ,  1/R = 1/r m,1 + 1/r m,2 with r m,1 and r m,2 being the  monomer radii, and 1 Wada et al. 2009).All of these are a function of the material properties-Young's modulus (E), Poisson's ratio (ν), surface energy (γ), and the material density (ρ i ) discussed in Section 2. The kinetic energy of a colliding aggregate is KE = (1/2)N total m m (v rel /2) 2 , where N total is the total number of particles composing both aggregates and m m = (4/3)πr 3 m ρ i is the mass of a monomer.The critical impact energy is E crit = kN total E break , where k is a dimensionless factor which takes into account whether the aggregate is fluffy or compact.Wada et al. (2009) find that for fluffy aggregates (BCCA) k ∼ 10, while for compact aggregates (BPCA) k ∼ 30.The critical relative velocity is then found by balancing the kinetic energy with the critical impact energy and solving.
Rolling friction is not included in our model as it may not impact the outcome of collisions as demonstrated by Arakawa et al. (2022a).We also note that viscous energy dissipation through the particle may not necessarily break apart every monomer connection in an aggregate at this particular critical fragmentation velocity as other tangential forces and impact parameter may determine the collision outcome depending on the assumed viscous dissipation timescale (e.g., Arakawa et al. 2022b).This framework does however provide an intuition for the compositional dependencies on fragmentation.
Figure 2 compares the calculated critical velocities for non-aggregate and aggregate particles, highlighting the impact of monomer size.The critical velocities for aggregate silicates composed of SiO 2 or Mg 2 SiO 4 , with monomer sizes in the range of roughly 0.1 to 1 microns, are comparable to those of non-aggregate strong and weak rock as demonstrated in Figure 2. As the monomer size increases the critical velocity for aggregates monotonically decreases, indicating that the choice of monomer size is an important parameter in modelling fragmentation (results discussed in 8.1).As mentioned at the end of Section 2, we choose a fiducial monomer size of 0.1 µm, consistent with particles in the ISM (Oberg & Bergin 2021) The difference in surface energy is also demonstrated as the critical velocities are consistently lower for SiO 2 than Mg 2 SiO 4 , where SiO 2 has a much lower surface energy (see Table 1).Figure 3 compares the differing critical fragmentation velocities between fluffy BCCA and compact BPCA aggregate particles with 0.1µm monomers and varying compositions.Compact BPCA H2O ice is the strongest aggregate.Because strength is determined by surface interactions between monomers, silicates coated in ice will have strengths determined by that ice.Weidenschilling (1977a) and Hayashi (1981): Protoplanetary disks are defined as either passively or actively heated depending on whether accretion onto the host-star is present.In a passively heated disk, stellar irradiation is the dominant heat source for the entire disk, where incoming radiation from the host star is absorbed in local regions of the disk and then re-emitted as a blackbody.An actively heated disk includes viscous heating caused by viscous midplane accretion of the disk onto the host star.If viscous accretion heating is included, it will be the dominant heat source in the inner disk while the mid to outer disk remains dominated by irradiation heating.Throughout a disk's lifetime, the disk can either remain fully passive, fully active, or switch from one to the other (Armitage 2017).
From Chiang & Goldreich (1997), the irradiation temperature follows a power-law depending on the stellar mass (M * ) and luminosity (L * ) where σ SB is the Stefan-Boltzmann constant, k is the Boltzmann constant, G is the gravitational constant, µ = 2.3m H is the mean molecular weight assuming a hydrogen-helium composition, and m H is the mass of a hydrogen atom.Using values from Powell et al. (2019), TW Hya has M * = 0.8M ⊙ , L * = 0.28L ⊙ , and T o = 82 K.The MMSN has T o = 120 K Chiang & Goldreich (1997).
The temperature due to viscous heating is determined by the disk's vertical optical depth (τ vert ) and gas surface density (Σ, Eq. 5,6) where κ = 0.5 cm 2 g −1 is the opacity, Ṁ = 10 where for a passive disk T accretion is set to 0. The midplane temperature profile is a key component of the model as all governing processes involve the temperature.In particular, the temperature can significantly alter ice line locations as viscous accretion heating pushes ice lines further out in the disk, especially for species with higher sublimation temperatures.

Ice lines
Ice lines are the radial and vertical locations in the disk where molecules freeze out onto silicate particles.We calculate the radial ice line locations for H 2 O, CO 2 , and CO, which are abundant volatiles in protoplanetary disks (see Oberg & Bergin 2021 for a review of disk chemistry).
An ice line is determined by the temperature (and to some extent disk surface density) at which a molecule's adsorption and desorption fluxes are balanced (see Hollenbach et al. 2008;Öberg et al. 2011;Powell et al. 2017;Oberg & Wordsworth 2019).This temperature is compared to the disk temperature profile in order to find the specific radius at which the fluxes are in steady state.The expressions for the fluxes (number of molecules per area per time) are: where n i is the molecular gas number density and c s = kT /µ is the isothermal sound speed.The number of adsorption sites per cm 2 on the particle is N s,i = 10 15 sites per cm 2 .The vibrational frequency of the molecules in the surface potential well is ν vib = 1.6×10 11 E i /µ i s −1 , with E i being the adsorption binding energy of the molecule in units of Kelvin and µ i being the molecular weight in grams.We take the fraction of occupied adsorption sites, f s,i , to be unity.The grain temperature, T grain , is assumed to be the same as the midplane temperature (see Eq. ( 11)).The adsorption and desorption fluxes depend strongly on the species molecular properties, which are listed in Table 2. Ice lines are key in determining which species are present throughout the disk.

RELATIVE VELOCITIES OF SOLID PARTICLES IN PROTOPLANETARY DISKS
The motion of small particles in a protoplanetary disk is primarily determined by the motion of the surrounding gas and by how strongly the particles move with the gas via gas drag.We calculate particle-particle relative velocities, following the works of Whipple (1972), Weidenschilling (1977b), Chiang & Youdin (2010), and Perets & Murray-Clay (2011).Particles throughout protoplanetary disks can be governed by three different drag regimes: Epstein, Stokes, and ram pressure.When the particle radius, s, is less than the mean free path, λ, of the gas so that s < (9/4)λ, the particle is in the Epstein drag regime.Typically, Epstein drag applies to particles in the outer disk that are well-coupled to the gas.For larger particles with s ≥ (9/4)λ, the governing drag regime must be determined using the Reynolds number, Re = 2sv rel /ν, where the kinematic viscosity ν = 0.5λv th , and vth = 8/πc s is the thermal velocity of the gas.When Re < 1 the particle is in the Stokes drag regime, while particles with Re > 800 are in the ram pressure drag regime.We approximate particles in the intermediate Reynolds number range by linearly extending the Stokes and ram pressure regimes and find the transition between the two to be where the two analytic coefficient of drag expressions intersect, Re ∼ 54.At this intersection, the approximation deviates from the intermediate regime by a factor of ∼ 4 (see Appendix A for details).
The drag forces in each regime are given by: where ρ g is the density of the gas and v is the relative velocity between the particle and the gas.
A particle is considered well-coupled to the surrounding gas if its orbital period is longer than the time it takes gas drag to stop the particle's motion.This is called the stopping time, and is found by dividing the particle momentum by the drag force.The dimensionless stopping time τ ≡ t s Ω is used throughout the model, and is also known as the Stokes number.Well-coupled particles have τ < 1.The density of the particle is ρ s = ϕ s ρ i , where ϕ s is the particle filling factor and ρ i is the material density (see Appendix B).
Marginally coupled particles, approximated by τ = 1, can be examined to build intuition for how the different drag regimes affect the sizes at which growing and drifting particles fragment due to collisions.These particles move at the highest velocities relative to the gas and are thus most likely to experience high energy particleparticle collisions that lead to fragmentation.The particle size at which τ = 1 as a function of orbital radius (Figure 5) can be found by solving for the sizes in Equation ( 15) and applying the drag conditions discussed above.Figure 5 compares the τ = 1 particles for each drag regime, displays the drag-dependant τ = 1 solution, and demonstrates how different particle compositions change the solution.A particle in the inner disk is dominated by the ram pressure regime, the central disk is dominated by the Stokes regime, and the outer disk is dominated by the Epstein regime.This simple analytic framework provides some intuition for where in the disk and at which particle sizes fragmentation may occur.

Particle-Particle Relative Velocities
Particle-particle relative velocities depend on particle size (s) and orbital radius (r) via the dimensionless stopping time (τ ).We apply analytic relative velocity expressions which are valid for all particle sizes and drag regimes throughout the disk.We numerically compute the stopping time for each particle through an iterative process which uses the particle's size and position in the Relative Velocity [cm s 1 ] Figure 6.Relative velocities are higher, and fragmentation is easier, for actively-heated disks with stronger turbulence.particle-particle relative velocities are displayed as a function of orbital radius and particle size for a passive (top) and active (bottom) disk.The columns correspond to varying turbulence parameterized by an α of 10 −2 , 10 −3 , and 10 −4 from left to right.Critical fragmentation velocities for BCCA and BPCA H2O solid aggregate particles are displayed as black labelled lines.
Varying profiles are shaped by transitions in drag regimes (see Figure 5).Particle size refers to the radius of the target particle, which in this case collides with a projectile half its size.
disk to determine the relevant drag force governing its motion.
There are two components to the relative velocity between two particles (v rel ): the relative laminar drift velocity and the relative velocity that arises due to turbulent motion of gas in the disk.The total relative velocity is the vector summation of these two components.
Laminar drift arises in disks as the gas orbits at sub-Keplerian velocities due to a radial pressure gradient.Solid particles which would otherwise orbit at the Keplerian velocity feel a headwind due to gas drag causing the particles to lose angular momentum and drift inward.The laminar drift velocity of an individual particle includes both radial (v r ) and azimuthal (v ϕ ) disk components which incorporate the inward drift and the orbital velocity (for a review, see Chiang & Youdin 2010): where η = 0.5c 2 s /v 2 K = 0.5H/r is the gas-pressure support parameter (H = c s /Ω is the disk gas scale height) and v K = GM * /r is the Keplerian velocity (Weidenschilling 1977b).With these components, the relative laminar drift velocity between two particles is then where the subscripts 1 and 2 refer to evaluation of Equations ( 17) and ( 18) for particle 1 and 2, respectively.We compute the relative velocity due to particle interactions with turbulent gas following the framework presented in Ormel & Cuzzi (2007), specifically Equations ( 16) and (21d) (see Appendix A of Powell et al. (2019) for details).When using Equation ( 16), we assume that the overturn time of the largest eddy is t L = 1 and separately calculate the various times based on the dimensionless stopping time of both colliding particles.For typical particle size distributions, collisions of comparably-sized particles both lead to the greatest rate of particle growth and the highest likelihood of collisional disruption.The gas velocity is taken to be v gas = √ αv th , where α is the Shakura & Sunyaev (1973) accretion disk turbulence parameter.Expected protoplanetary disk values for α lie in the range of 10 −5 to 10 −2 (Andrews 2020).Note that in our model, the α that parametrizes turbulent motion is a free parameter.In particular, given uncertainties in disk accretion models, we do not require α to have the same value that would be needed to model viscous accretion and we separately choose Ṁ in Equation ( 9).For reference, using the values given for TW Hya in Section 4 evaluated at r = 3AU, the accretion equation Ṁ = 3πΣν t , with turbulent viscosity ν t = αc s H, implies α = 3 × 10 −4 .The impact of turbulence on particle-particle relative velocity is displayed for passively and actively heated disks in Figure 6 with α's of 10 −2 , 10 −3 , and 10 −4 .Ultimately, the larger the turbulence, the larger the relative velocity, and with larger relative velocities fragmentation is more likely.For the remainder of this work, we choose a fiducial value of α = 10 −3 .
We note that particle settling may result in a minimum level of gas turbulence and hence a minimum expected α.Though we do not explicitly enforce this limit in our model, we comment on its magnitude here.Turbulence due to Kelvin-Helmholtz shear instability limits particle settling to the midplane of protoplanetary disks (e.g., Weidenschilling 1980;Sekiya 1998;Sekiya & Ishitsu 2000, 2001;Youdin & Shu 2002;Chiang 2008).The instability arises when the Richardson number (Chandrasekhar 1961) drops to a critical value of Ri ∼ 1 (Johansen et al. 2006a), corresponding to a maximum particle scale height (for a large midplane dust-togas ratio) of H p ≈ ηRi 1/2 r ∼ ηr (Chiang 2008;Gerbig et al. 2020).In our model, the particle scale height is calculated using (Ormel & Kobayashi 2012) Setting H p to its minimum allowed value yields α ∼ 0.5ητ /(1 − η).To drive the shear instability, particles must be sufficiently coupled to the gas to affect gas motion, so this effect generates maximum turbulence when particles have τ ∼ 1 and (since η ≪ 1), α ∼ 0.5η.As an example: at r = 3AU for the TW Hya parameters used above, η = 4 × 10 −4 for a passively heated disk and η = 10 −3 for an actively heated disk with Ṁ = 10 −8 M ⊙ /yr.Particles undergo collisional fragmentation if their relative velocity (v rel , derived in Section 5) in the disk reaches or exceeds the species-and structure-dependant critical fragmentation velocity (v crit , derived in Section 3).By defining fragmentation as where these two velocities are the same, we can solve for the sizes at which collisions will result in fragmentation at each orbital radius.In the inner disk, there exist two fragmentation sizes with v rel = v crit for each orbital radius, bounding a fragmentation region where v rel ≥ v crit centered on particles with τ = 1.Farther out in the disk, the range of particle sizes that fragment becomes smaller.For some material and disk properties, fragmentation ceases altogether in the outer disk.Using this framework, we can create regions of fragmentation in TW Hya for aggregate particles (based on Equation 4) composed of SiO 2 , Mg 2 SiO 4 , H 2 O, CO 2 , and CO, as well as for strong and weak non-aggregate solid particles (based on Equation 2), shown in Figures 7 and 8 respectively.The fragmentation regions are recreated for the MMSN in Appendix C for comparison.We note that even if particles do not reach critical fragmentation velocities, and are not affected by desorption and drift, they may still be prevented from growing by the erosion or bouncing barriers which are discussed in context of this work in Section 8.
While Figure 7 compares the fragmentation regions for a variety of compositions throughout the entire disk, particles composed of H 2 O, CO 2 , and CO ice are expected to desorb at their respective ice lines and are not expected to exist in the solid phase inward of their ice lines.The ice lines are represented as vertical lines Fragmentation regions in particle size-orbital radius phase space for SiO2 (yellow), Mg2SiO4 (green), H2O (blue), CO2 (red), and CO (purple) aggregate particles in TW Hya following the prescription from Wada et al. (2007Wada et al. ( , 2009) ) (Eq. 4).The left panels are for a passive disk and the right panels are for an active disk.The opaque region is for BPCA particles, while the fainter region is for BCCA particles.Ice line locations for H2O (blue), CO2 (red), and CO (purple) are shown as vertical lines.Regions for H2O (third row) directly correspond to the black outline highlighted in the middle column of Figure 6.  .with corresponding colors that match with composition color.The opaque colors represent BPCA particles and the more transparent colors represent BCCA particles.Silicate particles, as well as BCCA aggregates composed of CO 2 or CO ice, have the most widespread fragmentation regions although they are expected to be coated in ice beyond the H 2 O ice line.Nearly all BPCA ices only have fragmentation regions within their respective ice lines, meaning that all compact particles-and grains which are coated in these ices-do not undergo collisional fragmentation.Desorption, and potentially other barriers as discussed in Section 8, will ultimately determine whether these particles will grow to large sizes.All growing silicate particles in the inner disk, within the H 2 O ice line, are fragmenting.The silicates also match the expected fragmentation region for solid nonaggregate rocks using the critical disruption energy from Stewart & Leinhardt (2009) demonstrated in Figure 8.

FAVORABLE REGIONS FOR PLANETESIMAL GROWTH
We now focus on the interplay between collisional fragmentation and particle growth.An initially-small particle grows via coagulation and drifts toward the star due to gas drag.If growth and drift carry the particle into a fragmentation region, as illustrated in Figures 7  and 8, fragmentation occurs, preventing growth to larger sizes.If drift carries the particle past an ice line, desorption can cause it to lose the relevant ice species from its mantle and, potentially, to fall apart and begin growing anew in the absence of that solid species.We determine the fate of a particle by simultaneously integrating the equations for collisional growth, particle drift, and desorption, and then comparing our results with the frag-mentation regions computed in Section 6 and ice line locations computed as described in Section 4.2. Figure 9 illustrates an example outcome for a BCCA particle composed of H 2 O ice-including the τ = 1 profile (Figure 5), the relevant fragmentation region (passive BCCA H 2 O region from Figure 7), ice lines (derived in Section 4.2), and particle evolution paths (derived in Appendix D).The particle evolution paths are directly compared to regions of fragmentation, as well as ice lines, to check whether the particle will collisionally fragment or desorb during its evolution.
Analytically, a particle is expected to grow and drift along a path where the coagulational-growth and particle-drift timescales are equal (Equations D3 and D4), which is similarly demonstrated in Appendix 1 of Tsukamoto et al. (2017) where they find that the growth and drift timescales eventually converge.Numerically, we solve for the particle's mass and disk radial position over 1 Myr which is governed by the growth, drift, and desorption rates (dm/dt and dr/dt from Equation D5).In Figure 9, the analytic solution is shown as a dotted black line, while the solid black lines are the numerical solutions to the coupled differential equations.The three numerical solutions are initialized with disk radii of r = 20, 40, and 80 AU, with corresponding sizes of s = 1, 0.1, and 1 cm (the different particles are represented with a circle, triangle, and square respectively).
The general evolution of a particle in a disk includes the following phases: (1) the particle initially grows rapidly until the growth and drift timescales are comparable, (2) then the particle will grow and drift at the same rate until the particle's stopping time approaches unity (τ ∼ 1), and (3) once the particle reaches τ = 1 at a certain orbital radius, the growth rate increases while for a passively heated disk with f d = 10 −3 , α = 10 −3 , and ϕs = 0.3.Fragmentation extends inwards of the H2O ice line (barred blue), however the ice will sublimate and instead silicate fragmentation should be considered.The shape of the fragmentation region follows the H2O τ = 1 line from Figure 5 (blue line).The dashed black line is where the analytic growth and drift timescales (defined in Appendix D) are equal.After an initial phase of growth, particle evolution paths in the outer disk (solid black curves) follow the lower branch of this curve, offset by a multiplicative factor that arises from more accurate treatment of the integration.Each particle evolution (starting at the circle, square, and triangle) is shown for 1 Myr.Note that all three converge to the same evolution path where the growth rate is balanced by the drift rate.
the drift rate decreases causing the particle to grow in place.
Phase (3), where growth becomes much faster than drift, is likely for porous aggregates, however particles of little porosity where ϕ s → 1 may evolve differently.Okuzumi et al. (2012) find that porous aggregates are not hindered by radial drift in the inner disk and can continue to grow whereas compact solid particles (ϕ s = 1) are hindered, although they do not consider fragmentation.In all models, radial drift limits particle growth in the outer disk.We assume that particles are porous, keeping the filling factor the same to emphasize how particle structure and composition affect fragmentation.We use a value of ϕ s = 0.3, for both BPCA and BCCA aggregates, to represent particles in protoplanetary disks (Zhang et al. 2023).BCCA particles will typically have lower filling factors than 0.3 depending on the particle size and number of fractal dimensions as demonstrated in Tazaki et al. (2019).At lower particle filling factors, particle relative velocities are larger, causing fragmentation to be more likely (see Appendix B, Figure 13).During phase (2), as particles are growing towards τ ∼ 1, collisions will have likely compactified (or fragmented) BCCA aggregates, making BPCA aggregates more relevant for phase (3).
While phase (2) can be well described by the analytic solution, the numerical solution is needed to accurately model particle behavior once τ ∼ 1 in order to trace the desorption of an ice particle if it drifts interior to it's ice line.The particle will cease growing if it evolves into the corresponding species fragmentation region, but not necessarily if it drifts interior to its corresponding ice line.Since the particle may be composed of multiple species, the particle may fall apart at the ice line, however the smaller fragments will continue to grow and drift through the same phases described above.Because growth curves converge to a path along which the growth and drift timescales are comparable, the extent to which particles are disrupted at ice lines before reforming with a new mantle composition (for example the CO ice line shown in purple in Figure 9) does not substantially affect our results.In the H 2 O BCCA example of Figure 9, neither fragmentation nor desorption prohibit particle growth indicating that the outer disk may be a region where collisional growth is favorable.
As the particles drift inwards the overall dust-to-gas mass ratio, or the ratio between the solid and gaseous disk surface densities, will also evolve.At early disk lifetimes the dust-to-gas ratio is expected to be around 10 −2 corresponding to ISM values (Bohlin et al. 1978).At later disk lifetimes once the disk has had time to evolve, the dust-to-gas ratio can decrease by an orderof-magnitude to a value of around 10 −3 , or even lower (e.g., Birnstiel et al. 2010Birnstiel et al. , 2012;;Powell et al. 2019).For our modeling purposes we keep the dust-to-gas ratio constant throughout the disk and with time even though it varies for both (Alexander & Armitage 2007).The varying f d is relevant for the growth timescale (see Equation D3), thus we model particle evolution for both of these values to evaluate how composition-based fragmentation regions can change within a disk lifetime.
With all of these pieces, the final component of our model incorporates all of the species (SiO 2 , Mg 2 SiO 4 , H 2 O, CO 2 , and CO) into one fragmentation picture.In the inner disk, interior to the H 2 O ice line, the fragmentation regions which dominate are those for SiO 2 and Mg 2 SiO 4 .Beyond the H 2 O ice line, silicates are coated in ice (see Section 2), and the dominating fragmentation region is determined by the composition of the outermost layer of the ice mantle, monomer size, and how compact the aggregate is.Regions dominated by H 2 O, CO 2 , and CO are cut off inwards of their ice lines since these species are expected to desorb at the orbital radii of their ice lines leaving behind a silicate grain.Collisional fragmentation is not likely to occur for ice-coated BPCA particles as these fragmentation regions are mostly interior to their ice lines as illustrated in Figure 10, which also includes the numerical evolution paths for growing particles.BCCA particles, however, are likely to undergo fragmentation.While BCCA aggregate collisions may result in compaction at small sizes, as they get larger their growth will be limited by the meter-size barriers.
The results of Figure 10 demonstrate that there do exist regions of the disk where growth can be efficient and will not be impeded by fragmentation or desorption depending on the assumed midplane temperature.In general, for a passively-heated disk, growth can efficiently happen at early and later disk lifetimes-all three H 2 O particles evolve to sizes which are no longer limited by fragmentation.For an actively heated disk, which has a hotter midplane, at later times (f d = 10 −3 ), the H 2 O ice line is pushed further out in the disk such that particles will fragment before growing to large sizes.In either case however, BPCA particles starting beyond the CO 2 ice line may experience efficient runaway collisional growth beyond the fragmentation and drift barriers for f d = 10 −2 .Based on particle evolution paths in context of fragmentation and desorption, the favorable growth region in TW Hya is between the CO 2 and CO ice lines, while for the MMSN it is between the H 2 O and CO 2 ice lines.The differing regions between the MMSN and TW Hya are indicative that disk properties are an important factor in shaping favorable growth and fragmentation regions.

DISCUSSION
In this work, we develop a flexible framework for particle evolution that can be updated readily with improved constraints on fragmentation from future observations, laboratory experiments, and theoretical studies of velocity-dependant growth barriers.In this Section we discuss model sensitivities, other potential barriers to particle growth, observational tests and implications, and necessary constraints from laboratory data that could be used to validate or improve the reliability of this modeling framework.

Sensitivity to Monomer Size for Aggregate Fragmentation
As illustrated in Figure 2, the critical velocity for destruction of aggregate particles is strongly dependent on the size of their component monomers, r m .Our fiducial r m = 0.1µm is comparable to the sizes of ISM grains and is hence a minimum reasonable value (Oberg & Bergin 2021), yielding a maximum critical fragmentation velocity.In Figure 11, we recreate the top left panel of Figure 10 (TW Tya passive disk with BPCA and BCCA aggregate particles and α = 10 −4 ) for differing monomer sizes r m = 1, 10, and 100µm.Micron-sized monomers may narrowly find favorable growth regions within a disk in the case of BPCA aggregates, but for larger monomer sizes of 10 and 100µm no favorable regions exist.
As discussed in Section 2, monomer sizes are uncertain, and indeed may vary throughout the protoplanetary disk lifetime.In particular, very early in the disk lifetime, significant supersaturation of condensible gases likely allows the largest particles to accumulate substantial ice mantles via adsorption (Powell et al. 2022).Our choice of small monomers is most likely to apply after condensible volatiles are depleted by the growth and drift of these ice-rich particles.Because the surfaces of larger particles have less curvature, they nucleate condensation of volatiles more effectively (a process known as the Kelvin effect), meaning that small grains can remain ice-free even at early times (Powell et al. 2022).Hence, a time-evolving exploration of appropriate monomer sizes requires a complete size distribution of particles rather than the consideration of typical particle sizes applied in this work.Given our results, future work is merited to explore appropriate monomer sizes.This work will require self-consistent computation of the evolution of disk volatile abundances alongside particle coagulation, compactification, and volatile adsorption over the full distribution of particle sizes.
In the absence of such a full model, we appeal to Figure 1(b) of Powell et al. (2022) to demonstrate why we consider our fiducial choice to be plausible.The figure shows the abundances of ice-free and ice-coated particles at 30AU in TW Hya as a function of particle size at a late disk age of 5 Myr.Particles grow ice-free from s ∼0.1-100µm (first stage of Figure 1).At these small sizes, relative velocities are small (see Figure 6) and collisional growth likely leads to production of fluffy aggregates (e.g., Smirnov 1990;Meakin 1991;Kempf et al. 1999;Blum & Wurm 2000;Krause & Blum 2004;Wada et al. 2009;Paszun & Dominik 2006;Kataoka et al. 2013).As particles approach s ∼ 1mm, they become coated in ice and then drift inward.Particles that accreted ice early have already drifted away.The illustrated icy particles have adsorbed moderate ice coatings on top of aggregates having monomers of the initial small particle size, taken to be 0.1µm inherited from the ISM.These icy particles are consistent with case (a) of Figure 1.The black lines are numerically evolved particle evolution paths over 1 Myr for three H2O ice particles which begin at a certain size and orbital radius: circle (1 cm, 20 AU), triangle (0.1 cm, 40 AU), square (1 cm, 80 AU).The solid lines are calculated with a dust-to-gas ratio of f d = 10 −2 and the dotted-dashed lines are for f d = 10 −3 .The f d = 10 −3 case for the active disk fragments and desorbs before growing to large sizes.Fragmentation is pointed out as the red cross.The other particle evolution paths, however, do not undergo fragmentation or desorption and are able to grow to sizes no longer affected by the meter-size barriers.
Fragmentation is one of the most significant barriers to early stages of planet formation, however, there may be other barriers to continued small-particle growth, such as the bouncing and aeolian-erosion barriers.We note that the velocity of collisions where the bouncing barrier is relevant is not yet well constrained and may strongly depend on particle properties (e.g., Blum & Münch 1993;Langkowski et al. 2008;Windmark et al. 2012).Nietiadi et al. (2020) find that for small particles covered in ice mantles, the bouncing barrier critical veloci-ties range from 10 -100 m s −1 and are likely to increase for particles of larger size, meaning that the bouncing barrier limits collisional growth at higher relative velocities than the fragmentation velocities.Sirono & Ueno (2017) find that sintered BPCA H 2 O ice particles are likely growth-limited by the bouncing barrier rather than the fragmentation barrier as the monomer connections are stronger.They also find that sintered BCCA aggregates become more susceptible to fragmentation as compared to the unsintered case.Sintered particles are likely to exist within the disk environment, but for our purposes we only consider the unsintered case since the critical velocities from Sirono & Ueno (2017) fall within the critical relative velocity range in our study.Future work should include a distinction between the differing critical displacements for breaking monomer connections within a sintered versus unsintered aggregate particle.
The work presented here may be unaffected by sintering and the bouncing barrier as the ice particles that survive and undergo runaway growth never reach sufficiently high relative velocities where the bouncing or fragmentation barrier are efficient in limiting particle growth.Improved constraints on the bouncing barrier are necessary for a thorough understanding of particle evolution in disks.
The aeolian-erosion barrier (e.g., Paraskov et al. 2006;Schräpler & Blum 2011;Rozner et al. 2020) can efficiently erode larger pebbles and boulders, ranging from ∼ 10 − 1000 meters, down to centimeter sizes.Similar to collisional fragmentation, the aeolian-erosion barrier depends on a grain size (comparable to the monomer size in our discussion of collisional fragmentation), and the species's surface energy.Rozner et al. (2020) show that when grain sizes are 0.1 cm, the aeolian-erosion barrier can halt particle growth in disks at ∼10 cm.However, if we assume that the grain size is comparable to the nominal monomer size considered in this work (0.1 µm), the velocity threshold for the onset of aeolian-erosion increases by a factor of 100 such that aeolian-erosion will not significantly limit growth in protoplanetary disks.Future work may expand the model to include the aeolian-erosion barrier in order to evaluate the interplay between growth, drift, and aeolian-erosion for a variety of monomer sizes during early stages of planet formation.
In general, a complete dynamical model of particle evolution should include growth, drift, sublimation, and fragmentation, as well as other potential barriers to growth.While this is beyond the scope of our work, we highlight their importance for future studies.

Disk Substructure and Streaming Instability
Sustained particle growth in the favorable growth regions may potentially give rise to millimeter emission substructure.Rapidly growing particles in the diskdependant favorable growth region (between H 2 O and CO 2 ice lines for the MMSN, and between CO 2 and CO ice lines for TW Hya) can reach large sizes without drifting inwards, thus causing a pileup of particles resulting in a dearth of millimeter emission just interior to the favorable growth region.This may result in gaps in emission that would be present even in the case of young disks with larger dust-to-gas ratios (e.g., Williams & Cieza 2011;Andrews 2015).Several disks presented in Huang et al. ( 2018) from the DSHARP ALMA survey have significant substructure near their respective CO 2 ice line.Observations often look at older unobscured disks (f d = 10 −3 ), however, our results demonstrate that even at younger times particles can overcome the fragmentation barrier.Thus the work presented here may also be able to explain how significant millimeter substructure arises in young disks such as HL Tau (ALMA Partnership et al. 2015).Furthermore, If favorable growth regions are preventing particles from drifting interior to the H 2 O ice line, the abundance of gaseous H 2 O in the inner disk will be depleted.This H 2 O depletion is seen in several disk observations with the Herschel Space Observatory (e.g., Bergin et al. 2010Bergin et al. , 2013;;Kamp et al. 2013;Du et al. 2015;Salinas et al. 2016;Du et al. 2017), and may be explained with the favorable growth regions presented in our work.
Favorable growth regions are locations in the disk where particles can pileup as a result of their drift timescale becoming much longer than their growth timescale.We note that these regions will have higher dust-to-gas ratios, making them favorable for not only collisional growth but also for planetesimal formation through direct gravitational collapse.In particular, these pileups of growing particles could potentially instigate Resonant Drag Instabilities, particularly because particles in the favorable growth region naturally have τ ∼ 1, the approximate particle size most favorable to these instabilities and have higher densities of dusty and icy material.

Laboratory Constraints
Physical processes governing the growth and fragmentation of particles in disks strongly depend on laboratory derived material properties.In particular, a species surface energy, which is typically not well-constrained, determines the fragmentation properties of particles composed of that species.Improved surface energy constraints for all potential refractory and volatile species present in disks are thus necessary for an accurate picture of particle fragmentation.Furthermore, the surface processes during particle impacts are also not yet wellconstrained.In this work, we have assumed dry surface physics, however, collision-induced heating can potentially melt the ice on a particle's surface which can decrease the likelihood of particle fragmentation (Nietiadi et al. 2020).This may resemble sintered aggregates, although more work needs to be done in distinguishing the various surface connections, especially for particles coated in a variety or mixture of ices (Sirono & Ueno 2017).Disk particles with mixed compositions may also have collisional properties that vary from those presented in this work due to differing structures between the grain-ice and ice-ice layers (e.g., Fogarty et al. 2010).For example, when Musiolik et al. (2016b) tested mixtures of ice (CO 2 + H 2 O) the experiments found an order-of-magnitude increase in sticking velocity, indicating that ice mixtures can result in decreased particle fragmentation.Future work that describes the surface and collisional properties of mixed and partially-melted ice particles with a variety of disk-relevant compositions will be useful in improving models of particle evolution in disks.

SUMMARY & CONCLUSIONS
We develop a particle growth, drift, and fragmentation model which self-consistently solves for a particles' stopping time and includes ice sublimation, to test how different particle compositions and properties can withstand fragmentation.We find that ice particles in the outer disk, in between the CO 2 and CO ice lines for TW Hya and in between the H 2 O and CO 2 ice lines for the MMSN, may efficiently grow to large sizes without fragmenting or sublimating.
Our model produces regions of fragmentation in particle size-disk orbital radius phase space by comparing a particles' relative velocity with material and structural dependant critical fragmentation velocities.We test our results by varying the turbulence (α), dust-togas ratio (f d ), size ratio between colliding particles (χ), particle filling factor (ϕ s ), monomer size (r m ), disk temperature profile, and in particular particle composition.This study evaluates the fragmentation regions of the refractory species SiO 2 and Mg 2 SiO 4 , the volatile species H 2 O, CO 2 , and CO in their ice form, as well as a more general prescription for weak and strong non-aggregate rock.
For the fiducial case of TW Hya (α = 10 −3 , χ = 0.5, ϕ s = 0.3, r m = 0.1 µm) our model shows the following behavior: • In passive disks, all compact BPCA ice-coated H 2 O, CO 2 , and CO aggregate particles beyond ∼ 1AU are able to grow without fragmenting or desorbing.This is the case for both the MMSN and TW Hya, and is independent of disk age.
• In active disks, and at early times when f d = 10 −2 , compact BPCA ice-coated aggregate particles are also able to grow unimpeded beyond ∼ 10AU.
• Fluffier BCCA aggregate particles have fragmentation regions present throughout the disk making such particles more susceptible to fragmentation as compared to growth.BCCA H 2 O ice particles, however, are still able to grow beyond the meter-size in all scenarios except at late times (f d = 10 −3 ) in active disks.
• Uncoated silicate particles are still expected to collisionally fragment throughout the disk.In particular, SiO 2 using the adapted critical velocity from Wada et al. (2007Wada et al. ( , 2009) ) and weak rock using the adapted formalism from Stewart & Leinhardt (2009) will always have relative velocities reaching critical fragmentation velocities.Mg 2 SiO 4 and strong rock are less susceptible to fragmentation, and have fragmentation regions which reach to ∼1 AU and ∼10 AU for passive and active disks respectively.
• The critical fragmentation velocity decreases and particles become more susceptible to fragmentation as turbulence increases.H 2 O ice will fragment throughout the disk with α = 10 −2 .While for a passive TW Hya with α = 10 −4 , H 2 O ice will never reach collisional fragmentation velocities.
• Aggregates composed of larger monomers (r m > 10µm as opposed to 0.1µm) will likely undergo fragmentation.
These results indicate that particle growth may happen efficiently through collisions, without fragmenting or sublimating, beyond the H 2 O ice line for the MMSN and beyond the CO 2 ice line for TW Hya.D.P. acknowledges support from NASA (the National Aeronautics and Space Administration) through the NASA Hubble Fellowship grant HST-HF2-51490.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555.R.M.C. acknowledges support from NSF CAREER grant number AST-1555385 and support from NASA Interdisciplinary Consortia for Astrobiology Research (ICAR) grant 80NSSC21K0597.E.S.Y. would like to thank Karin Öberg for insightful discussions, and is grateful to the mentors, family, and friends who have supported and encouraged me throughout this work.We also thank the anonymous referee for their very helpful and thorough review.For our modeling purposes, we do not include the intermediate gas drag regime which spans Reynolds numbers 1 < Re < 800, in between the Stokes and ram pressure regimes.Instead, we extend the Stokes and ram pressure regimes into this Re phase space and use their intersection as the transition between the two.In general, the drag force is F D = 0.5C D πr 2 ρ g v 2 rel (Brown & Lawler 2003;Cheng 2009;Perets & Murray-Clay 2011), where C D is the coefficient of drag.For s > (9/4)λ, the Stokes and ram regimes have   .

Figure 3 .
Figure3.Different aggregate compositions define the overall strength of particles.Compact BPCA H2O ice is the strongest aggregate.Because strength is determined by surface interactions between monomers, silicates coated in ice will have strengths determined by that ice.

Figure 4 .
Figure4.For an active disk, the inner disk follows the viscous accretion heating profile (red) while the outer disk follows the passive stellar irradiation heating profile (orange).The temperature profiles for both the MMSN and TW Hya are displayed with dashed and solid lines respectively.

Figure 5 .
Figure 5. Left: Comparison of the Epstein, Stokes, and ram pressure drag regimes for τ = 1 particles.The dashed line represents the τ = 1 solution and serves as an intuitive proxy for the fragmentation region dependencies on drag regime.Right: The τ = 1 solution varies with composition due to different material densities assuming ϕs = 0.3 .
Figure7.Fragmentation regions in particle size-orbital radius phase space for SiO2 (yellow), Mg2SiO4 (green), H2O (blue), CO2 (red), and CO (purple) aggregate particles in TW Hya following the prescription fromWada et al. (2007Wada et al. ( , 2009) )  (Eq. 4).The left panels are for a passive disk and the right panels are for an active disk.The opaque region is for BPCA particles, while the fainter region is for BCCA particles.Ice line locations for H2O (blue), CO2 (red), and CO (purple) are shown as vertical lines.Regions for H2O (third row) directly correspond to the black outline highlighted in the middle column of Figure6. .

Figure 8 .
Figure8.Fragmentation regions calculated for solid rocks using the prescription fromStewart & Leinhardt (2009) (Eq.2), to be compared with the fragmentation regions from Figure7.Distinction between weak (light blue) and strong (dark blue) rock properties is discussed in Section 3. The left panels are for a passive disk and the right panels are for an active disk.These regions most closely resemble SiO2 and Mg2SiO4 aggregate particles from Figure7.

Figure 9 .
Figure9.H2O BCCA fragmentation region (shaded blue) for a passively heated disk with f d = 10 −3 , α = 10 −3 , and ϕs = 0.3.Fragmentation extends inwards of the H2O ice line (barred blue), however the ice will sublimate and instead silicate fragmentation should be considered.The shape of the fragmentation region follows the H2O τ = 1 line from Figure5(blue line).The dashed black line is where the analytic growth and drift timescales (defined in Appendix D) are equal.After an initial phase of growth, particle evolution paths in the outer disk (solid black curves) follow the lower branch of this curve, offset by a multiplicative factor that arises from more accurate treatment of the integration.Each particle evolution (starting at the circle, square, and triangle) is shown for 1 Myr.Note that all three converge to the same evolution path where the growth rate is balanced by the drift rate.

Figure 10 .
Figure 10.Composite fragmentation regions with desorption of ices taken into consideration for TW Hya (top) and the MMSN (bottom)-leaving SiO2, Mg2SiO4, BCCA H2O, BCCA CO2, and BCCA CO fragmentation.The left panels are for a passive disk and the right panels are for an active disk.

Figure 11 .
Figure11.Aggregates composed of large monomers (1, 10 and 100 µm) have larger regions of fragmentation that can span the entire disk over many sizes for all particle compositions.As in Figures7 and 10, the lighter and darker regions are BCCA and BPCA particles respectively, while colors correspond as yellow for SiO2, green for Mg2SiO4, blue for H2O, red for CO2, and purple for CO.These Figures are for a passively heated TW Hya with f d = 10 −3 , ϕs = 0.3, χ=0.5, and α = 10 −4 .

APPENDIXA.
DERIVING THE REYNOLDS NUMBER FOR DISTINGUISHING BETWEEN THE STOKES AND RAM PRESSURE DRAG REGIMES Figure 12 compares the C D in each regime as a function of Re.The intersection between Stokes and ram pressure occurs at Re ∼ 54.1, where C D is only smaller by a factor of ∼ 4.25 from the intermediate value.

Figure 12 .
Figure 12.Comparison of the analytic expressions of the coefficients of drag CD for the intermediate, Stokes, and ram pressure drag laws.The vertical dotted line is the switch between Stokes and ram pressure regimes, and is found to be at Re ∼ 54.1.

Table 2 .
Ice Line Molecular Properties

Table 3 .
(Stewart & Leinhardt 2009)terial Properties(Stewart & Leinhardt 2009) Wada et al. (2007ical velocity of aggregates for all compositions strongly depends on the monomer size.The SiO2 and Mg2SiO4 regions are the critical fragmentation velocities for silicate aggregate particles followingWada et al. (2007Wada et al. ( , 2009)  )with the lower bound and upper bound at each monomer size defined by BCCA and BPCA clusters respectively.The Strong & Weak Rock region is the critical fragmentation ve- Stewart & Leinhardt (2009)regate meter-sized rock following theStewart & Leinhardt (2009)energy calculations using values from Table3.
−8 M ⊙ yr −1 is the standard observed mass accretion rate, and Ω = GM * /r 3 is the orbital angular frequency (e.g., Table 4 includes a list of the fiducial parameters and their values used in calculating relative velocities and corresponding fragmentation regions in Figures 6, 7, 8, 9, and 10.See Appendix B for a discussion of the effects of varying size ratio (χ) and filling factor (ϕ s ) on the relative velocity.

Table 4 .
Fiducial model parameters