Global Modeling of Nebulae With Particle Growth, Drift, and Evaporation Fronts. II. The Influence of Porosity on Solids Evolution

Incremental particle growth in turbulent protoplanetary nebulae is limited by a combination of barriers that can slow or stall growth. Moreover, particles that grow massive enough to decouple from the gas are subject to inward radial drift which could lead to the depletion of most disk solids before planetesimals can form. Compact particle growth is probably not realistic. Rather, it is more likely that grains grow as fractal aggregates which may overcome this so-called radial drift barrier because they remain more coupled to the gas than compact particles of equal mass. We model fractal aggregate growth and compaction in a viscously evolving solar-like nebula for a range of turbulent intensities $\alpha_{\rm{t}} = 10^{-5}-10^{-2}$. We do find that radial drift is less influential for porous aggregates over much of their growth phase; however, outside the water snowline fractal aggregates can grow to much larger masses with larger Stokes numbers more quickly than compact particles, leading to rapid inward radial drift. As a result, disk solids outside the snowline out to $\sim 10-20$ AU are depleted earlier than in compact growth models, but outside $\sim 20$ AU material is retained much longer because aggregate Stokes numbers there remain lower initially. Nevertheless, we conclude even fractal models will lose most disk solids without the intervention of some leap-frog planetesimal forming mechanism such as the Streaming Instability (SI), though conditions for the SI are generally never satisfied, except for a brief period %for a brief stage around $\sim 0.2$ Myr at the snowline for $\alpha_{\rm{t}}=10^{-5}$.


INTRODUCTION
Observations over the last several decades indicate that not only does planet formation appear to be extremely robust, it can lead to very diverse outcomes (Lissauer et al. 2011;Winn & Fabrycky 2015;Triaud et al. 2017;Wittenmyer et al. 2020, and others). Yet, despite their ubiquity, a clear picture of how these systems came to be remains lacking because we still do not possess a coherent picture of how growth proceeds from sub-micron dust grains to planetary building blocks. Primary accretion of the first 100km size planetesimals, in particular, is not understood, but is a key initial condition for subsequent stages of planetary growth in our own solar system and beyond (e.g., Morbidelli et al. 2009;Kenyon & Bromley 2010;Schlichting et al. 2013;Chambers 2014;Johansen et al. 2014). The physics of secondary accretion -embryos that grow by pair-wise mergers of large bodies (e.g., Kary & Lissauer 1994;Kokubo & Ida 2002;Chambers et al. 2010;Bodenheimer et al. 2018), or by pebble accretion (e.g., Ormel & Klahr 2010;Lambrechts & Johansen 2012Bitsch et al. 2019) -is more well characterized, but the 100 − 200 km size pre-existing bodies adopted as the initial conditions for these models are merely assumed to be formed previously by some other means.
The single most critical factor in determining the path taken by primary accretion is the poorly known intensity of global turbulence. Grain growth begins by low-velocity "perfect" sticking of sub-micron monomers that are well-coupled to the gas, but relative velocities between growing particles in turbulence increase sufficiently that they become subject to a gauntlet of barriers that by themselves, or in combination, frustrate growth beyond a certain threshold. The most elementary barrier is simply fragmentation that can effectively grind particle growth to a halt. However, even before this, grains can encounter the so-called bouncing barrier where collisional energies are insufficient to lead to any fragmentation or erosion, but large enough to prevent sticking Zsom et al. 2010). The bouncing barrier does not halt growth entirely, but it can significantly slow growth by limiting the range of sizes that a grain can grow from (Windmark et al. 2012a;Estrada et al. 2016). Once particles grow to sufficient size that they begin to decouple from the gas and drift radially inward due to a headwind drag from the more slowly rotating, pressure-supported gas (Weidenschilling 1977;Cuzzi & Weidenschilling 2006;Johansen et al. 2014), they encounter a third, so-called radial drift barrier because they drift in faster than they can grow. This is especially problematic in turbulence, which decreases the local density of the midplane particle layer (Dubrulle et al. 1995). For compact particle growth models, radial drift tends to dominate outside the snow line (e.g., Brauer et al. 2008;Birnstiel et al. 2010;Estrada et al. 2016;Drażkowska et al. 2016). Even if incremental growth can somehow escape these barriers and continue to larger sizes, turbulence poses yet another barrier to growth when gas density fluctuations gravitationally excite disruptive collisions between 1 − 10 km size bodies (Ida et al. 2008;Nelson & Gressel 2010;Gressel et al. 2011;Ormel & Okuzumi 2013).
Indeed, after much debate over the years, it is increasingly thought that protoplanetary nebulae in the early stages of their evolution are at least weakly-to-moderately-turbulent in the regions of the disk ( 100 AU) in which particle growth is of the greatest interest (see e.g., Turner et al. 2014;Lyra & Umurhan 2019, for reviews). Even as Magneto-Hydrodynamic (MHD) models of turbulence have trended towards very weak turbulence even at high altitudes (Bai & Stone 2013;Turner et al. 2014;Gressel et al. 2015, and others), recent theoretical advances point to a number of purely hydrodynamical mechanisms (Nelson et al. 2013;Marcus et al. 2013Marcus et al. , 2015Lyra 2014;Stoll et al. 2017; Barranco et al. 2018, see Lyra & Umurhan 2019) that appear to be able to drive turbulence almost anywhere in the disk, at the level of α t ∼ 10 −5 − 10 −3 , where α t is the so-called turbulence coefficient which parameterizes the magnitude of turbulence strength.
This gauntlet of barriers to incremental growth in turbulence has led to the idea that 10 − 100 km bodies are assembled directly from a sufficiently dense reservoir of growth-frustrated smaller particles via some collective effect, such as the Streaming Instability (SI, Youdin & Goodman 2005;Youdin & Johansen 2007) or Turbulent Concentration (TC; Cuzzi et al. 2010;Hartlep & Cuzzi 2020) that under the right conditions can leapfrog all of these barriers. However, with few exceptions (e.g., Johansen et al. 2007, who assumed that meter-sized particles could exist in a fully MRIturbulent nebula), SI has only been modeled in globally nonturbulent disks and the viability of the SI under plausible turbulent conditions has not been seriously studied. Recently, Umurhan et al. (2020, see also Chen & Lin 2020 have established new criteria for the operation of the SI in arbitrary turbulence for arbitrary particle sizes. Comparing their results with sizes for compact particles determined in previous versions of this work (Estrada et al. 2016), they found that self-consistent conditions of particle size and turbulent intensity for planetesimal formation by SI are not generally obtained. Meanwhile, it has been found that TC may be able to produce planetesimals from swarms of growth-frustrated "pebbles" close to the realistic size range (Hartlep & Cuzzi 2020).
It is evident that, if the nebula is even weakly turbulent, accurate estimates of the limits on size and density of particles that can grow simply by sticking are needed to provide self-consistent initial conditions for these "leapfrog" models of planetesimal formation. Finding ways to overcome the radial drift barrier would be especially important: if particles could grow faster than they can drift, then it might be possible for them to grow to larger sizes without being lost to the inner disk regions, or even the central star. These particles would still have to contend with possible fragmentation, but if they could grow large enough and in sufficient number, and/or achieve some enhanced degree of local solids density, then they might more easily trigger planetesimal formation by SI or TC. To that end, in this paper we incorporate a model for porous particle growth and compaction in our global nebula evolution code (Estrada & Cuzzi 2008;Estrada et al. 2016) to explore the consequences of fractal aggregates on the solids evolution within the protoplanetary gas disk.
Similar models have been developed previously, but none have applied fractal growth and compaction to a model that incorporates all of the relevant physics and self-consistently accounts for the temporal and spatial evolution of the composition and size distribution of particles in a dynamically evolving disk as we do. Okuzumi et al. (2012) modeled the formation of planetesimals outside the snowline via fractal, porous particle growth and collisional compression in a static Minimum Mass Solar Nebula (MMSN, Hayashi 1981), but assumed perfect sticking and ignored fragmentation.
They found that although some local compaction occurs, the porosity of a growing aggregate continues to increase as a fractal, reaching planetesimal mass with internal densities of ∼ 10 −4 g cm −3 or less. In Section 2.4.5, we describe how we improve on these assumptions. Following along these lines, Krijt et al. (2015) used the Okuzumi et al. model in a static MMSN, but added fragmentation and erosion along with non-collisional compaction effects (Kataoka et al. 2013a). Using a Monte Carlo model, they show that erosion and restructuring by collisions with smaller particles stall growth well below planetesimal masses. Most recently, Homma & Nakamoto (2018) have applied the Okuzumi et al. model (and that of Kataoka et al. 2013a,b) in an evolving nebula with infall from the molecular cloud. Like Okuzumi et al., Homma & Nakamoto also assume perfect sticking and ignore fragmentation and erosion, and only consider the growth of aggregates outside the snowline. It is found that in this hotter nebula, where the snowline evolves outwards (∼ 7−15 AU) due to the infall, growth still stalls well below planetesimal masses even with perfect sticking. Drażkowska & Dullemond (2018) have also recently explored a model with infall in which their snowline evolves outwards with time during the buildup stage. They find that planetesimals mostly form later, but under certain conditions can form in the buildup stage, but these authors assume compact particle growth and adopt different turbulent strengths for disk evolution and turbulent mixing of the dust.
In this paper we employ the same fractal growth and compaction recipe (Suyama et al. 2012;Okuzumi et al. 2012;Kataoka et al. 2013a). However, none of these previous models consider multiple species or have a self-consistent calculation of disk opacity and temperature which depends on the evolving size-distribution (Estrada et al. 2016, hereafter, Paper I). Realistic disk composition will contain a mixture of silicates, iron metal, iron sulfide, organics, and water and other ices (see Table 1) each with its own associated evaporation front (EF). Our model is novel in the treatment of these, and is ideally suited to study the global implications of fractal growth over a range of (albeit imperfectly constrained) parameters. In Section 2 we summarize the basic workings of our global nebula evolution model, and expand on changes made to the code for this paper, including how we implement fractal growth and compaction into our numerical scheme. In Section 3 we describe the results of our simulations, comparing models in which we include fractal growth to their compact particle equivalent over a range of α t . We also present models that vary other disk parameters in Appendix B. Generally, while lacking infall, the nebula density and temperature, and the central star properties and timescales modeled here are characteristic of the Class 0/I stage of nebula evolution (< 1 Myr) during which time the first planetesimals form (Kruijer et al. 2017). In Section 4 we discuss implications for planetesimal formation, as well as observations of particle porosity. In our companion paper, we discuss the compositional evolution of our disk models in more detail (Estrada & Cuzzi 2022, hereafter Paper III). In Section 5 we summarize our results.

NEBULA MODEL
The simulations presented herein use our parallelized 1 + 1D radial nebula code in which we simultaneously treat particle growth and radial migration of solids, while evolving the dynamical and thermal evolution of the protoplanetary gas disk. Our code includes the self-consistent growth and radial drift of particles of all sizes, accounts for vertical settling and diffusion of the smaller grains, radial diffusion and advection of solid and vapor phases of multiple species of refractories and volatiles, and contains a self-consistent calculation of the opacity and disk temperature that allows us to track the evaporation and condensation of all species as they are transported throughout the disk. In the following sections, we only briefly summarize our code and will mostly focus on the description of additional physics or improvements that have been implemented for this work.
One difference between our implementation here and in Paper I is that in our previous work we employed an asynchronous time-stepping scheme. Briefly, each radial bin had its own timestep, with the innermost radial bin of the grid defining the minimum particle growth timestep ∆t min , which is a fraction of the local orbital period. This allowed for significant computational savings because only the innermost bin is called every ∆t min , while other bins were called at the appropriate time interval. We utilized a much larger global timestep, typically ∆t global = 100×∆t min to "synchronize" all radial bins to the current absolute simulation time, and then executed spatially global calculations like solving the gas and particle evolution equations, and determining the disk temperature. In this paper, we now use ∆t = ∆t min for all radial bins which leads to higher accuracy and faster run times, but at the expense of increased computational resources. However, we do continue to use a global timestep because as we showed previously (see Appendix C of Paper I) these global processes are not changing quickly over ∆t global compared to growth times, even for the highest turbulence strength we employ here. The reader is referred to Paper I for a more detailed description of our code.

Gas Disk Evolution
The initial conditions for the gas disk models used in this work are derived from the analytical expressions of Lynden-Bell & Pringle (1974), as generalized by Hartmann et al. (1998) in terms of an initial disk mass M disk and radial scale factor R 0 (see Paper I). This gas surface density distribution is fairly similar to that of Desch (2007) in the 0.3 − 0.5 Myr timeframe of both models, and represents a denser inner nebula than the MMSN. The time-dependent evolution of the gas surface density Σ and gas radial velocity v g in one dimension are then obtained by solving (Pringle 1981) The turbulent eddy kinematic viscosity ν = α t cH is parametrized in terms of a turbulent intensity α t , and depends on the evolving disk temperature T both through the gas sound speed c, and the nebula gas density scale height H. The gas mass accretion rate is given byṀ = −2πRΣv g whereṀ > 0 indicates flow towards the star, andṀ < 0 indicates mass flux outwards. The turnaround radius, at which the mass flux changes sign, is generally an increasing function of time, and at t = 0 is R t R 0 /2 for our initial conditions (Hartmann et al. 1998). The water snowline begins around ∼ 10 − 15 AU in our models, so they can perhaps best be associated with the post buildup stage due to infall (Drażkowska & Dullemond 2018;Homma & Nakamoto 2018).

Evolution of Solid and Vapor Species
Our model is capable of tracking multiple species in both vapor and solid phases. In Table 1, we list the condensibles used in this work along with their corresponding condensation temperatures T i , compact particle density ρ i and initial global mass fraction in the solid statex i relative to H and He. The initial mass fractionsx i are constructed using data from Table 2 of Lodders (2003) but differ in detail from Lodders' values. To get our values forx i , we started with the initial "refractory organics" or "CHONs" fraction from Pollack et al. (1994) as a guide for C (see also Jessberger et al. 1988;Jessberger & Kissel 1991;Lawler & Brownlee 1992;Mumma & Charnley 2011), and then partitioned the remaining C in the ratio of 1:1:1 between the "supervolatiles" CO, CH 4 and CO 2 . The silicates are determined by placing all the Mg into a combination of orthopyroxene (MgSiO 3 ) and olivine (Mg 2 SiO 4 ), while the refractory iron (metal) fraction is what remains after all of the S is placed in FeS. Some of the remaining O is used in determining the Ca-Al "refractory" fraction which is a mixture of Al 2 O 3 , CaO and CaTiO 3 , while the rest is placed in water. By contrast, Lodders (2003) has no refractory organics, assumes no O-bearing supervolatiles, allows water to capture almost all the O that is not consumed by silicates, and assumes the C is in methane ice or in methane hydrate that travels with the water ice. The total metallicity for our disk models isZ = ix i 0.014. We subdivide the particles into "dust" and larger "migrator" phases -essentially particles smaller than or larger than the fragmentation mass (see Sec. 2.4,and Paper I). These and the vapor phase fraction of each species are denoted by subscripts d, m and v respectively. As the disk evolves, these fractions change with time as particles grow and are transported throughout the nebula. The local instantaneous mass fraction, or concentration x i of each phase of species i, which is the ratio of the surface density of the species to the gas surface density, varies with time and is given by The total phase fraction of all species f v,d,m is just the sum over i in Eq. (3), e.g. for the smaller "dust" fraction, Similarly, the locally varying total mass fraction of condensables is the sum over all i, . The evolution of the gas surface density (Sec. 2.1) includes the vapor fraction of all species, but we do not include the contribution of the vapor phase in the molecular weight of the gas which would affect the sound speed, and thus viscosity (see, e.g., Schoonenberg & Ormel 2017;Charnoz et al. 2021). Large enhancements in the vapor phase could lead to a significant decrease in c at EFs (see Paper III) which then might lead to a pressure bump that can produce even stronger enhancements in solids, and expand the conditions under which, for instance, leap-frog planetesimal formation mechanisms may be satisfied (see Sec. 4.1 for more discussion). We will revisit the influence of the gas molecular weight in a future paper.

Particle Stopping Times
Particle growth in the protoplanetary disk begins by sticking of sub-micron grains which are dynamically coupled to the nebula gas. As grains grow into larger particles and agglomerates, they begin to decouple from the gas and collide at higher velocities. This influence of the nebula gas on the motion of particles is determined by the Stokes number where t ed is some characteristic eddy turnover time which we take to be the integral scale, or the turnover time of the largest eddy Ω −1 for global turbulence (e.g., Cuzzi et al. 2001;Carballido et al. 2010). Here, Ω = GM /R 3 is the orbital frequency at semi-major axis R. The stopping time t s is the time needed for gas drag to dissipate the momentum of a particle of mass m relative to a gas with local volume density ρ. The drag force a particle of radius r feels depends on its size relative to the molecular mean-free path λ mfp : where λ mfp = µ m /ρc, with µ m = 1.3 × 10 −4 g cm −1 s −1 is the dynamic viscosity of the gas, and A is the projected area of an aggregate. Equation (5) describes two distinct flow regimes, the Epstein flow regime and the Stokes regime. Particles can be well or poorly coupled to the gas in either regime, in principle. In the Stokes regime, particles' stopping times are affected by a drag coefficient C d which depends on the particle-to-gas relative velocity ∆v pg through the particle Reynolds number Re p = 2rρ∆v pg /µ m . For the simulations in this paper, Stokes flow particles typically have Re p 1 so that C d = 24/Re p (Weidenschilling 1977).

Radial Diffusion-advection and Vertical Diffusion
The radial motions of both solid and vapor species are determined for each species and phase from the advectiondiffusion equation (Desch et al. 2017) wherev is the net, inertial space velocity (advection + drift), D is the diffusivity and S i represents sources and sinks for solids and vapor of species i which include growth, radial transport and destruction of migrating material (see Sec. 2.4.3,and Paper I). For the vapor phase,v = v g and D = ν, whereas for the dust particles,v and D (which are both functions of height z) are determined at any radius R in the midplane layer or "subdisk" (see below) through a mass volume density-weighted mean over particle masses less than or equal to the fragmentation mass where ρ d = k ρ d k , and with h d k the particle scale height as defined in Equation 10 below. The radial drift of particles relative to the gas due to the local pressure gradient (Nakagawa et al. 1986;Takeuchi & Lin 2002) v rad has two components. The first is imposed by the radial motion of the gas moving with advective velocity v g , while the second is the radial velocity of the particle with respect to the gas. It is important to note that the sign of the radial drift velocity depends on particle mass -massive particles tend to drift radially inwards (v < 0), while less massive ones can be radially advected outwards with the gas (v > 0). We account for this in our code by determining the particle mass (wherev ≈ 0) that separates the population of particles moving inward from that moving outward, and then using the respective fractional masses of these populations as a weighting factor in solving the advection-diffusion equation for each population (see Paper I). For particle masses beyond the fragmentation mass, i.e., migrators or "lucky particles", we explicitly follow the evolution (including growth and radial transport) of their mass distribution in our code (see Sec. 2.4, and Paper I).
The vertical distribution of the mass volume density of particles ρ d k (z) and ρ m l (z) are accounted for through a balance between settling and diffusion that combines elements of several studies (Dubrulle et al. 1995;Cuzzi & Weidenschilling 2006;Youdin & Lithwick 2007). The scale height for particles of mass m in Equation 8 is given by (see Appendix B, Paper I) and similarly for migrators. The particle mass whose scale height sets the subdisk height, which defines the midplane volume where particle growth is followed in our code (Sec. 2.4), is taken to be one half of the mass of the mass-dominant particle at any given R (see Paper I). Finally, the total mass volume density at R of solid material in the subdisk is the sum over all particle mass bins in the dust and migrator populations indexed by k and l, ρ solids = k ρ d k + l ρ m l .

Disk Thermal Evolution
Because the disk temperature depends on the evolving particle size distribution through the opacity, a self-consistent calculation is essential to capturing the disk's dynamical evolution. We assume that the nebula is heated through a combination of internal viscous dissipation (primarily near the midplane), and external illumination by the stellar luminosity L . The equation we use for determining the midplane temperature is given by Nakamoto & Nakagawa (1994): which must be solved iteratively because the Rosseland and Planck mean optical depths τ R and τ P , respectively, are temperature dependent and thus affect the evolving particle size distribution and the solids and vapor fractions of all species. In Eq. (11), the first term on the RHS is a local, vertically integrated viscous dissipation rate, while the second term on the RHS accounts for the stellar flux on each disk face, and φ is a grazing incidence angle depending on disk geometry. In our model, we assume a flared disk geometry with a general radial variation given by Chiang & Goldreich (1997, see also Kenyon & Hartmann (1987; Ruden & Pollack (1991)) where R AU is radial distance in AU. As we did in Paper I, we employ a time variable luminosity using the model for a 1 M star (D'Antona & Mazzitelli 1994;Siess et al. 2000). For the latter, the initial luminosity L is roughly 12 L at the beginning of a simulation, which drops to ≈ 3 L after 0.5 Myr. One issue identified with these models is that they are based on oversimplified initial conditions that introduce uncertainty in the evolutionary track for stellar ages t 1 Myr (Baraffe et al. 1998(Baraffe et al. , 2002. However, the luminosity reduction factor may be at most ∼ 10% (see Palla & Stahler 1999) which we consider reasonable for the purposes of our simulations. In future work, especially when considering different metallicities as well as stellar masses, we will consider more recent evolutionary track models (e.g., di Criscienzo et al. 2009;Tognelli et al. 2011).
The optical depths τ in Eq. (11) are functions of the opacity κ as τ = κΣ/2. We define the Rosseland and Planck mean opacities from the basic wavelength-dependent opacity κ λ weighting in the standard way To determine the κ λ , we utilize the opacity model of Cuzzi et al. (2014, see also Sengupta et al. 2019) which includes realistic material refractive indices for the species listed in Table 1, and particle porosity. We note that we currently use the optical constants for (crystalline) silicates from Pollack et al. (1994). However, more recent applications (e.g., Birnstiel et al. 2018) tend to use (amorphous) astronomical silicates (Draine 2003) which are more absorbing, especially at shorter wavelengths by up to an order of magnitude. We will consider these in future applications since several recently discovered hydrodynamical instabilities which can lead to sustained turbulence (Lyra & Umurhan 2019) are operationally sensitive to the magnitude of the disk solids opacity. We have not included the Rosseland and Planck gas opacities in our models as they really only start to become important for temperatures above 2000 K which are not achieved in our simulations, though we have included them in a followup paper (Sengupta et al. 2022, which also includes a model for disk winds) using a table lookup derived from the work of Freedman et al. (2014). In our model, we specifically treat the evaporation fronts (EFs) -those locations in the disk where phase changes between solids and vapor can occur -of all species. The evaporation of inwardly migrating material, and the subsequent recondensation of outwardly diffused vapor onto grains, can lead to significant enhancements in solids (outside) and vapor phase (inside) an EF, as well as altering the composition of solids and vapors with implications for accretion, chemistry and mineralogy. We do not treat EFs as sharp boundaries, but allow for the "phase" change to occur linearly over a small temperature range T i − ∆T EF ≤ T ≤ T i + ∆T EF , with ∆T EF = 0.05 K. Allowing for this gradual radial transition prevents unrealistic drops in opacity and temperature just inside an EF, and effectively mimics buffered temperature changes as material is evaporated or condensed over a range of nebula altitudes (see Paper I, for more discussion). This approach works well for compact particles, but it should be acknowledged that for large particles, and especially for the large fluffy aggregates in our fractal growth simulations, it may not capture a scenario in which an aggregate's surface layers insulate volatile material in the interior. An improved model would require we treat the kinetics of evaporation and condensation at EFs (e.g., see Ros & Johansen 2013;Schoonenberg & Ormel 2017) which we will implement in future work.
For simplicity, we have also assumed that evaporation and condensation are reversible processes, but this is almost certainly not the case for the refractory organics. We would expect that these CHONs would likely break down irreversibly into CO, CO 2 , NH 3 and "carbon chains" like acetylene. These products would have an effect on the magnitude of enhancements in both refractory carbon solids outside, and their vapor fraction inside the organics EF (see further discussion in Paper III). The depletion of C in this way would be more consistent with the observation that the most primitive carbonaceous chondrites contain only ∼ 10% of the carbon found in comets (see Woodward et al. 2021). On the other hand, organics may be stickier than what we assume in this work (Kouchi et al. 2002;Homma et al. 2019, though see Bischoff et al. 2020) which might facilitate faster growth to larger sizes, and/or also allow for retention of some organics sequestered in the interiors of the larger particles and aggregates.

Particle Growth
We handle particle growth using the moments method solution to the coagulation equation for the "dust" portion f (m) of the particle size distribution (Estrada & Cuzzi 2008) smaller than the fragmentation barrier mass m f (Sec. 2.4.2), and explicitly follow growth of "migrators" for masses larger than m f . The moments method greatly speeds up calculations and allows us to focus on the more interesting evolution of the larger sizes. In the following sections, we summarize the collisional physics used in our code, and describe how we add fractal, porous growth to our model.

Relative Velocities
We include both stochastic and deterministic relative velocities in our code. The former are due to either collisions between grains and gas molecules, or interaction of particles with turbulent eddies within the gas, while the latter are radial and azimuthal drift, and vertical settling, that result from the friction caused by particle-gas interactions within the mean gas flow (Nakagawa et al. 1986;Weidenschilling & Cuzzi 1993, and section 2.2.1).
For very small grains, thermal or Brownian motion dominates relative velocities between particles of masses m and m. The mean relative velocity is where T is the local temperature and k B is the Boltzmann constant. Generally these stochastic motions are only important for the smallest particles in our distribution, and at the highest temperatures found in the inner disk regions. For the turbulence-induced relative velocities between particles of arbitrary size, we employ the closed-form expressions of Ormel & Cuzzi (2007) in which ∆v I and ∆v II are velocity contributions from so-called Class I and Class II eddies (see Völk et al. 1980). In terms of the Stokes numbers of m and m , and the turbulent large eddy gas velocity v t = α 1/2 t c, where Re = ρν/µ m is the gas Reynolds number. Equations (16) and (17) illustrate the different coupling that exists between particles and eddies of different sizes, and St * 12 = MAX(St * , St * ) is the larger of the two particle Stokes numbers at the boundary between these two eddy classes. We solve for St * 12 using Eq. (22) of Ormel & Cuzzi (2007) which takes into account eddy-crossing effects. When these effects become important, it can influence the Stokes number at which fragmentation occurs (see Sec. 2.4.2).
We note that for very small particles, when St < Re −1/2 ≡ St η , ∆v t = ∆v I = 0 for same-sized particles. Here St η = t η /t ed identifies the Kolmogorov scale, where t η is the turnover time of the smallest eddy. However, as has been pointed out (e.g., Okuzumi et al. 2011), the dispersion in the mass-to-area ratio m/A between aggregate to aggregate gives rise to a small relative velocity between them. This turns out to be particularly important in the treatment of fractal aggregates with fractal dimension D 2. To account for this, we replace the area A in the differential mass-toarea ratio ∆(m/A) between aggregates with the mean projected areaĀ of the aggregates (as in Okuzumi et al. 2009, their Eq. 47) and the size of its standard deviation normalized by the mean mass-to-area-ratio (e.g., Okuzumi et al. 2012): where we adopt δ = 0.1 as in Okuzumi et al. (2011). For larger particles, systematic pressure-induced velocities arise from drag forces (primarily azimuthal) from the nebula gas, which rotates more slowly due to the radial gas pressure gradient. We solve for the radial and azimuthal components of the velocity for the particles and the gas using a set of equations generalized from Nakagawa et al. (1986) for a particle size distribution using matrix inversion (see appendix A.4, Paper I). The mean particle-to-gas deterministic relative velocity can be constructed in a similar fashion. In the case where the local solids fraction is high 1 Ishihara et al. (2018) recently confirmed a finding by Pan & Padoan (2015), but at much higher Reynolds number, that the Ormel & Cuzzi model predictions for relative interparticle velocities seem to be a factor of 1.9 too large for particles with St < 1. Pan & Padoan (2015) made some suggestions about the reasons for this, which are connected to simplifying assumptions by Ormel & Cuzzi (2007) and earlier works on the subject, and which we have assessed in a preliminary way and found plausible. Fortunately it is easy to correct the Ormel & Cuzzi model expressions to agree with actual numerical simulations. Moreover, Pan & Padoan (2015) also note that the actual property of interest, the collisional velocity that most people use the relative velocities to represent, is actually slightly larger, decreasing the discrepancy of the Ormel & Cuzzi expressions to a factor of 1.3 or less. Consequently here we use the expressions in Ormel & Cuzzi (2007). enough to produce a feedback effect on the gas, the relative velocities must be determined iteratively. The normalized gas pressure gradient 2 where v K = ΩR is the local Kepler velocity, can be quite strong at the outer edge of the disk as ρ can decrease quite sharply there. This can lead to rapid inward radial drift of even very small particles. In most cases when the local solids-to-gas mass volume density ratio = ρ solids /ρ is not too large, the particle drift velocity can be well approximated by (Weidenschilling 1977;Cuzzi & Weidenschilling 2006) We make use of this approximation in Sec. 2.4.3. The rms relative systematic velocity between m and m due to this drag, as well as differential settling, is with the vertical settling velocity w = −Ω 2 zt s = −Ωz · St. A detailed description of how we calculate these are given in Appendix A4 of Paper I. We then obtain the traditional mean particle-to-particle collisional velocities by combining all of the relative velocity components in quadrature (e.g., Tanaka et al. 2005;Okuzumi et al. 2012 Our models employ a Gaussian relative velocity PDF centered on these mean values (Sec. 2.4.3). This was acknowledged to be a simplification in Paper I, since it is not strictly correct to include the fluctuating turbulent relative velocity with the deterministic velocities in this way. Garaud et al. (2013)  Small grains are efficient at sticking at smaller sizes owing to their low relative velocities. But as they grow, their relative velocities increase, gradually decreasing their sticking efficiency S and eventually halting further growth altogether. We account for these barriers in our collisional kernel using mass-and-velocity-dependent sticking coefficients for m ≤ m. One such barrier to growth derived from experimental work is the so-called bouncing barrier Zsom et al. 2010) in which the relative velocity between two grains is too high for any sticking, but too low for fragmentation. If bouncing is included in the kernel, then gives the first barrier encountered during growth as S b → 0 where the threshold velocity v b for bouncing collisions for either compact or porous aggregates is mass-dependent. For instance, for similar-sized compact silicate particles, Güttler et al. (2010) found v b (C 0 /m) 1/2 , with C 0 = 10 −7 g cm 2 s −2 . For compact grains with density ρ p , and when turbulence dominates their motion, the Stokes number St b at which bouncing of equal sized particles occurs under typical nebula conditions is approximately where S b = 0 defines the bouncing mass m b that corresponds to St b . For fractal, porous grains (Sec. 2.4.4), we adopt a different form for v b with a knee derived from Güttler et al. (2010) for similar-sized aggregates: with C 1 = 1.55 × 10 −7 g cm 2 s −2 , and C 2 = 8 × 10 −8 g cm 4 s −4 .The bouncing barrier does not halt growth, but merely slows it, by restricting the size range from which a given aggregate can grow to particles smaller than m b . Collisions become fragmenting, however, when the relative velocities between two particles of similar size exceed a certain threshold energy Q f (Stewart & Leinhardt 2009;Beitz et al. 2012). In a similar way to bouncing, the fragmentation coefficient determines the fragmentation mass m f when S f → 0. Unless turbulence is very weak, the highest relative velocities tend to be between similar-sized particles, so that fragmentation will first occur there. The fragmentation Stokes number in turbulence can be similarly derived from Eq. (17) when St = St and dropping terms in Re −1/2 which are negligible in the "fully intermediate regime" (Ormel & Cuzzi 2007): We include the factor y * = St * /St here explicitly to distinguish between cases where eddy-crossing effects are important. When St 1, y * ≈ 1.6 so that St f takes on the more familiar form of the last term on the RHS of Eq. (27). However, eddy-crossing effects can become important when St/α t η −1 (Ormel & Cuzzi 2007, also see Jacquet et al. 2012) so that y * can be greater than 1.6, and thus decrease the size at which fragmentation occurs. On the other hand, when St ∼ Re −1/2 , Eq. (27) can underestimate the fragmentation size because terms in Re −1/2 cannot be ignored. Our models cover the full range, and both situations can and do occur.
It can be seen from Eq. (27) that in regions of the disk where the temperature does not vary significantly (i.e., in the outer disk) St f is nearly constant, though the value of Q f may vary with composition which we allow for in our models. For silicate particles we take Q f = 10 4 cm 2 s −2 (equivalent to a threshold velocity of 1 m s −1 , e.g., Blum & Münch 1993;Stewart & Leinhardt 2009), and Q f = 10 6 cm 2 s −2 for water ice grains which are apparently stronger, and may grow larger and more rapidly than silicate grains (Wada et al. 2009(Wada et al. , 2013Okuzumi et al. 2012), though perhaps only over a limited range of temperatures near the evaporation point (Musiolik & Wurm 2019). For the more volatile ices (and for water-ice in two of our models which we refer to as "cold H 2 O ice", when T 140 K, see Table 2 and Appendix B), we adopt Q f to be the same as for silicate particles (Musiolik et al. 2016). Given that our aggregate particles are composed of different fractions of icy and non-icy species, we do a mass fraction-weighted to determine the fragmentation energy at any disk location. A similar approach is taken for the bouncing threshold (Paper I).
Finally, the outcomes of high-velocity collisions have been observed experimentally to lead to growth (Wurm et al. 2005;Kothe et al. 2010) when a mutual collision between different-sized particles can erode but not fragment the larger particle, and is enough to fragment the smaller of the two. In such a circumstance, it is possible for there to be deposition of mass on the larger particle if the efficiency of accretion exceeds that of erosion. This is the so-called mass transfer regime. When one considers a PDF of velocities, as we do for growth beyond the fragmentation barrier (Sec. 2.4.3), mass transfer can allow so-called "lucky particles" to grow much larger than the mass-dominant particle size of the distribution. We use the model of Windmark et al. (2012b,c, also see Beitz et al. 2012 to account for mass transfer, and the description of our implementation is fully described in Paper I.

Destruction Probabilities Beyond the Fragmentation Barrier
For particle growth beyond the fragmentation mass m f , we employ a statistical scheme using a PDF of relative velocities to determine the probability that a migrator particle (i.e. m > m f ) is destroyed due to fragmenting collisions with equal-size and smaller particles in the distribution. For compact particles, the probability that a migrator of mass m will be fragmented during a time ∆t by some m ≤ m is given by where ζ(m , m) is an integral over the PDF of relative velocities. In this work, as we did in Paper I, we primarily use a Gaussian PDF of relative velocities with rms velocity equal to the mean (Carballido et al. 2010;Hubbard 2012;Pan & Padoan 2013) where the integration is over the range equal to, or greater than the critical impact velocity 3 v c = Q f (m + m)/m . In Eq. (29, and also C1 in the Appendix C), v is used as an integration variable and should not be confused with the azimuthal velocity component of a particle as defined in Sec. 2.4.1. Although other workers have utilized similar Gaussian PDFs, it has been argued that a Maxwellian distribution may be more appropriate, especially when relative velocities are dominated by Brownian or turbulent motion (Galvagni et al. 2011;Windmark et al. 2012a;Garaud et al. 2013). To that end, we also explore models that utilize a Maxwellian PDF (Garaud et al. 2013) in Appendix C in order to compare with our usage of a Gaussian PDF.

Fractal Growth
Initially, our models begin with a "dust" powerlaw size distribution in mass such that particle number density f (m) ∝ m −q where q = 11/6, composed of compact, spherical grains with a minimum size of r 0 = 0.1 µm, and thus a monomer mass is m 0 = (4/3)πρ p r 3 0 where we define ρ p as the compacted density, or the density of a compact particle 4 . In this work, an initial monomer density that contains all of the solid phase species in Table 1 is ρ p = 1.52 g cm −3 (though note ρ p varies both temporally and spatially), which yields a monomer mass m 0 = 6.4 × 10 −15 g. As particles grow, they remain spherical aggregates regardless of whether growth is compact, or fractal -we cannot make a distinction between them based on shape 5 . As demonstrated in Estrada & Cuzzi (2008) with regards to the moments method we utilize for the particle distribution with m < m f , fractal growth can be incorporated into the collisional kernel without any loss of generality. Porous particles of mass m and radius r can be defined in terms of a fractal dimension 2 ≤ D ≤ 3 as m = m 0 (r/r 0 ) D and r = r 0 (m/m 0 ) 1/D . The density, or similarly the volume, of a fractal, porous particle can also be defined in terms of a enlargement factor or inverse filling factor ψ (e.g., Krijt et al. 2015) as where V comp is the volume the compact (D = 3) particle of equivalent mass would occupy, and ρ int is the internal density of a porous aggregate. A familiar description of a porous aggregate is its porosity ϕ which is 1 minus the volume filling fraction: ϕ = 1 − ρ int /ρ p . Thus, ψ = 1/(1 − ϕ), and the relationship between the fractal dimension and enlargement factor is The stopping time in the Epstein regime (Eq. 5) depends on the midplane gas density ρ = Σ/ √ 2πH and can be rewritten in terms of the Stokes number and surface mass density (Estrada & Cuzzi 2008) 6 as: It is interesting to note that D = 2 represents a special case that an aggregate of any size and mass has the same Stokes number as a monomer, St = √ 2πρ p r 0 /Σ, and would remain strongly coupled to the gas. However, aggregates will eventually suffer compaction through collisions, and their St will grow. The main effect of porosity is that porous 3 There is a typo of a factor of 2 in the definition of the critical impact velocity in Paper I. 4 That is, if the particle internal density is constant, the number of particles dN (m, dm) lying in a mass bin between (m, m + dm) is f (m)dm = N 1 m −11/6 dm and expressed in terms of particle radius is dN (r, dr) = N 2 r −7/2 dr. 5 In our models, we assume that the lower bound of our particle-size distribution r 0 remains a monomer which is equivalent to assuming that fragmentation leads to redistribution of solids in the form of a power law down to this minimum size (see Paper I). However, the moments method (Estrada & Cuzzi 2008) does allow for the lower bound to be evolved as well using an additional moment. The expectation would be that this would even further restrict growth because much of the feedstock for continued particle growth comes from the smaller sizes, especially if the bouncing barrier mass is exceeded.
aggregates (D < 3) will have smaller St than completely compact particles of the same mass. Thus, they will not drift as quickly, will be diffused more easily, and may have the potential to overcome the radial drift barrier.

Compaction Model
Fractal growth of particles by low-velocity sticking of compact monomers of mass m 0 and radius r 0 initially leads to very low density, fluffy aggregates, but as relative velocities increase, these aggregates will be subject to compaction due to collisions with other aggregates. Modeling of this process is generally done with Monte Carlo coagulation codes that allow individual collisions to be followed directly in order to quantitatively assess the amount of compaction that occurs (e.g., Zsom et al. 2010;Krijt et al. 2015;Lorek et al. 2018) in a growing aggregate using a recipe for collisional compaction. However, the computational space is limited in these models, because adding porosity and compaction makes the already cumbersome Smoluchoswki equation multi-dimensional and thus far more numerically intensive to solve. Apart from being impractical for a global model such as ours, it is probably too elaborate given the uncertainties. Okuzumi et al. (2012, see also Homma & Nakamoto 2018) solve the coagulation equation explicitly for the particle mass distribution f (m), and have addressed compaction by deriving a corresponding equation for f (m)V (m) using a volume averaging approximation (Okuzumi et al. 2009), which reduces the computational burden of solving a multidimensional equation to solving two one-dimensional ones. In this approximation, V (m) represents an average or characteristic volume for aggregates of mass m. However, this technique for including compaction does not lend itself readily to our method of moments approach, even if pure fractal growth does (Sec. 2.4.4).
Since we are most concerned with the rate of growth of the largest particle in the distribution, we are only interested in its average volume. Thus we apply a statistical approach where we calculate the change in volume by an estimate of the number of collisions an aggregate experiences due to other aggregates of different sizes in a time step ∆t. The dust population then is characterized by a single D(t) which characterizes the mass dominant particle. On the other hand, because we follow the growth of migrator particles explicitly beyond the fragmentation barrier, migrators can have individually different D. We describe our procedure below. In what follows, unprimed quantities refer to the growing aggregate, and primed to the impactor.
The primary goal of a recipe for collisional compaction is to determine the enlargement factor ψ of a growing aggregate m after a collision. Much like our treatment of bouncing and fragmentation (Eqns. [23] and [26]), at the heart of any growth with compaction scheme is consideration of the collisional impact energy between an aggregate of mass m and other aggregates with m ≤ m: which is the total kinetic energy in the center-of-momentum frame. The amount of compaction an aggregate suffers depends on the ratio between the impact energy and the so-called "rolling energy" E roll , which is defined as the energy required for a single monomer to roll over another by an angle of 90 o (Dominik & Tielens 1997;Blum & Wurm 2000). Physically, low energy collisions (E imp E roll ) lead to perfect sticking with little or no internal restructuring of an aggregate, while for increasing relative velocities (E imp > E roll ), collisions begin to modify the aggregate's internal structure. For equal-sized monomers, Dominik & Tielens (1997) and Blum & Wurm (2000) give an expression for E roll which we write here as where the rolling force F roll = (8.5 ± 1.6) × 10 −5 g cm s −2 and the surface energy density γ 0 = 14 ± 2 g s −2 are measured for uncoated SiO 2 spheres (Heim et al. 1999). The factor γ/γ 0 (e.g., ) allows us a simple way to scale the rolling energy for different materials. For this work, we choose the most recent values of γ 25 and 200 g s −2 (Gundlach et al. 2011;Krijt et al. 2014) for silicates and water ice, respectively. As in the case for fragmentation and bouncing, because our aggregates are composed of both ice and non-ice, we do a weighted average to determine the rolling energy at any disk location (see also Krijt et al. 2018). In the limit E imp E roll , Okuzumi et al. (2009) derived an empirical formula for the merging of two aggregates with volumes V p and V p via a hit-and-stick collision For V p ∼ V p , V HS 3V p which is basically equivalent to pure fractal growth where D 2, and aggregates grow without any net compaction, with a nonuniform internal structure. In the opposite limit (E imp E roll ) for a head-on collision of equal-sized aggregates, the "merged" volume is now explicitly a function of the rolling energy Wada et al. 2008;Okuzumi et al. 2012): Here, V 0 = m 0 /ρ 0 is the volume of a monomer, b = 0.15 is a non-dimensional fitting parameter and N is the total number of monomers contained in the merged volume. For the intermediate case when Okuzumi et al. (2012) adopt and update the analytical fit of Suyama et al. (2012) to complete the porous compaction recipe for the new value of the enlargement factor: HS > V 5/6 p + V 5/6 p and ξ > 1, In order to implement this scheme in our code, we approximate the compaction of a particle statistically, due to the number of collisions it suffers with other aggregates in a given ∆t. The number of collisions that a particle of mass m ≥ m experiences due to collisions with particles of mass m in a time ∆t is estimated using where n(m ) is the number density of particles of m (e.g., see Eq. 75 of Paper I). We then sum over all collisions with aggregates m using Eq. (37) modifying V p after each collision. We do not include all particle masses as this would quickly become cumbersome for a very large difference in mass between the aggregates (see below). Okuzumi et al.
(2012, also see Krijt et al. 2015) note several caveats in using their prescription, in particular, that their formulation is only tested for similar-sized aggregates with mass ratio m /m 1/16. For these authors, as in the case here, this is not a major concern since the majority contribution to growth in our models occurs between similar-sized particles.
For m < m f , the end result of tallying all collisions using Eq. (37) is a new D derived for the largest particle 7 m = m L ≤ m f after each growth step, which is then used in the next growth step for the entire population with m < m f . In this way the particle distribution (m < m f ) has a single power-law dependence for the enlargement factor ψ(m) = ψ d = (m/m 0 ) 3/D−1 consistent with our assumption of a power-law size distribution for the particle f (m) for masses less than the fragmentation mass. On the other hand, the f (m) for migrators (m > m f ) does not generally follow a power law, so that migrators can have both different D and ψ since the growth of the migrator branch of the particle distribution is done explicitly (Paper I).
As aggregates grow to larger sizes, collisions between these and similar-sized aggregates can occur on a timescale longer than a simulation timestep. This is taken into account naturally when solving the full coagulation equation(s), but when the timestep ∆t < t coll , important collisions would be omitted that otherwise might occur at later times. Likewise, even in the case where there are many collisions between different-sized aggregates, N coll most likely will not be an integer. We account for this by rounding N coll down to the next integer value n coll , and we count the final collision (or in the case when N coll < 1, n coll = 0, the only collision), as a single collision, calculate what the contribution ∆ψ would be to ψ, but then do a simple weighting in which ψ(n coll + 1) = ψ(n coll ) + ∆ψ(N coll − n coll ). 7 In the moments method, the largest particle in the particle mass distribution is, by definition, and thus is equivalent to the mass-weighted mean of the distribution (Paper I).  37), and with both collisional and non-collisional compaction effects (solid curve) using Eqns. (39-40), and asssuming αt = 10 −3 . In both cases, aggregates begin as monomers which grow fractally (D ∼ 2) with decreasing internal density until either Eimp E roll and/or gas ram pressure begins to compact the aggregate. When noncollisional compaction effects are included, aggregates will generally continue to compact as they grow, whereas if only collisional compaction is considered, aggregates can remain very underdense even at planetesimal masses.
This assures that there are contributions to the volume change of an aggregate of mass m from the full range of particles with mass m ≤ m we consider during each time step.
In this work, we also consider the compaction for porous aggregates due to gas ram pressure, or via their own self-gravity. In general, our particles do not get large enough for self-gravity to be a factor, because they don't grow sufficiently massive or porous, but we include it for completeness. The external pressure that an aggregate of low internal density can withstand is (Kataoka et al. 2013a;Krijt et al. 2015): which is compared to the pressures associated with the gas and self-gravity (Kataoka et al. 2013b): Implementation of these conditions is quite straightforward in our code. As in the case with collisional compaction, we conduct a check on the largest particle in the particle distribution which may be as large as the fragmentation mass, and for all migrator particles in the distribution beyond the fragmentation barrier. If P gas > P c , we adjust ψ until these quantities are equal, and likewise in the case that P grav is larger. As a demonstration that the above approach works sufficiently well, in Figure 1 we plot the evolution of the internal density ρ int of the largest aggregate under collisional compaction (dashed curve), and compaction due to collisions, ram pressure and gravity (solid curve). For simplicity and comparison, we have chosen a static, Hayashi (1981) MMSN disk model at 5 AU withZ = 0.01, a monomer particle density of ρ p = 1.4 g cm −3 , and for the purpose of collisional velocities we adopt α t = 10 −3 (cf., Okuzumi et al. 2012;Krijt et al. 2015;Homma & Nakamoto 2018). In those models, bouncing, fragmentation and erosion are not included, so that growth to larger sizes for this test is done only using our moments method and thus strictly remains a power law where the largest particle size in the dust distribution is equivalent to the mass-weighted mean (Paper I). We find that a mass ratio of m /m 6.25 × 10 −6 works fairly well to reproduce the qualitative behavior in the evolution of the internal density ρ int = ρ p /ψ as seen by other workers.
Growth proceeds fractally (D ∼ 2) up until the point where E imp ≈ E roll which occurs roughly at m 10 −5 g, whence compaction begins. However, we find that there is a slight delay as to when the "turnover" occurs where the aggregate's internal density remains roughly constant in the collisional-compaction-only case. The slight increase in internal density occurring near m ∼ 1 g is also where r λ mfp , and is likely due to the bridging function we use between

Model
Type Epstein and Stokes regimes (Podolak et al. 1988, see Paper I), while a decrease in ρ int as r (9/4)λ mfp (m ∼ 10 3 g) is because the aggregate has transitioned to the Stokes regime. The sharp increase in ρ int seen at m ∼ 10 4 − 10 5 g is the point where St St η . Here the aggregate first begins to experience type II eddies (section 2.4.1), which cause a sharp increase in the turbulent relative velocity to occur, leading to further compaction . At around m 10 10 g, St 1 and the radial drift velocity begins to decrease, causing the internal density to decrease again.
On the other hand, for the case that also includes non-collisional compaction effects (Fig. 1, solid curve), we note that the point at which E imp = E roll roughly coincides with the point where gas ram-pressure compaction occurs. This is likely because ∆v pg is not too different from ∆v pp in the intermediate regime. From then, the evolution is dictated by the flow regime the aggregate is in. Initially in the Epstein regime, the internal density increases systematically as it becomes more and more decoupled from the gas with ρ int ∝ (∆v pg ) 1/3 . The point at which the particle enters the Stokes regime (r ∼ 9/4λ mfp ) occurs at a higher mass than the collisional-compaction-only case. Once in the Stokes regime, the three different sub regimes are apparent -the internal density increases more slowly when Re p < 1, ρ int ∝ (∆v pg ) 1/3 r −1/3 , but begins to increase more sharply in the intermediate (Re p 1) sub-regime where ρ int ∝ (∆v pg ) 7/15 r −1/5 , and finally ρ int ∝ (∆v pg ) 2/3 when Re p 1. The internal density decreases once St 1, as the particle-to-gas relative velocity decreases (m ∼ 10 10 g) until self-gravity compaction kicks in around m ∼ 10 12 g.

RESULTS
As we did in Paper I, we choose as our baseline model a nebula with a stellar mass of M = 1 M , and an initial disk mass of M disk = 0.2 M . We adopt a fiducial turbulent parameter of α t = 10 −3 , but we explore a range of values from 10 −5 − 10 −2 . The value of the scaling parameter R 0 (Sec. 2.1) is a free parameter which sets the initial surface density profile in terms of the initial disk mass. The rationale for the chosen value of R 0 in the literature has not always been clear, however. A value of R 0 = 10 AU was preferred by Hartmann et al. (1998) based on young star statistics (a value used by Paper I), but it has been set to even larger values of 20 − 60 AU (Ciesla & Cuzzi 2006;Garaud 2007;Brauer et al. 2008;Hughes & Armitage 2012;Yang & Ciesla 2012). Using a different functional form, Cuzzi et al. (2003) chose the equivalent parameter to match the specific angular momentum of the solar system. In this paper we nominally choose R 0 = 20 AU, giving general agreement with the radially compact solar nebula of Desch (2007), but explore a model with R 0 = 60 AU for comparison, and we also consider models with different initial disk masses to show how these choices influence our model simulations. The latter models are compared with our fiducial model in Appendix B. All of our models include the growth barriers, mass transfer and erosion discussed in Sec. 2.4.2. In total, we have conducted 15 simulations for both fractal and compact particle growth which are summarized in Table 2. We focus in this paper on the evolution of the particle mass and porosity distribution, and the ambient properties of the disk. We discuss the evolution of disk bulk composition through the redistribution of refractory and volatile species for these same models in our companion Paper III. When discussing the evolving particle mass distribution over time, we will in some plots overlay approximations of the different barrier masses which we derive from their Stokes number expressions (Paper I, and also Sec. 2.4.2). The first is the bouncing barrier, which we estimate by equating the threshold bouncing velocity to the turbulent velocity: These particles are small and always in the Epstein regime, so a general expression is where for compact particles (D = 3) C b = C 0 and n b = 3/4. For fractal aggregates, for the limits on m defined in Eq. (25). Similarly for the fragmentation barrier, we equate Q f ∼ c 2 α t St and find an expression for m f in the Epstein regime In some instances, the fragmentation may occur in the Stokes regime. In such cases, using the appropriate expression for St (for instance, see Eq. 58 of Paper I) will lead to a better approximation. In terms of mass, the particle where the transition to the Stokes regime occurs is m λ = m 0 (9λ mfp /4r 0 ) D . We note that these barriers correspond to collisions by same-sized aggregates. Finally, the radial drift barrier can crudely be approximated by equating the radial drift time t dr = R/u ∼ (1 + St 2 )/2βΩSt of an aggregate to its growth time t gr = m/ṁ, whereṁ ∼ πr 2 ρ solids ∆v ppS with S < 1 some average sticking coefficient. We nominally chooseS = 0.5, though with bouncing included, a smaller value may be more appropriate which would increase growth times. Determining the appropriate drift mass m d over the various drag regimes can be complicated because t dr = t gr leads to a transcendental equation in St which we solve iteratively up to a maximum Stokes number of unity where t dr has a minimum value. When the growth and drift times intersect, they do so first at a value of St < 1, beyond which the growth time becomes longer. In order to reintersect the drift curve and once again have shorter growth times than drift, a particle would then have to grow sufficiently to reach St > 1 which is generally not possible because drift limits their growth. Hence the radial drift barrier is not a single, well-defined mass, but an entire range of particle growth space in which the particle masses exceed m d , but remain smaller than m f (which is a well-defined, local barrier). More discussion of m d follows in the next section.

Fiducial Model
In this section we compare our fiducial model with α t = 10 −3 for both compact (sa3g) and fractal (fa3g) growth assuming the Gaussian velocity PDF (in Appendix C, we compare with a Maxwellian PDF). In Figure 2, we show the evolution of the particle mass distribution as a function of semi-major axis at 10 4 , 10 5 and 5 × 10 5 years. We chose to cut off the evolutions at 5 × 10 5 years because (a) the particle properties were not evolving quickly (in great part due to the bouncing barrier, e.g., see Zsom et al. 2010), and (b) it appears that planetesimal formation (which is not taken into account in our current model) had proceeded to a significant extent in the inner nebula, and probably at the snowline, by that time (Kruijer et al. 2017). The corresponding color scale shows the solids mass volume density m 2 f (m) per bin of width ∆m as a function of particle mass and disk radius. Plotted in each panel are the estimates of the bouncing barrier mass (brown dot-dashed), fragmentation mass (brown dashed), mass of a particle with the same size as the gas mean free path (brown dotted), and the radial drift mass (brown solid curve) above which particles drift in faster than they can grow (i.e., t gr > t dr ). When present, we also plot the actual fragmentation mass (black dashed), and the mass dominant particle, if it is beyond the fragmentation barrier (black solid curve). When the fragmentation barrier has not been reached, the mass dominant particle is simply the most massive particle in the distribution. The minimum mass in the distribution corresponds to the monomer mass. These curves are most easily studied on-screen with some magnification.  (sa3g) and fractal (fa3g) growth at, from left-to-right, 10 4 , 10 5 and 5 × 10 5 years. Plotted in each panel is a snapshot of the particle mass distribution as a function of semi-major axis. The color scale on the right indicates the particle solids mass volume density per bin of width ∆m, m 2 f (m). The brown curves overlaying the distributions are the mass that corresponds to the mean free path m λ (dotted), and estimates for the bouncing m b , fragmentation m f , and radial drift m d curve masses (above which growth times are longer than drift times, see text and Sec. 2.4.2) as indicated in the leftmost panels. When and where present, the dashed black curve corresponds to locations where the actual fragmentation barrier has been reached. Growth beyond the fragmentation barrier occurs to varying degrees, up to "lucky particles" which are few in number and contain little mass. When the mass distribution extends beyond the fragmentation barrier, the mass dominant particle (particle size containing most of the mass) is shown by the black solid curve. We note that in the fractal models, the m d increases very steeply to masses well beyond the plotted range. We have chosen to have a break in the curves (m f as well) so that they do not cut through the upper right label and legend. The particle size distributions for models sa3g and fa3g after 0.5 Myr can be seen in Appendix A.
Following the particle mass distributions in time from left to right, it can be seen that the fragmentation mass (black dashed curve) in the inner disk 8 regions is reached relatively quickly for both the compact and fractal cases extending out to ∼ 4 − 5 AU after only 10 4 years. After the fragmentation barrier is reached, m f evolves slowly due both to evolving nebula conditions and the inward loss of solids via slow radial drift. The fragmentation barrier is much more difficult to achieve further out in the disk, even after 0.5 Myr, due to a much higher Q f outside the water snowline (which migrates inwards from outside 10 AU to ∼ 4 AU after 0.5 Myr as the disk cools; see Fig. 3, left panel). With a higher Q f , particles can grow to much larger sizes and Stokes numbers (see Fig. 6). This increases their radial drift velocities, allowing even fractal aggregates to drift inwards more rapidly from the snowline outwards to ∼ 20 AU. This leads to a significant drop in the solids surface density in this region (see Fig. 3, right panel). The growth rate in this region begins to significantly slow as the fragmentation barrier (dashed brown curve) is approached, not only due to the decrease in the solids volume density at all sizes (seen as a discernible vertical inflection from yellow to orange), but also because of the bouncing barrier which largely limits particle or aggregate growth from masses m < m b (though in the fractal case, bouncing can still lead to compaction). As a consequence, once the fragmentation barrier is reached, no further growth occurs.
On the other hand, in the inner disk, growth does proceed beyond the fragmentation barrier leading to a population of "lucky particles" or migrators the largest of which generally contain negligible mass (Paper I). At 10 4 years, both compact and fractal growth cases are characterized by a steep dropoff in solids density with masses larger than the fragmentation mass, which are gradually destroyed due to our use of a velocity PDF. As time goes on, however, continued incremental growth by compact lucky particles eventually leads to a secondary peak in the particle mass density at around 1 g (∼ 0.3 − 0.4 cm in radius), as seen at 0.5 Myr (more detail revealing additional growth peaks can be seen in the size distribution, Appendix A). This is due to mass transfer combined with the velocity PDF (see, e.g., Paper I). No such secondary peak is seen in the fractal growth case as time goes on. At first, it may seem as if this behavior is related to whether aggregates are in the Epstein or Stokes drag regime. From comparison with the m λ curves, most of the distribution in the compact growth case is in the Epstein regime, while the fractal aggregates are in the Stokes regime. However, we believe it is simply due to the fractal particles having a much larger m f than the corresponding compact growth case. By entrapping more mass, fractal aggregates limit the abundance of small particles available for "lucky" particles to grow beyond m f . As was seen outside the snowline, growth is mostly restricted by the bouncing barrier, which slows growth enough that a secondary peak cannot be achieved (cf. Fig.  7). It is interesting to note that the largest mass achieved in the inner disk is roughly the same in both compact and fractal cases, suggesting that the growth limit imposed by the combination of bouncing and radial drift is independent of particle porosity.
What we believe the drag regime does determine is the slope of the fragmentation limit curve as seen in these simulations. In the compact growth case (Epstein regime), the fragmentation curve is relatively flat in the inner disk. In the Epstein regime (see Eq. [42]), and for constant Q f , the compact growth m f ∝ (Σ/T ) 3 suggesting the outward decrease in Σ and T (Fig. 3) are similar (though EFs and changes in particle density can affect it to a lesser extent). In the Stokes regime for Re p < 1, m f ∝ R 3 /T 3/2 which should always increase outwards. The same similarly holds true for fractal growth.
The curves for the drift mass m d require some explanation. Quite generally, the "drift barrier" has been associated with St = 1 particles because it is at Stokes unity where a particle's inward drift velocity is maximum (though see, e.g., Birnstiel et al. 2012). However, the drift barrier is more well defined as t gr = t dr , a function of ambient conditions in the nebula (e.g., the local metallicity), showing that the drift barrier is not a sharp transition. For instance, particles can still drift inwards even if their growth times are shorter (t gr < t dr ); conversely, particles can still continue to grow even if their drift times are shorter (t gr > t dr ). Moreover, there are ambient conditions under which the "barrier" can be undefined (t gr and t dr never intersect): the growth time can always be shorter than the drift time for all particle or aggregate masses in the distribution. One can also find a situation in which the growth times are always longer (e.g., for very low ρ solids ). These conditions are reflected in regions where the m d curve is not plotted. Naturally, because the drift barrier is a function of ambient conditions (in particular, changes in the local Z), the m d curve evolves with time. Generally, the drift barrier is only a factor in the icy outer nebula and there, mostly for compact particles. The fractal model cases are a bit more complicated because the local maximum mass fractal dimension D is used in calculating t gr , and thus t gr depends on how compacted the aggregates are. As a result m d can increase very steeply to extremely large masses over a short radial range (this also applies to m f ).
The above conditions are demonstrated by the brown solid curves of Fig. 2. In the compact particle growth model, the drift curve m d is present outside the snowline over the course of the simulation. Growth times are initially shorter than drift times out to ∼ 150 AU and particles grow, but as the simulation progresses, the drift mass curve begins to flatten somewhat and by 0.5 Myr, has evolved below the pre-existing largest particles from ∼ 100 − 150 AU. Much of the material outside the water ice EF, despite being mostly below the radial drift curve, has drifted inwards (see Fig. 5) inside the snowline by this time. For porous aggregates the m d curve seen near ∼ 200 AU increases sharply to masses well beyond the plotted range, so growth times for all aggregates out to this distance are much shorter than drift times. Initially, aggregates are not very compacted, but as they are compacting the m d curve should also decrease. In this α t = 10 −3 model, that decrease is very modest after 0.5 Myr (but will be much more evident in our other α t models, see Sec. 3.2.1 and 3.2.2). We see in the lower panels of Fig. 2 that all aggregates from the water snowline (roughly 4 AU by 0.5 Myr) outwards to ∼ 15 − 20 AU have masses that are well below m d (which is off the top of the plot), so they are still growing rapidly, but through most of this region their growth is limited by fragmentation. Because of the large masses and St of these aggregates, much of the mass in solids outwards of the  Fig. 2. "Flat" regions seen in some parts of the profile are EFs, which migrate inwards considerably over the course of the simulation. The blue curve is a standard MMSN temperature profile (Hayashi 1981) for comparison, which beyond 10 AU is steeper than our profile that has a ∼ R −3/7 dependence due to our flared disk geometry. Inside the water snowline, the mean slope of T has roughly a ∼ R −1 dependence. Right panel: Total solids surface density at the same times. Both compact and fractal growth models are characterized by a sharp decrease in surface density outside the snowline and beyond, as growth proceeds quickly to large particles there which drift in. This is visible by the shift from yellow to orange color in Fig. 2, for fractals at 10 5 yr, and for both by 5 × 10 5 yr. However, the fractal case retains much more material outside starting outside of ∼ 20 AU because aggregates are increasingly more underdense. Note that in the fiducial model, the water snowline migrates inwards from outside 10 AU to ∼ 4 AU after 0.5 Myr. Kinks in Σ solids seen at ∼ 50 and ∼ 125 AU at 0.5 Myr most notable in the compact growth case correspond to enhancements at the CO2 and CH4 EFs. snowline has drifted inwards (indicated by the density step decrease to darker orange colors as one moves outwards from the snowline). In fact, Σ solids for the fractal model is even lower than the compact case between ∼ 10 − 15 AU; however, the reverse is true outside ∼ 20 AU (see below). Conversely, outside of ∼ 150 − 300 AU where ρ solids is very low, the growth times are much longer than the drift times even for monomers, regardless of whether growth is fractal or compact. Inside the water snowline the Σ solids and Z after 0.5 Myr remain similar. In this region no curve for m d is shown for either model because the growth times are always shorter than the drift times for these nebula conditions. Figure 3 (left panel) shows the midplane disk temperature for our fiducial models at the same evolution times as in Fig. 2. The blue curve gives the temperature profile for a standard MMSN (Hayashi 1981) for comparison. Locations in the disk where the radial temperature profile flattens correspond to EFs. Early on at 10 4 years, EFs for water (160 K) to Fe (1810 K) are present. The EFs for the supervolatiles initially lie much further out in the disk, e.g., CO 2 is located at ∼ 105 AU. As the disk cools, these EFs are seen to evolve inward significantly over time. For example, the water snowline moves inwards from ∼ 9 − 10 AU to inside 4 AU by 0.5 Myr, in which time the silicate EF (1450 K) has evolved to within the inner boundary of our computational grid. Though not noticeable here, the CO 2 EF has evolved inwards to ∼ 50 AU, but kinks associated with enhancements in solids at this EF can be seen in the Σ solids profile of Fig. 3 (right panel), which includes both dust and migrators. More detailed analysis of the compositional evolution is discussed in Paper III.
The relatively small differences in the temperature profiles between the two models can be intimately tied to the evolution of the particle mass distributions through their Rosseland and Planck mean opacities, which are plotted in Figure 4. Early in the simulation, there is not much difference between the compact and fractal models. This changes though, and by 10 5 years the opacity in the fractal growth model outside the water snowline drops drastically (this region can be associated with the humps in the particle mass distributions seen in Fig. 2). This large dip reflects that in the fractal case, growth to large sizes is occurring much more quickly as a result of their enhanced cross sections. The larger (and higher St) fractal aggregates have a more rapid inward drift than the compact particles just outside the snowline, so the local Z is significantly lower at 10 5 years (solid curves, Fig. 5; and also Fig. 3, right panel). The sharp opacity drop also leads to a slight, but notably lower temperature (Fig. 3, left panel). By 0.5 Myr, the compact  Fig. 2. In both cases the decrease in opacity outside the snowline (4-5AU) is due to larger particles and material drifting inwards and lowering the local effective Z. However, the fractal case retains much more mass in the outer disk from 40 to 100 AU.  Fig. 2 as a function of semi-major axis. The sharp decrease in the total solids surface density outside the snowline (located at ∼ 4 AU after 0.5 Myr) and beyond seen in Fig. 3 (right panel) is reflected more strongly here as an order of magnitude drop in Z compared to the initial disk metallicity ofZ 0.014. and fractal particles outside the snowline and inside 10 AU have similar opacities, but outside ∼ 20 AU the opacity increases in the fractal case as there is more material retained in the outer disk. Just interior to the snow line (between ∼ 2 − 3 AU), despite the surface densities being similar (Fig. 3, right panel), the opacity is lower in the fractal case (Fig. 4). This is because, while thermal opacities are generally higher for porous particles relative to compact grains of the same size (there are more of them; Cuzzi et al. 2014), here the aggregates in the fractal case are much larger than in the compact case, and the local Z is also somewhat lower, so the overall effect is a lower opacity. In Section 4.2, we briefly discuss scattering properties of porous grains as it relates to the observations. In Figure 6 we explore further the distribution of the masses and Stokes numbers in these fiducial models. The fractal aggregates always have larger masses. At 10 5 years, the largest Stokes numbers in the inner disk are comparable between the compact and fractal models, but outside the snowline the most massive particles in the fractal case have achieved higher St than in the compact growth case even as early as 10 5 years. The resulting higher drift velocities Figure 6. Distribution of Stokes numbers associated with the evolution of the particle mass distributions for the compact (sa3g) and fractal (fa3g) fiducial growth models shown in Fig. 2. Stokes numbers of the distribution are indicated by the color map at right, and labeled (grey) curves of constant St are plotted to aid in visualization. As noted in the text, the more rapid growth of fractal particles beyond the snowline leads to more massive, and higher St particles by 10 5 years compared to compact growth (middle panels). Note that outside of 200 − 400 AU, even small grains have very large St 1 making them essentially unaffected by radial drift. The total mass fraction is very low so growth times are exceedingly long. With careful inspection, one can see the outline of the water ice EF (e.g., ∼ 4 AU at 0.5 Myr), as well as the outline of m λ in the fractal growth case (by reference, the dotted brown curves in Fig. 2; the reader may find it easier to view this onscreen and magnified).
allow the region just outside the snowline to be depleted sooner in the fractal case than for the compact growth case (the dip seen in Fig. 5, in the blue curve at 10 AU). By 0.5 Myr outside the snowline, the compact case has essentially caught up to the fractal case with maximum St 0.03 − 0.04. The compact particle case also shows that the very outer disk (outside of R 20 AU) is becoming more and more depleted in material as mass dominant particles have achieved higher St, and therefore radial drift velocities, relative to the fractal case. Beginning outside of ∼ 20 AU, the model sa3g shows that mass dominant particles have St few × 10 −3 , whereas the Stokes numbers for fa3g drop sharply to values of ∼ 10 −4 and lower further out, explaining why Z and Σ solids there remain high.
Fractal growth in mass is in general faster than compact particle growth over all sizes owing to the larger cross section of aggregates; however, the Stokes numbers of fractal aggregates grow more slowly than compact particles in the earlier growth stages, and in fact remain similar to monomers until compaction begins to set in (see Sec. 2.4.4, and Fig. 1). Once significant compaction occurs, the fractal aggregate St can grow more quickly than compact particles. This is evident in the region outside the snowline to ∼ 20 AU. The St of fractal aggregates and compact particles are similar at 0.5 Myr, and also at 10 4 years, but clearly at 0.1 Myr the fractal aggregates have already achieved much higher St values. This is reflected in the much lower opacity in model fa3g (see cyan curves at 0.1 Myr, Fig. 4). By 0.5 Myr, the opacities (and surface densities, black curves Fig. 3, right panel) are similar for both models. Outside of ∼ 20 AU, the difference in Stokes numbers in Fig. 6 between sa3g and fa3g then is largely due to slower dynamical times and decreasing ρ solids with distance (both which influence growth rate) such that fractal aggregates are not yet compacted enough to be comparable to the compact growth model particles. Thus a characteristic of fractal growth is that aggregates can be retained further out in the nebula for longer periods of time, to provide source material for potential planetesimal formation there (see black solid curve in Fig. 5). Figure 6 also helps explain an effect seen in previous plots, of a transition that occurs around ∼ 200 − 400 AU in these models where a stark depletion (though relatively subtle in total mass) in material is evident, both in the opacity (Fig. 4) and metallicity (Fig. 5) plots. This transition is due to the steep decrease in the gas surface density. On the one hand, the steep "edge" of the gas disk leads to a strong pressure gradient facilitating very rapid inward drift. On the other hand, even small aggregates quickly achieve large Stokes numbers at these low gas densities and can quickly become immobile, remaining so until the gas density increases as the disk spreads outwards at later times. All panels in Fig. 6 demonstrate this evolution. Inside the transition, one can see a spike of increased growth produced by the increased solids fraction having St ∼ 10 −3 − 10 −2 . Outside the transition, despite particles being immobile with St > 0.1, there is very little material, and thus growth times are exceedingly long regardless of whether growth is compact or fractal.
Overall, a general result of these fiducial α = 10 −3 simulations is that fractal, porous growth is characterized by massdominant aggregates that are as much as ∼ 2 − 4 orders of magnitude more massive than their compact counterparts, even as their internal densities are orders of magnitude smaller (Sec. 3.3). The masses of "lucky" particles, however, are not too different between cases. We find that regardless of whether growth is fractal or compact, the mass-dominant aggregates are in the fragmentation regime in the inner disk regions, but for the larger Q f used for icy particles beyond the snowline, aggregate or particle masses and St can be much larger and much more prone to rapid radial drift. Slowed growth as a result of the higher fragmentation mass threshold (which limits the amount of fine dust), and steady loss of material inwards, conspires to keep aggregates in the drift-dominated regime. Both compact and fractal particles quickly evacuate the region between the snowline and ∼ 40 AU, but porous aggregates survive longer, and in greater abundances, outside of ∼ 40 AU.

Strong Turbulence αt = 10 −2
In Figure 7 we present snapshots of simulations comparing compact (sa2g) and fractal (fa2g) growth models using a higher α t = 10 −2 at 10 4 , 10 5 and 5 × 10 5 years. Due to the higher relative velocities and collisional frequencies induced by the stronger turbulence, growth times are more rapid and both the bouncing and fragmentation barriers are reached at much smaller sizes, with the latter generally lying within the Epstein regime in both models. As a result of lower m b and m f , and similar to what was found in Paper I, lucky particles are even more "lucky" under stronger turbulence relative to our fiducial model for three reasons: first, the fragmentation barrier occurs at smaller m f and small particles are constantly replenished from which lucky particles mostly grow; second, many decades of possible growth are left before the barriers imposed by the combination of bouncing and radial drift; and third, higher relative velocities between the largest and smallest particles increases the efficiency of mass transfer (Windmark et al. 2012a). Much like our fiducial case α t = 10 −3 , the masses of the largest lucky particles are not too different between the compact and fractal models.
For this α t , the compact growth case (top row) shows that a strong secondary peak appears inside the water ice EF relatively early on (after 10 4 yrs), which persists throughout the simulation (another tertiary growth peak is more easily seen in the size distribution, see Appendix A). This is due to mass transfer (see Sec. 2.4.2) and our velocity PDF, significantly increasing the mass of the mass-dominant particles beyond the fragmentation mass (black solid curve). This is consistent with the behavior seen in the compact growth simulations with α t = 10 −3 , where significant growth beyond the fragmentation barrier occurs with a well-defined secondary peak appearing by 10 5 years (Fig. 2). We see a slightly smaller secondary peak in the fractal case as well (third row, Fig. 7; also see Appendix A), whereas no distinct secondary peak was seen in the fiducial model (the largest achieved particle masses were not far beyond m f ). We note that the fragmentation curves are flatter at this α t , (also see below) consistent with being for the most part in the Epstein regime where m f depends on the behavior of both T and Σ (cf. Fig. 2, and Sec. 3.1) 9 .
Beyond the water snowline, the compact particle growth case is fragmentation limited out to 50AU by 0.1 Myr, and over almost the entire extent of the disk by 0.5 Myr (though lucky particles can be found as far out as ∼ 15 AU), as the disk gas rapidly evolves outwards. In contrast, the fractal case is fragmentation limited to about ∼ 100 AU or so; outside this location, growth to larger sizes is still proceeding slowly, and may still reach the fragmentation limit at later times. The strong diffusion and advection associated with α t = 10 −2 coupled with lower mass particles prevents Figure 7. Comparison of models with αt = 10 −2 for both compact (top panels) and fractal (bottom panels) growth at, from left-to-right, 10 4 , 10 5 and 5 × 10 5 years. Plotted in each panel is a snapshot of the particle mass distribution as a function of semi-major axis. 1st and 3rd row: Particle solids mass volume density per bin of width ∆m. 2nd and 4th row: Stokes numbers. The particle size distributions for models sa2g and fa2g after 0.5 Myr can be seen in Appendix A. the sharp drop in solids surface density seen outside the snowline in the fiducial case, and indeed Σ solids (see cyan curves, Fig, 9) decreases in a similar manner to the gas over much of the disk.
The mass-dominant particles for both compact and fractal particles mostly lie below the radial drift mass curve where it is defined (as in the fiducial model). By 10 4 years the m d curve drops especially sharply near ∼ 200 AU, crossing to below the mass-dominant particle mass. Unlike the fiducial cases, however, this location migrates to much larger R with time as the nebula gas advects outwards. Thus particles outside 200 AU, which are initially growing slowly or not at all, eventually find themselves below m d and the increasing gas density allows for both significant growth and inward drift, even very far out in the disk. Also by 10 4 years in the compact growth case, m d is defined everywhere inside the snowline, but as time progresses, the growth times become faster than the radial drift times due primarily to the increase in Z in the inner disk, causing the m d curve to vanish. In the fractal case, the m d curve is only mildly less restrictive than the fiducial case, at least initially, with growth times always faster than drift inside ∼ 200 AU (after 10 4 yrs). By 0.5 Myr, the nebula gas density has increased so much even at 1000 AU that St < 1 everywhere. The radial drift mass curve has lowered and flattened significantly as well in both models, but m d remains still far larger than the largest particles or aggregates, and mostly well above the fragmentation curve.
With regards to growth, the α t = 10 −2 cases indicate (more clearly in the fractal growth model) that the peak in the particle mass distribution (which occurs outside the water ice EF) has been reached sometime earlier in the simulation than for α t = 10 −3 . The fiducial models have the largest peak masses after 0.5 Myr, while the peak masses in fa2g (and sa2g) occur between 0.1 − 0.2 Myr. The peak in the mass distribution for the α t = 10 −2 models (which outside the snowline for the mass dominant particle is m f ) notably begins to decrease (in fact, by over an order of magnitude) due to the rapid outward expansion of the disk gas. At the same time (though not discernible from Fig. 8), the fragmentation St f has grown larger from 0.1 to 0.5 Myr with the increase in St f ∝ T −1 (Eq. 27) almost all due to the disk's lower temperature later in the simulation. The decrease in the fragmentation mass m f ∝ (Σ/T ) D/(D−2) (Eq. 42, see also Birnstiel et al. 2012) is because Σ has decreased much faster relative to T over the same timescale. This behavior is different from the fiducial model where the peak in the particle mass distribution has steadily increased, because the gas disk has not spread out nearly as much and Σ remains relatively large. However, it suggests that at later times the peak masses for the α t = 10 −3 models will eventually become smaller as the disk spreads and the gas density locally decreases. The α t = 10 −2 case is also interesting in that the stopping times for both the compact and fractal particles are closer to the Kolmogorov eddy time -that is, St ∼ St η = Re −1/2 (Sec. 2.4.1). As discussed in Sec. 2.4.2, this means that the approximation for the fragmentation Stokes number given by Eq. (27) is inaccurate. Although y * = 1.6 (Ormel & Cuzzi 2007), one can no longer ignore terms in Re −1/2 in Eq. (15), which increases the estimate of the fragmentation mass. This has been taken into account in the curves shown in Fig. 7. Finally, the difference in peak masses of the particle mass distributions between the compact and fractal growth cases (here we mean m f in the inner disk and the largest particles in the outer disk) is more pronounced at α t = 10 −2 than for the fiducial α t = 10 −3 case -now ∼ 4 − 8 orders of magnitude by 0.5 Myr. This is due to a combination of more underdense aggregates (Sec. 3.3) and overall smaller m f , which leads to smaller radial drifts. This is also evident in the Stokes numbers (2nd and 4th rows of Fig. 7). The Stokes numbers achieved in the compact and fractal growth models are quite similar to each other, but about an order of magnitude smaller than those for the fiducial case in Fig.  6. 3.2.2. Weak Turbulence αt = 10 −4 and αt = 10 −5 The situation is quite different when choosing a smaller α t = 10 −4 , the compact (sa4g) and fractal (fa4g) growth models of which are shown in Figure 8 (the corresponding size distributions at 0.5 Myr can be found in Appendix A). In this case, the bouncing and fragmentation barriers occur at much larger masses (sizes) which are mostly in the Stokes drag regime out to 10 − 20 AU. Inside the snowline, and as a result of the larger fragmentation mass (which consolidates more mass in fewer particles), we see considerably less growth beyond the fragmentation barrier in both models (rows 1 and 3) than for higher values of α t . This is in spite of the larger bouncing barrier threshold which by itself would allow, at least early on, a wider range of feedstock to grow from. We believe the reason for this is that the rate of growth has slowed sufficiently (as feedstock at smaller sizes becomes more depleted) that a rough growth "ceiling" has been reached where particles see little or no benefit from the velocity PDF or mass transfer, both of which are key for lucky particle growth.
Unlike the model with α t = 10 −2 , the outer disk beyond the water snowline in the compact growth case for α t = 10 −4 effectively becomes drift limited after 0.5 Myr. The radial drift curve for m d has gradually evolved (over much of the Figure 8. Comparison of models with αt = 10 −4 for both compact (top panels) and fractal (bottom panels) growth at, from left-to-right, 10 4 , 10 5 and 5 × 10 5 years. Plotted in each panel is a snapshot of the particle mass distribution as a function of semi-major axis. 1st and 3rd row: Particle solids mass volume density per bin of width ∆m. 2nd and 4th row: Stokes numbers. Though it is more apparent in the compact particle models, the St of particles outside the snowline tend to be relatively constant over a large radial extent (Paper I). For these models, the water snowline is located near ∼ 3 AU after 0.5 Myr. The particle size distributions for models sa4g and fa4g after 0.5 Myr can be seen in Appendix A. Figure 9. Total solids surface density for the sa2g, sa4g and sa5g (dashed curves) and fa2g, fa4g and fa5g (solid curves) models after 0.5 Myr. Peaks in Σ solids can be associated with solids enhancement outside EFs. Also plotted in grey are the fiducial models for direct comparison. region) well below the most massive particles in the distribution. For instance, at 0.1 Myr, the mass dominant compact particles are well below m d , but by 0.3 Myr (not shown) the m d curve has dipped below the largest particles outside ∼ 40 AU, and under these conditions we can say growth is "drift limited". Note that the St associated with m d is now 1. The reason of course is because the growth times have increased greatly due to the much lower Z in the region. As mentioned in Section 3.1, particles can still continue to grow beyond m d , but material begins to drain (drift inwards) from the region more quickly than it grows. Particles Stokes numbers would need to exceed unity to reverse this trend, but they can never do so (Sec. 2.4.2). Indeed, after 0.5 Myr much of the disk material has migrated into the inner disk (see Fig. 9) and further growth is quashed well before reaching the fragmentation barrier. Some enhancements in Σ solids are seen outside the EFs for CO 2 ( 50 AU) and CH 4 ( 125 AU), which observationally may show up as supervolatile bands, particularly in young Class 0/I disks (e.g., Segura-Cox et al. 2020; see Paper III, for more discussion). The systematic draining of the outer disk material into the inner disk is an expected result for low-α t , compact particle growth models (see Paper I, and references therein).
In the α t = 10 −4 fractal case, much like we saw for the fiducial model (fa3g), rapid growth just outside the H 2 O EF leads to a sharp decrease in surface density in a band out to ∼ 10 − 20 AU, much earlier on (third row, 10 5 years), but much more material remains outside 20AU compared to the compact growth model (see Fig. 9, solid and dashed black curves). Initially, in the outer disk most of the mass is not fragmentation limited, though by 0.5 Myr the most massive aggregates have approached m f out to ∼ 15 AU. At that same time a significant fraction of the particle mass distribution lies above the m d curve (between ∼ 8 − 20 AU), due to a dip in the solids surface density which significantly lengthens the growth times. Interestingly, the maximum particle masses achieved in the compact and fractal cases in the inner disk are roughly similar to each other (note that the fragmentation masses themselves only differ by ∼ 10), but the maximum masses achieved outside the snowline differ by 4 orders of magnitude with the fractal case reaching particle masses as high as ∼ 10 9 g. Although large, this is still not massive enough for compaction by self-gravity. Using our nominal particle density (1.52 g cm −3 , Sec. 2.4.5), this is a ∼ 75 m size object with mass equivalent to roughly a 5 m radius compact particle, whereas in the compact growth case, the largest radius particle is ∼ 25 cm (see Appendix A, Fig. A1).
The similarity in the inner disk maximum particle masses can also be explained with the help of the distribution of Stokes numbers shown in Fig. 8 (2nd and 4th rows). Both the compact and fractal case are in the Stokes regime, with m f > m λ and Re p < 1, so neither the largest compact nor fractal aggregate particles have any dependence on the local T and Σ, and in fact the Stokes numbers are roughly the same for the two cases. Outside the snowline after 0.5 Myr, a typical mass-dominant particle Stokes number for the compact case is St ∼ 0.02 − 0.04, but for the fractal case St ∼ 0.1 − 0.15. The disparity is because the maximum achieved mass in the compact case is severely drift limited, whereas the fractal case is mostly fragmentation limited, and drift only starts to become a limiting factor much later in the simulation. Figure 10. Comparison of models with αt = 10 −5 for both compact (top panel) and fractal (bottom panes) growth at 5 × 10 5 years. Plotted in each panel is a snapshot of the particle mass distribution as a function of semi-major axis. Descriptions the same as in Fig. 2. The fractal case still maintains a region outside of ∼ 20 AU with a significant surface density of solids (see Fig. 9) that is increasingly below the fragmentation curve, but becomes drift limited again around ∼ 80 AU. Material just outside the (steeply decreasing) radial drift curve in the rightmost panel of Figure 8 (3rd row) has significantly larger Stokes numbers than the material interior to it (see the inflection at ∼ 10 3 g in the 4th row). Thus material is drifting in rapidly there, and corresponds nicely to a steep drop in Σ solids (Fig. 9, black solid curve). This is not seen in the higher α t models where disk spreading (and smaller St) maintains all masses well below the radial drift curve (Fig.  7, 3rd row). Even though aggregates have m < m d in the region between ∼ 20 − 80 AU, we expect that even this region will become depleted at later times. In both models, outside of ∼ 300 AU, physically small particles remain trapped because of their large St (yellow colors), but these particles may eventually grow and drift as the gas disk slowly spreads outwards, decreasing their St as the local gas density increases.
It should be noted that the high Stokes numbers outside the snowline combined with low α t = 10 −4 turbulence means that in this model eddy-crossing effects are just starting to become important (Sec. 2.4.2). Indeed, the useful parameter St/α t is ∼ 300 (see Sec. 2.4.2) for the compact case, and ∼ 1200 in the fractal case, both comparable or larger than the local headwind parameter β −1 , a necessary condition for y * > 1.6 (Fig. 11, black curves; see Sec. 2.4.2 and also Ormel & Cuzzi 2007;Jacquet et al. 2012;Jacquet 2014, for more discussion on the significance of St/α t ). For the fractal case, we note that the fragmentation curve lies slightly above the particle mass distribution. This is because the eddy-crossing correction leads to a significantly different y * and collision speed such that the simple expression in Eq. 27 used in these fragmentation curves is an overestimate.
We now turn to compact and fractal growth models for α t = 10 −5 . The rapid loss of material via radial drift is even stronger as can be seen in Figure 10 where we show the solids mass volume density (left) and St (right) for the compact particle (top) and fractal growth (bottom) simulations after 0.5 Myr. In this case St/α t ∼ 1 − 3 × 10 4 for the compact and fractal models, considerably larger than β −1 (Fig. 11, salmon curves) compared to the α t = 10 −4 case where St/α t is smaller. Here, by 0.5 Myr the compact particle growth simulation has lost almost all of the outer Figure 11. Pressure gradient β as a function of semi-major axis after 0.5 Myr for all models discussed thus far. Solid curves correspond to fractal models, and dashed to compact growth models. The blue curve corresponds to the MMSN nebula. Though the gradient shows considerable variation in the most inner disk regions, at no point is β < 0, so appreciable pressure bumps do not develop over the simulation time of these models. Figure 12. Midplane temperature profiles after 0.5 Myr for all models discussed thus far. Solid curves correspond to fractal models, and solid curves to compact models. The blue curve corresponds to the MMSN nebula. disk and much of the inner disk material, leaving distinct bands of trapped solids at the water and supervolatile EFs. The innermost region ( 0.8 AU) is a rich band of organics. This strong band arises, however, because we have treated organics evaporation and recondensation as a reversible process, which is likely not accurate (see footnote 11 and Paper III, for more discussion of compositional evolution). In the fractal growth case, the solids surface density profile changes considerably with time, looking similar to the α t = 10 −4 case at 0.2 Myr outside the snowline (0.2 Myr snapshots not shown for either α t ), but by 0.5 Myr most of this material has been lost inwards due to even faster drift than for α t = 10 −4 , and only leaving material remaining trapped in the bands outside their corresponding EFs (Fig. 9, salmon solid curve). In the region outside the snowline, fractal aggregates have achieved masses 10 9 g, with St ∼ 0.2 − 0.3 (both about a factor of two larger than the α t = 10 −4 case), but the local Z has decreased strongly to Z ∼ 2 × 10 −4 .
In Figure 12 we plot the midplane temperature profiles after 0.5 Myr for all the models discussed so far. The temperature profiles for the fiducial model (grey curves) are included for comparison. As one might expect with the stronger diffusion and advection associated with the higher α t = 10 −2 case (cyan curves), the temperature profile is less steep overall than in the α t = 10 −4 , 10 −3 models. In the latter cases, a steeper gradient in the temperature builds up inside the snowline (which lies further inward than in the α t = 10 −2 case) due to the inward drift of material there that increases Σ solids (see Fig. 9), and thus the local Z and opacity. The pressure gradient is very much affected by the steep gradients in T , as can be seen by the variations in Fig. 11. However, because aggregates are smaller and more coupled to the gas flow for α t = 10 −2 , there is less variation in both Σ solids and β over the disk (cyan curves in Figs. 9 and 11). It should be noted that the variations in β across the model set are not sufficient to lead to pressure bumps that could potentially trap particles. When comparing these temperatures across the range of α t , the fiducial model (grey curves) is hotter inside the snowline than both the higher and lower turbulence cases after 0.5 Myr, which may seem counterintuitive. As it turns out, an optimal combination of solids enhancement across the snowline (which is highest in these models after 0.5 Myr), the mass/size of the dominant aggregate, and the ambient gas surface density keeps the region inside the snowline hotter longer for the fiducial model. In the α t = 10 −2 case, the dominant aggregate masses are smaller than the fiducial model, increasing the opacity, but Z is not much enhanced by drift, and Σ is much lower because of fast viscous evolution. In the α t = 10 −4 case as in the fiducial model, the enhancement in Z in the inner disk is also high due to rapid radial drift, and Σ remains high due to slower viscous evolution, but the larger dominant aggregate masses allowed by α t = 10 −4 lead to lower opacities. Thus α t = 10 −3 seems to be a "sweet spot" from the standpoint of retaining a hot, massive inner nebula. In the α t = 10 −5 case (salmon curves), Z is very low because of inward drift out of the inner nebula such that the opacity is very small, other than in regions just outside an EF. In Paper III we discuss the EFs in more detail.

Evolution of Internal Aggregate Densities
In Figure 13 we plot the particle internal density ρ int = ρ p /ψ (recall ψ is the enlargement factor, Sec. 2.4.4) for our fiducial model (Sec. 3.1, middle row), model fa2g (top row) and model fa4g (bottom row) to better visualize the underdense nature of the fractal aggregates across these simulations. In all models, aggregates initially grow fractally with D 2 but eventually begin to compact as collisions become sufficiently energetic and E roll is exceeded, and/or more often, gas ram pressure sets in (Sec. 2.4.5). In these plots only the monomer mass (the minimum mass in the distribution) has a compact particle density ρ p , which is a function of semi-major axis. The locations of EFs where changes in aggregate density occur are in most cases discernible upon inspection (e.g, the water snowline is located between ∼ 3 − 4 AU in these models after 0.5 Myr), even for CO (e.g., middle row at ∼ 350 AU). Paper III discusses the radial distribution of aggregate bulk composition in more detail.
In the inner disk, compaction occurs much more quickly due to the faster dynamical times, higher gas surface density and a lower rolling energy threshold. Once the fragmentation barrier is reached in the inner disk, internal densities vary within a range of ρ int ∼ 10 −3 −10 −2 g cm −3 . Outside the snowline where E roll can be about an order of magnitude larger (Sec. 2.4.5, though see Sec. B regarding "cold ice"), the most massive aggregates are more underdense with ρ int ∼ 10 −4 − 10 −3 g cm −3 (or even smaller). The maximum achieved aggregate masses inside the snowline differ by ∼ 1 − 2 orders of magnitude between these fractal models, but the outer disk has maximum aggregate masses ranging from ∼ 10 5 − 10 9 g. For the given densities, this means that the effective radius of these fractal aggregates ranges from ∼ 3 − 75 m, with the upper value corresponding to model fa4g at ∼ 14 AU. Even the high end of this range of aggregate mass is still too small for self-gravity to be important (note that for α t = 10 −5 , the largest aggregate has a mass roughly twice this value). By equating P c = P grav (see Sec. 2.4.5), we can estimate what the mass needs to be for self-gravity compaction if one knows the current fractal dimension D: At 0.5 Myr and at 14 AU in the α t = 10 −4 model, D 2.6 and E roll ∼ 7.2 × 10 −9 g cm 2 s −2 which predicts m grav ∼ 7 × 10 10 g. However, it is unlikely that icy aggregates will ever achieve such a mass, as they are not only fragmentation limited, but radial drift limited as well (Fig. 8, 3rd row, 3rd panel).
While the aggregates outside the snowline appear quite large, and have very large cross sections, their rapid growth to such large sizes has come at the cost of rapidly decreasing Σ solids of their "feedstock", which slows down their growth rate. Despite fractal growth leading quickly to very underdense aggregates, which drift more slowly than a compact particle of equivalent mass, the local nebula conditions don't allow for continued incremental growth to larger sizes, and conditions are such that the local Z is significantly reduced. In Section 4.2, we discuss how these highly porous aggregates relate to the observations. Figure 13. Evolution of the internal density ρint = ρp/ψ of aggregates for models fa2g, fa3g and fa4g shown in Figs. 2, 7 and 8. Despite being relatively massive (∼ 10 5 − 10 9 g), the large aggregates outside the snowline are quite underdense, though they continue to be compacted. In the inner disk, particles were initially underdense, but begin to be compacted earlier on and more vigorously. Transitions in density due to compositional changes at EFs can be seen as vertical bands in the plots. Darker vertical bands (lower density particles) initially seen at ∼ 250 AU are solid methane (see Paper III). This region of the disk also exhibits rapid inward particle drift where St ∼ 1 (a location that moves outwards with time due to viscous spreading of the disk), which manifests as the peak in particle mass seen at 150-200 AU until 10 5 years. The particle size distributions for these models after 0.5 Myr can be seen in Appendix A.

Implications for Planetesimal Formation
Though we do not currently include an algorithm for planetesimal formation in our code, we can say with confidence that planetesimals do not likely form via incremental growth over the phase space we have explored. We can also explore whether conditions for leapfrog planetesimal formation by SI can be satisfied. These conditions are determined almost entirely by a combination of St and α t through the midplane solids-to-gas volume density ratio = ρ solids /ρ ≈ Z (α t + St)/α t (Umurhan et al. 2020, and Paper I). To that end, Figure 14 shows a more detailed temporal evolution of (R) from 0.001 to 0.5 Myr, for a subset of compact particle (left panels) and fractal aggregate (right panels) growth models for turbulent intensities of α t = 10 −3 , 10 −4 and 10 −5 . Models for α t = 10 −2 are not shown because for this large diffusivity, the well-coupled nature of the particles and aggregates to the nebula gas prevents any significant enhancements anywhere in the disk. Indeed, the distinction in the evolution of nebular properties between fractal aggregate and compact particle growth models (apart from their difference in particle mass) begins to disappear with increasing α t over the epoch simulated here.
The top panels show that, at first blush, the compact and fractal growth simulations for α t = 10 −3 (our fiducial case) initially look very similar. At a closer look, the more rapid growth of fractal aggregates in the fractal growth model fa3g just beyond (and to some degree, inside) the evolving snowline leads to small differences in both the inner and outer disks (Sec. 3.1). However, by 5 × 10 4 years (the H 2 O EF is located at ∼ 7 AU) a deviation is seen where model fa3g has achieved a higher due to faster growth (and thus drift) of material from 20 AU. By 0.1 Myr though, the compact growth model sa3g has largely caught up, and the variations of inside the snowline look relatively similar after that, except that the enhancements at the H 2 O EF for model sa3g are somewhat larger. The reason for this has to do with the much lower St of fractal aggregates beyond 20 AU compared to compact particles, which allows mass to be retained in the outer disk for longer periods. Thus, in model fa3g, less material has been processed across the snowline because material beyond ∼ 20 AU has not had time to drift in yet. On the other hand, model sa3gwhose particle St are about an order of magnitude larger than in model fa3g outside of 20 AU has had a higher, more continuous inward mass flux of material, producing a larger enhancement both outside the snowline, but especially outside the organics EF (the higher of the two notable peaks, e.g., at 1.95 AU after 0.5 Myr compared to 2.9 AU for water, see also Paper III) which is further inward. Nevertheless, it is expected that the material in model fa3g will eventually drift in at later times. Overall the differences between the fractal and compact growth cases are small and = 1 is never approached. The simulations for α t = 10 −4 (middle panels) show a sharper deviation between compact and fractal growth. Lower turbulence values allow growth to larger St, and thus faster radial drift rates with more pronounced systematic loss of material from the outer disk to the inner disk, as is evident in model sa4g. Local peaks in near 100 AU and evolving inward caused by rapid inward particle drift are due to a combination of a strong pressure gradient and increasing St of material further outward. By 0.5 Myr (blue curve), most of the material outside the snowline (located at ∼ 2.7 AU) has been depleted in model sa4g, except for some distinct bands of material trapped at the H 2 O EF and the supervolatile EFs at 50-80AU (see Paper III, and Sec. 3.2.2). The more rapid drift and growth to larger sizes further out in the disk in the compact growth model explains the larger contrast in the enhancement in outside the snowline between these α t = 10 −4 models (and compared to sa3g and fa3g), which occurs simply as a result of more mass being processed across the water snowline in sa4g. In the fractal model, solid material is mostly being tapped from inside ∼ 20 AU, so model fa4g still contains a significant amount of solid mass outside of ∼ 20 AU (see Fig. 9) owing to smaller St numbers and slower drift. As a result, much larger enhancements are seen at the water snowline ( 0.1) in the compact growth case with the largest peak around 0.2 Myr, and even after 0.5 Myr a strong peak in is maintained, but in the fractal model the magnitude of the peaks are notably less in comparison. Moreover, in the compact growth model there is a large reservoir of water vapor inside the snowline (see Fig. 6 of Paper III) that contributes to keeping the peak in solids strong (and becoming overwhelmingly water ice composition) via slow outward diffusion, but also as the disk cools and the H 2 O EF evolves inwards. The reservoir of water vapor in fa4g inside the snowline is also substantial, but lower in magnitude than sa4g (and still of relatively diverse composition) again due to less mass being processed across the EF. Finally, as in the fiducial α t = 10 −3 model, the strongest enhancement occurs outside the organics EF (which is a questionable feature) 10 .
Generally speaking, as α t decreases in our models, the resulting growth to larger aggregate or compact particle masses with larger St, eventually leads to more rapid drift and depletion of much, if not all of the region outside the water snowline. However, because particles (and their St) are significantly smaller inside the snowline, drift is slower there, allowing for a much larger buildup at the (questionable) organics EF in the compact case, or even the silicates EF as can be noted in the fractal model for simulation times 0.3 Myr. Indeed, inside the water snowline, one finds sustained regions with 0.1 over time, but St remain relatively small. Outside the water snowline, the can also increase to values exceeding 0.1 (in sa4g) early on as Stokes number growth outpaces any decrease in the local Z, but as the masses (and St) increase and with them the radial drift speeds, the local Z drops precipitously (and eventually Figure 14. Evolution of the midplane solids-to-gas ratio = ρ solids /ρ over time, for compact particle (left panels) and fractal aggregate (right panels) growth models with αt = 10 −3 , 10 −4 and 10 −5 . The location of the water ice snowline at the corresponding simulation times is indicated by the vertical dotted lines. Only the lowest turbulence model sees a very short period where > 1 at the snowline, with St 0.1. Models with αt = 10 −2 are not shown because there is no significant enhancement in anywhere. Planetesimal formation by Streaming Instability requires midplane > 1 (Paper I; Umurhan et al. 2020). as well) to values considerably below 0.1 at the same time. That is to say, there is never a situation where planetesimal formation by SI can occur at any of the EFs in the inner disk or at the water snowline in either model (see Umurhan et al. 2020, and below for more discussion).
One interesting wrinkle appears when α t = 10 −5 (bottom panels). These extremely low global turbulence models allow growth to much larger masses, which show an even faster radial drift. However, we find for both the compact particle (sa5g) and fractal aggregate (fa5g) growth models that there is a brief period of time near ∼ 0.2 Myr that lasts ∼ 2 × 10 3 − 2 × 10 4 years during which a buildup in the solids surface density at the water snowline leads to > 1. The corresponding Stokes numbers (particle radii) are ∼ 0.05 (38 cm) and ∼ 0.35 (150 m radius). These two conditions together would satisfy the condition for SI. Like the models for α t = 10 −4 though even more so, the spike in is due to the sustained high mass influx of material drifting in from the outer disk, and partially because of the inward migration of the water snowline as the disk cools. In sa5g and fa5g, the relatively faster cooling time allows for the even stronger enhancement of water vapor inside snowline to condense outside the snowline. For instance, in model fa5g, the local solids mass density ratio increases by nearly two orders of magnitude (black dashed curve), to a composition of nearly pure water ice (the compacted densities are near ∼ 0.9 g cm −3 ). This pure ice composition is discussed in Paper III, where we explore the bulk composition of planetesimals if they were to form anywhere in Figure 15. Theoretically predicted SI growth timescales in a turbulent disk (colors) with superposed particle St from our models (symbols). The growth timescales are for the fastest growing SI modes, and are based on the model of Umurhan et al. (2020) and shown in shades of color in units of local orbit times. A sharp diagonal boundary in (αt, St) phase space that cleanly separates the blue and green regions represents the = 1 line, with smaller values above the line. The values shown are specifically valid for a local Z = 0.01 just outside the snowline (i.e., T ≈ 160 K). Various growth barriers and turbulence constraints discussed in section 4.1 are shown as lines with arrows indicating putative forbidden regions in parameter space. The plotted St symbols correspond to particle radii dominating the particle mass, from our models for four different values of αt. Solid filled symbols represent non-porous particles while stippled symbols represent porous fractal aggregates. Two classes of particles are shown at αt = 10 −3 : the yellow symbols are for nominal "sticky" water ice, and the orange colored species denote particles whose fragmentation energies (Q f ) are relatively low, representing cold water ices as discussed in section B (e.g., Musiolik & Wurm 2019). The particle St predicted in our growth results generally fall within the < 1 region of parameter space wherein SI growth times are predicted to be very long, in agreement with all published numerical simulations (Umurhan et al. 2020). Only the fractal growth model (fa5g) for an extremely low αt = 10 −5 falls below the = 1 line. The compact growth model (sa5g) can also satisfy this condition, but only at higher local Z (the = 1 curve shifts upwards for increasing values of the metallicity, see Umurhan et al. 2020).
the disk. After this brief period where spikes, decreases falling below unity. By 0.3 Myr, most of the outer disk in model sa5g has collapsed into the inner region and is being lost to the central star, whereas significant solids mass remains in the disk even after 0.4 Myr in model fa5g. By 0.5 Myr though, most of the disk in the fractal aggregate model is depleted as well.
Some recent models of compact particle growth have indeed assumed a turbulent intensity as low as 10 −5 , but only for the particles, while allowing for far larger turbulence driving disk evolution (Drażkowska & Dullemond 2018). However, these authors find planetesimal growth (based on satisfying the condition ≥ 1 and St ≥ 0.01) even for particle layer α t as high as 10 −3 , so while the behavior we see here could be connected to the fact that SI can form planetesimals at the snowline in those models, it is hard to judge the effects and viability of using different values for the turbulence for the particles and gas, and for such high values of α t where we do not see the conditions for SI satisfied 11 .
Finally, we note that if conditions for SI are satisfied at the snowline, the simulations may not properly reflect what happens afterwards since a period of planetesimal formation may produce bodies which are effectively immobile over long time scales. These planetesimals might sweep up inwardly drifting material, or form a planetary core, either of which could be a barrier to further inward mass flux. Thus the results for the extremely low value of α t = 10 −5 should only be cautiously extended beyond those times (see below). These results for different α t highlight the complexities to be found in these global evolutions, and the critical role played by α t .
What do these predictions imply about SI activity? In Figure 15 we plot predicted SI growth rates under putative disk turbulent conditions, along with the largest value of St admitted by our modeling for several different values of α t . The SI growth rates (colors) are based on the theory presented in Umurhan et al. (2020) wherein turbulence is represented by an enhanced viscosity and particle diffusivity quantified by α t . In this figure the growth rate of the fastest growing SI mode is shown as a function of α t and St for a disk model with Z = 0.01 and with a local opening angle δ ≡ H/R = 0.07 representative of thermodynamic conditions appropriate to the snowline (T 160K; Umurhan et al. 2020, also see their Table 6 and Fig. 12) 12 . All St points are from our models, and selected at times that closely satisfy these metallicity and temperature conditions. As discussed at length in both Umurhan et al. (2020) and Chen & Lin (2020), the diagonal = 1 line (where = ρ solids /ρ g ≈ Z 1 + St/α t ) divides the parameter space regarding SI activity. For local conditions where falls below 1 (the region in Fig. 15 above the = 1 line and with St 0.3) the SI is severely weakened, rendering it essentially inoperative. For > 1 (below the = 1 line and/or St > 0.3), the SI can grow relatively fast (1-10 orbit times), and planetesimals readily form in numerical simulations conducted assuming these conditions, as discussed in Umurhan et al. (2020). The graphic also shows other lines representing several familiar particle growth obstacles: cold and warm ice fragmentation, particle radial drift, and SI mode pattern drift (the prescriptions for these constraints are detailed in Sections 7.2-7.4 of Umurhan et al. 2020). The plot also shows a horizontal line representing the likely minimum level of globally driven HD disk turbulence with α t > 4 × 10 −5 (e.g., the VSI, Flock et al. 2017;Lyra & Umurhan 2019), which is likely to be widespread throughout the disk except possibly for singular narrow radial patches in which no linear instability is thought to operate (e.g., see Figs. 1 and 10 of Pfeil & Klahr 2019).
As Fig. 15 demonstrates, most of the self-consistent combinations of particle St and turbulent intensity α t that result from our growth modeling fall within the < 1 zone of the α t -St parameter space in which the SI does not lead to planetesimal formation. Only particle growth models with α t = 10 −5 show the possibility of > 1 (and then for only a brief time, as discussed earlier in this section) 13 , where SI growth can lead to planetesimal formation. But, such low levels of global turbulence are currently thought to be unexpected because of currently known routes to linear instability in protoplanetary disks (section 1).

Observations of Particle Porosity
There is some tension between the low internal densities of the fractal aggregates in our simulations (Sec. 3.3, and as found by others, Okuzumi et al. 2012;Krijt et al. 2015) and the ∼ 0.1 − 0.5 g cm −3 internal densities found in comets and cometary particles (e.g., Sierks et al. 2015;Fulle et al. 2016;Pätzold et al. 2016). Another source of tension is a possible discrepancy between the high porosities of our particles and ALMA mm-cm wavelength observations of protoplanetary disks, which most analyses treat as compact. No systematic attempt has yet been made to actually constrain particle porosity from ALMA observations, although it would seem desirable to try based on the maturity of the field (Carrasco-González et al. 2019;Macías et al. 2021;Sierra et al. 2021). This is partly because the radiative transfer problem involved has recently become more challenging. First, there is a growing recognition of the importance of scattering by wavelength-sized particles, in addition to pure emission (Kataoka et al. 2014(Kataoka et al. , 2015Zhu et al. 2019), especially for large disk optical depths. Second, treating the scattering properties of highly porous particles with the current approximate techniques may be problematic in some ways (Tazaki & Tanaka 2018). Birnstiel et al. (2018) published calculations of extinction efficiencies for porous particles using a traditional set of approximations, but did not systematically compare them with observations along with results for compact particles, while allowing for possibly degenerate parameters (specifically, the also-unknown actual particle size), so the case remains open. A proper analysis of mm-cm observations for improved constraints on particle porosity would be quite valuable. Figure 16. Variation of the monomer size for model fa3gQ (cold H2O ice model) after 0.1 Myr. In the context of our adopted fractal growth model , larger monomer sizes do lead to slightly more compact particles compared to the nominal model (solid curve), but porosities are still much larger than typical cometary IDPs.
In addition to improving the radiative transfer modeling, there are several other ways one might reconcile the porous particle models and the observations. The particle growth theory used here uses the most advanced current treatment of the porosity of growing, colliding aggregates, but even this treatment is based on certain assumptions and some parameters are uncertain. One of these is the assumed monomer size which we have taken to be 0.1 µm in this work. In Figure 16, we show that increasing the monomer size to values as large as 2 µm (e.g., Okuzumi & Tazaki 2019) in our adopted fractal growth model (Suyama et al. 2012;Okuzumi et al. 2012;Kataoka et al. 2013a) makes some difference, but not nearly enough to bring porosities to values consistent with most cometary IDPs.
Perhaps a more straightforward way to reconcile the observations of actual disk particles with the models of truly primary growth presented here is to consider compaction processes that might occur after sizeable planetesimals form. Aggregates of sufficient mass can compact due to self-gravity, for example. Using Eqns. (39) and (40) one can estimate that for an icy aggregate with such low internal density to compact sufficiently to an internal density of 0.1 g cm −3 would require a mass of ∼ 5 × 10 18 g, equivalent to a ∼ 20 − 25 km radius body. Thus formation of planetesimals via whatever means to at least this mass could be a pathway to resolve this conundrum. The debris from subsequent breakup of such planetesimals (which may make up the dust content at later times) should then retain lower porosities. So, if the particles in disks studied to date are really low-porosity, it might imply they are secondary, rather than primary, products. Observations of earlier disks (Types 0/I) would also be revealing.

SUMMARY
In this work, we have conducted simulations that compare compact and fractal aggregate particle growth models in globally evolving early stage protoplanetary disks (ages up to 0.5 Myr, covering the formation age of the first planetesimals). We have studied most closely the effects of varying the turbulent intensity α t , but have also explored a subset of the potential parameter space for other nebula or particle growth properties. We can draw some general observations from our results.
1. Fractal aggregates grow to much larger masses than their compact particle counterparts, and reach those masses in a much shorter time. This is mostly because of the larger cross sections that characterize fluffy, underdense aggregates. Inside the snowline, growth is generally fragmentation limited across all models, with fragmentation masses being larger by up to several orders of magnitude for fractal particles than compact ones of the same mass, because of the smaller St (and lower collision speed) of the porous particles. Growth beyond the fragmentation barrier, due to mass transfer and a velocity PDF (Windmark et al. 2012a;Garaud et al. 2013), diminishes with decreasing turbulent intensity for both growth models, with no "lucky particles" seen for α t = 10 −5 in either case. This comes about because of a reduced production rate (by fragmentation) of particle "feedstock" at the smaller sizes (m < m b ) suitable for sweepup by larger particles, so that beyond the fragmentation limit the mutual collisional destruction rate dominates over growth by sweepup. Especially beyond the snowline, where the fragmentation mass limit is much higher (with the exception of "cold ice" model fa3Qg), growth rarely reaches the fragmentation mass m f , except at the highest α t compact growth model (sa2g), where fragmentation can once again become an effective production source of small particles to feed "lucky" large ones. Overall, the turbulent intensity α t is one of the most important factors determining particle growth outcomes and thus global nebula radial structure.
2. Our simulations confirm that radial drift is a less influential factor for fractal aggregates than for compact particles, for a given α t , during much of their growth phase. The Stokes numbers for small compact particles initially increase more quickly compared to fractal aggregates (despite the faster growth rate of fractal particles overall) so porous particles have slower radial drift at early times. By 0.1 Myr or less though, in almost all models, fractal aggregates just outside the snowline where growth is most rapid achieve larger Stokes numbers than the corresponding compact particles and begin to suffer strong radial drift. Meanwhile, beyond ∼ 20 AU the Stokes numbers of porous aggregates typically remain smaller by up to an order of magnitude relative to compact particles even after 0.5 Myr. Thus, a characteristic of fractal growth models is that they retain significant amounts of mass in the outer disk for longer periods, while carving out broad, deep gaps in particle surface mass density from the snowline to tens of AU. Interestingly though, as the turbulent intensity increases to higher values (α t 10 −2 ), this distinction between fractal and compact growth models begins to lessen considerably.
3. Despite their larger masses (and sizes) and longer retention in the outer disk, we conclude that for our given compaction model and its assumptions, even porous fractal aggregates will also eventually be lost en masse to the inner disk and the central star due to radial drift, without the intervention of a process or processes of planetesimal formation that would allow solids to be retained for much longer. However, the longer retention time for the fractal aggregates would allow more time for planetesimals to potentially form, mitigating material loss compared to compact particle growth models.
4. We find that even after collisional compaction sets in, for either Gaussian or Maxwellian (Appendix C) collision velocity PDFs, it generally does not lead to increasing internal density as aggregates grow in mass, in agreement with previous work (e.g., Okuzumi et al. 2012;Krijt et al. 2015). Compaction is instead dominated by gas ram pressure, with internal particle density increasing along with increasing relative particle-to-gas velocity ρ int ∝ (∆v pg ) 1/3 (Sec. 2.4.5). Self-gravitational compaction is not achieved in our models because the largest aggregates fall just short of the required critical mass.
5. We tested models in which the fragmentation energy threshold of ice-rich particles Q f is significantly larger than silicates (as is usually assumed) only over a limited range of temperature about the snowline (Musiolik & Wurm 2019), which we refer to as cold H 2 O ice models (e.g., see Umurhan et al. 2020). These models (which also apply to supervolatile ices) retain material outside the snowline for much longer periods, because the icy aggregates and their Stokes numbers are then an order of magnitude smaller than under the usual "strong water ice" assumptions throughout the outer nebula, except for in a narrow band just outside the snowline. As a result, inward drift transport of condensibles is diminished, and both the depletion of solids outside the snowline, and the enhancement of solids Z inside the snowline, are muted (see Appendix B, and Paper III).
6. Simulations in which we vary the initial disk mass are qualitatively similar in how they evolve compared to the fiducial model, with predictably lower disk temperatures for decreasing disk mass. The more extended disk model R 0 = 60 AU (Appendix B) maintains a higher metallicity Z in the planet forming region (and out to ∼ 200 AU) because the model starts with more mass distributed in the outer disk. Differences in St are most apparent outside of ∼ 20 AU, but inside 20 AU, Stokes numbers are remarkably similar.
7. Evaporation fronts serve as locations where solids can become trapped over long periods, even after the bulk of solid material has largely been lost; this retention effect appears especially dramatic for α t ≤ 10 −4 . The effect is characterized by narrow bands of volatiles, each lying stably just outside its EF. In the case of supervolatiles, the bands could continue to be fed by outwardly advecting and/or diffusing gas, that is partly composed of the supervolatile's vapor phase, recondensing outside the respective EF (see Paper III, for more discussion). It is intriguing to compare these structures with narrow bands seen in multiple ALMA images of disks (e.g., Andrews et al. 2018), especially Class 0/I (e.g., Segura-Cox et al. 2020). This effect may be amplified by decomposition of refractory organics at their "EF", and conversion into simpler supervolatiles.
8. We find that particles in our disks beyond R 200 − 300 AU (where Σ drops off sharply, and the gas radial velocity v g increases rapidly) have increasingly large St so that the combination of radial drift and small ρ solids means growth times are exceedingly long, at least until the nebula gas expands outwards. However, this may be an artifact of our initial conditions, and the effect is lessened with larger R 0 or α t . 9. Quite generally, our simulations do not achieve conditions necessary for planetesimal formation by streaming instability (SI), with almost all models resulting in combinations of St, α t (and Z) that place them firmly in the "turbulent regime" where the local particle-to-gas mass volume density ratio = ρ solids /ρ < 1 (referred to as Zone II, Umurhan et al. 2020); indeed this basic criterion has been a known requirement for SI since the earliest studies (Youdin & Goodman 2005). We find that only for the lowest turbulent intensity (α t = 10 −5 ) models is there a relatively short period of time lasting ∼ 2 × 10 3 − 2 × 10 4 years, and only at the snow line, where > 1 (laminar regime, Zone I, Umurhan et al. 2020). This α t and corresponding Stokes numbers (∼ 0.05 and ∼ 0.4 for sa5g and fa5g, respectively) are consistent with critical values necessary for SI to lead directly to planetesimal formation as is usually assumed (Umurhan et al. 2020). However, for the higher turbulent intensities of α t ∼ 10 −4 − 10 −3 commonly found in simulations of critically triggered hydrodynamic turbulence (for a review, see Lyra & Umurhan 2019;Lesur et al. 2022), this moment of instability is never achieved, because particle growth to large enough sizes (masses) is not fast enough to keep up with the radial drift of material across the snowline. The local Z drops precipitously to such low values that the mechanism never becomes viable and/or the particles are too large and rare to drive the gas velocities as required for SI.
We thank the NASA Ames Research Center's Origins Group for numerous conversations, and specifically useful discussions with Denbanjan Sengupta, Uma Gorti and Maxime Ruaud as well as Carlos Carrasco-Gonzalez, Enrique Macias, and Anibal Sierra. We thank Thomas Hartlep and Diane Wooden for thorough internal reviews, and we thank an anonymous reviewer for very helpful comments to improve the paper's exposition. P.R.E. and J.  Figure A1 shows as a reference the particle size distributions for compact particle and fractal aggregate growth models with α t = 10 −4 , 10 −3 and 10 −2 after 0.5 Myr as a function of semi-major axis and solids mass volume density m 2 f (m). In these panels are also highlighted the locations of various EFs in the inner 10 AU of the disk after this simulation time. These are the water snowline (cyan), the organics EF (orange) and the troilite EF (green). As noted in the text, other EFs were present at earlier times but have since migrated beyond the computational grid inside 0.5 AU as the disk cools. "Flat" regions seen in these distributions correspond to particle sizes that are below the fragmentation size r f . Recall that in our model the "dust" distribution (defined by r < r f ) is modeled using the moments method (see Sec. 2.4, and Paper I), while particle growth is followed explicitly for r ≥ r f . The non-flat regions can thus be associated with "lucky" particle growth beyond the fragmentation barrier seen in the related mass distribution plots for the given α t (section 3).
The leftmost panels show snapshots for the sa2g (top) and fa2g (bottom) models presented in Sec. 3.2.1. For α t = 10 −2 , particle sizes are smallest compared to lower turbulence models because the fragmentation and bouncing barriers occur at much smaller sizes. For instance, r f ∼ 0.001 cm at 1 AU for sa2g, whereas it is ∼ 1 cm for fa2g. As a result, there is substantial growth beyond the fragmentation barrier in the inner disk for both models, and for the compact case, also outside the water snowline. This figure also shows clearly the secondary and even tertiary growth peaks that were seen in Fig. 7 (rightmost panel of top row), which are due to a combination of mass transfer and our velocity PDF (see Sec. 3.1). Model fa2g also has additional peaks beyond the secondary, more easily seen here than in  Figure A1. Particle size distributions for compact particle (top panels) and fractal aggregate growth models (bottom panels) plotted on roughly the same scale in solids mass volume density m 2 f (m) after 0.5 Myr for the indicated simulations. These wireframes (plotted out to 10 AU) can be compared with the rightmost panels of the particle mass distributions in the main text. Several EFs are indicated by the cyan (water), orange (organics) and green (troilite) curves. Regions that appear flat at this scale correspond to particle sizes that have not reached the fragmentation barrier. Recall that in the moments method, the "dust" population is modeled as a power law up to the fragmentation barrier r f . the corresponding particle mass distribution plot. The tertiary peaks decrease in magnitude, however, as particles are ground down -likely due to increasing destruction probabilities further inward in the disk. The largest particle sizes achieved are roughly 1 cm for the compact particle growth model, but as much as 2 − 3 m outside the water snowline in the fractal aggregate model. Note that even the secondary mass peaks contain far less mass than at r f (see also Paper I).
The middle panels show snapshots of our fiducial models sa3g (top) and fa3g (bottom), with α t = 10 −3 . These cases display a more stark contrast between fractal aggregate and compact growth models regarding lucky particles. The compact particle case still shows significant growth beyond the fragmentation barrier out to the water snowline with a secondary and tertiary growth peak (both clearly visible in Fig. 2, right upper panel), while in the fractal aggregate case no secondary peak develops and the mass and size of the mass-dominant aggregate (m M , r M ) remain near, but are slightly larger than, the values for fragmentation (black curve, lower right panel of Fig. 2). Due to the lower α t , the fragmentation sizes are considerably larger than the α t = 10 −2 model. At 1 AU, the fragmentation size for model sa3g is ∼ 0.03 cm, whereas it is ∼ 13 cm for model fa3g. The largest particle sizes achieved outside the water snowline are ∼ 15 cm for the compact particle case, but ∼ 27 m for the fractal aggregate model. Though more difficult to see in this orientation, a systematic decrease in the solid mass volume densities can be seen outside the water ice EF.
Finally, the rightmost panels show snapshots for the sa4g (top) and fa4g (bottom) models shown in Fig. 8 of Sec. 3.2.2. In this case there are no secondary growth peaks. In the compact particle model, growth beyond the fragmentation size is very limited and the mass dominant particle is still near the fragmentation barrier (upper right panel, Fig. 8), while lucky particle growth is completely stymied in the fractal aggregate model (see discussion, Sec. 3.2.2). Aggregates that grow beyond r f find only relatively brief existence before being destroyed. As expected with such a low α t , the fragmentation sizes and largest particle sizes achieved are larger than in the models discussed above. At 1 AU, the fragmentation size in model sa4g is ∼ 2.3 cm, and is ∼ 54 cm in model fa4g. Outside the water snowline, the largest sizes achieved are ∼ 38 cm and ∼ 64 m, respectively. As in the fiducial case, a systematic drop in the solid Figure B1. Fiducial model after 0.5 Myr (black solid curves) compared to a suite of models where we vary the disk mass MD, disk scaling parameter R0, and fragmentation threshold Q f as listed in Table 2. (a) Gas surface densities; (b) total solids surface densities; (c) midplane temperatures; (d) disk metallicities; (e) Rosseland mean opacities; and (f) Stokes numbers of the mass dominant aggregates, all as a function of semi-major axis. The solid blue curves that appear in panels (a-c) are for a standard MMSN. In each panel, the cyan region demarcates the range in location of the water snowline for these models (between ∼ 2.7 − 4 AU). mass volume densities can be seen outside the H 2 O EF that can be associated with the steep drop in solids surface density (Fig. 9).

B. VARIATION WITH SOME INITIAL DISK PARAMETERS
In this appendix, we compare our fiducial model for fractal aggregate growth with α t = 10 −3 (fa3g, Sec. 3.1) to simulations with different disk mass M D , disk size scaling parameter R 0 , and fragmentation threshold Q f to assess the effects on the particle mass distribution. These simulations vary only a subset of parameters, and of the parameters highlighted only a fraction of their parameter space is explored here. As summarized in Table 2, specifically we compare models with initial disk masses of 0.1M (fa3M01g) and 0.02M (fa3M002g), a disk scaling parameter R 0 = 60 AU (fa3R060g), and perhaps of greatest interest a model in which water ice is only "sticky" in a narrow range about the snowline (fa3Qg) as suggested from the recent work of Musiolik & Wurm (2019). This model is referred to as the "cold ice" model in Umurhan et al. (2020). The lowest mass model fa3M002g corresponds most closely with the standard MMSN 14 which has surface density ∝ R −3/2 , and temperature ∝ R −1/2 , while the gas surface densities in the higher mass models (including our fiducial model) though more massive, are generally in agreement with the models of Desch (2007) in radial extent hypothesized for the same timeframe. In Figure B1 (panels a-f) we plot a subset of modeled quantities at 0.5 Myr as a function of semi-major axis for the suite of fractal aggregate growth simulations described above. In all panels our fiducial model (Sec. 3.1) is given by the solid black curve. The solid blue curves that appear in panels a-c refer to a standard MMSN disk.
Focusing first on models in which the disk mass is varied (fa3M01g, fa3M002g), the differences between these and the fiducial model are intuitive -Σ is lower overall and evolves similarly to the fiducial case, although the lower mass model is much cooler in the inner disk (dotted curve in panel c, closer to MMSN). The variations in Σ as a function of semi-major axis seen between 1 − 4 AU in the higher mass models are due to strong gradients in T , associated with the water and organics snowlines, that lead to fluctuations in the local viscosity ν and thus relative velocity. This contributes (along with the smaller fragmentation threshold inside the snowline) to the sharp decrease in particle St (panel f) near ∼ 3 AU, which then coincides with a pileup in Σ solids (panel b) and thus an enhanced Z (panel d). These effects are muted in the lowest mass model (fa3M002g; dotted curve, panel c) at this point in the simulation because, owing to its initially low mass, this model cools quickly. The lower T in fa3M002g means that the fragmentation limit St f , which is in the Stokes regime (Re p < 1), is higher than for the fiducial case by a factor of several in the inner disk (recall from Eq. 27 that St f ∝ T −1 ), and leads to generally higher Stokes numbers (panel f) so that solid material drains away at a faster rate. Beyond the snowline on the other hand, particle masses in this low mass disk are also smaller than in the fiducial model fa3g, leaving St comparable (or slightly smaller), so particles experience similar radial drift speeds. Although the evolution of Σ solids appears to follow the evolutionary trend of higher mass models, as can be seen in panel e, Z is generally smaller everywhere compared to the fiducial model. In the inner nebula, model fa3M01g is also cooler than the fiducial case, but less so than fa3M002g so that differences with the fiducial model are less distinctive.
Model fa3R060g (in which the disk scaling parameter is set to 60 AU rather than 20 AU) initially distributes the mass such that the surface density is smaller in the inner disk (∼ 2 − 2.5 times less at 0.5 AU), and larger in the outer disk (∼ 10 times higher at 100 AU). This leads to Z (panel d) having its steep cutoff much further out in the disk ∼ 300 AU. In other models, this region is where the Stokes numbers begin to increase dramatically -causing rapid depletion of material within the region, and isolation of material outside of it where St 1 (e.g., right panel of Fig. 3, and Fig. 6, bottom panel). This behavior still occurs in this model at early times, but after 0.5 Myr, Stokes numbers in the outermost portion of the disk are mostly St 1, leading to relatively rapid radial drift with very little material left isolated (see Paper III). Generally the temperature and Stokes numbers are comparable to the fiducial model, but one notable difference is that in fa3R060g there is more material outside the snowline (which is reflected in the plot of Z in Fig. B1), so the opacity there is higher as well (panel e). This is simply because the disk starts out with much more mass at larger radial distance. The more extended disk retains more mass for the same evolutionary time, and in fact there is about ∼ 2.5 times more total mass remaining in this model compared to the fiducial case (Paper III).
Finally, we look at a cold H 2 O ice model in which Q f is much larger than for silicates only in a region about the snowline (fa3Qg, red solid curves), compared to all other models in which water ice is sticky everywhere it is solidthat is, everywhere outside the water ice EF. The surface density of solids outside the snowline (located at 4 AU in this and the fiducial model) does not decrease nearly as much as the fiducial model (which is factor ∼ 2.5 times lower at ∼ 5 AU) because aggregates are considerably smaller for this Q f (they are fragmentation limited), so their radial drift times are longer. Interior to the snowline, the surface density becomes lower than in the fiducial case because less material has been able to drift inwards across the snowline, so aggregates (and their St) are smaller there as well. The decrease in Σ solids outside the water snowline still correlates with the substantial decrease (minimum at ∼ 5 AU) in the opacity, but outside and just inside the opacity is higher which leads to the modest increase in temperatures seen about the snowline in panel (c). The clearest indication of the effects of the lower Q f across most of the outer nebula, however, can be seen in the distribution of Stokes numbers with semi-major axis (panel f) where only a very narrow region near the snowline displays the peak values of other models (St ∼ 0.03 − 0.04), while outside 5 − 6 AU the Stokes numbers drop to ∼ 10 −3 and then increase slightly going outwards. The behavior is similar to a constant Q f = 10 3 model in Paper I. The spike in St just outside the snowline in panel (f) still coincides with a dip in Z (panel d), though the drop is less dramatic than in the fiducial model. Like every other fractal model, this model also retains more mass in the outer disk, but because of the low Q f , the surface density of solids (and Z) remains considerably higher from ∼ 6 − 70 AU (panels b and d).

C. EFFECTS OF THE CHOICE OF VELOCITY PDF
Here, we examine the effect of the velocity PDF on growth beyond the fragmentation barrier. In Section 2.4.3 (also, see Sec. 2.4.1) we discussed the destruction probabilities we utilize to model growth beyond the fragmentation barrier which, combined with mass transfer, determine the largest mass a compact or aggregate particle can grow to beyond m f . We found in our simulations that the dominant limitation in growing beyond the fragmentation barrier is the mass m f at which the barrier is reached, which determines the amount of "dust" remaining for lucky particles with m > m f to grow from. The value of m f is largest for the smallest α t , but also is larger for fractal aggregates than for compact particles at the same α t (e.g., see Sec. 3.1). Thus, models with α t 10 −4 saw little if any growth beyond m f , while α t = 10 −2 saw the most because there is still ample feedstock (especially m < m b ) to grow from (see also Appendix A). Despite the difference in fragmentation mass, the largest particle mass m L achieved inside the snowline, and its variation with semi-major axis, are similar for most all of our simulations, roughly ranging from ∼ 10 − 1000 g. An obvious exception is the lowest α t = 10 −5 fractal aggregate model (Fig. 10, lower left panel). Note that when m L > m f (i.e., lucky particles), they generally carry a negligible amount of the mass, whereas outside the snowline, aggregates and particles generally do not reach or exceed m f and so m L is the mass dominant particle.
To compare directly with our Gaussian PDF approach, we also employ the model of Garaud et al. (2013) which more rigorously takes into account that particles both have deterministic and stochastic velocities. In our notation, where ∆v D (m , m) is the deterministic relative velocity as defined in Sec. 2.4.1 between m and m , and which ensures that the mean velocity of the distribution is what would be obtained when the stochastic velocities dominate the motions. The integrand in Eq. (C1) reduces to a Maxwellian distribution in the limit that ∆v D σ, whereas in the opposite limit when ∆v D σ, it reduces to a Gaussian with a significant tail for large values of the velocity. In the latter limit, the deterministic velocities are much larger than their rms velocities. Integration of Eq. (C1) gives In the limit where ∆v D is very small compared to σ (e.g., for same-sized particles), the 3rd term on the RHS of Eq. (C3) reduces to 2/π(v c /σ)exp[−v 2 c /2σ 2 ]. As noted by Garaud et al. (2013), the mean collisional velocity is not the same as the traditional construct (Eq. which we adopt for simulations that utilize the Garaud et al. scheme. Figure C1 shows a comparison between our Gaussian PDF for the relative velocities and destruction probabilities between particles (Eqns. 22 and 29, and designated by "g"), or a Maxwellian PDF (Eqns. C4 and C3, labeled as "m") for several simulations as indicated in Table 2. The Maxwellian approach is more rigorous in the limit that systematic Figure C1. Comparison of simulations that use Gaussian (g; Sec. 2.4.3) and Maxwellian (m) velocity PDFs . These highlight subtle differences in the particle mass distributions, and in particular the difference in the mass of the "lucky particles" mL, which are the largest mass to which particles or aggregates can grow beyond the fragmentation barrier m f (dashed curves). The solid curves are the dominant mass particle mM of the distribution. Sharp variations in mM and mL especially seen in the top two panels are associated with crossing EFs. For example, the organics EF is at ∼ 1.5 AU, and the water-ice EF is located at ∼ 4 AU, outside of which Q f increases so aggregates can begin to grow larger. velocities ∆v D ∆v t whereas in the opposite limit, the PDF is Gaussian-like but with a long tail. Our standard approach used here (and also in Paper I) is Gaussian in both limits, thus the possibility of significant differences arises at least in the limit where (stochastic) turbulence-induced velocities dominate. The top panel of Fig. C1 shows that the difference between the two fractal aggregate growth models with our highest α t (=10 −2 ) simulations after 0.5 Myr is not considerable. Most notable is that the mass m L of the largest or luckiest particle tends to be smaller when the PDF is Maxwellian. The corresponding mass of the mass-dominant particle m M appears smaller for the Maxwellian PDF in some cases, or comparable in others.
A more direct PDF comparison between compact and fractal growth models, for the fiducial case α t = 10 −3 , can be seen in the lower two panels of Fig. C1. The compact growth models show the most notable differences in m L and m M (the Maxwellian values again tending to be smaller) and the fractal growth simulations show smaller differences (bottom panel). We note that the difference in the fragmentation mass between models sa3g and sa3m inside of 1 AU is due to modest differences in the temperature profiles, while somewhat sharp variations in m M and m L are due to crossing EFs, notably the organics EF at ∼ 1.5 AU in the top two panels. These PDF differences effectively vanish for smaller α t as the turbulence-induced velocities become less effective on the more massive particles or aggregates at the fragmentation limit, and growth beyond the barrier mostly stalls. Thus we find, at least for these simulation times and disk parameters, that between the two PDFs, the differences are not significant enough to affect the qualitative nature of the results. Execution of the Garaud et al. scheme in our code is slightly more involved, however, which led to roughly 20% longer simulation times for the same parameter set.