Massive Protostellar Disks as a Hot Laboratory of Silicate Grain Evolution

Typical accretion disks around massive protostars are hot enough for water ice to sublimate. We here propose to utilize the massive protostellar disks for investigating the collisional evolution of silicate grains with no ice mantle, which is an essential process for the formation of rocky planetesimals in protoplanetary disks around lower-mass stars. We for the first time develop a model of massive protostellar disks that includes the coagulation, fragmentation, and radial drift of dust. We show that the maximum grain size in the disks is limited by collisional fragmentation rather than by radial drift. We derive analytic formulas that produce the radial distribution of the maximum grain size and dust surface density in the steady state. Applying the analytic formulas to the massive protostellar disk of GGD27-MM1, where the grain size is constrained from a millimeter polarimetric observation, we infer that the silicate grains in this disk fragment at collision velocities above ~ 10 m/s. The inferred fragmentation threshold velocity is lower than the maximum grain collision velocity in typical protoplanetary disks around low-mass stars, implying that coagulation alone may not lead to the formation of rocky planetesimals in those disks. With future measurements of grain sizes in massive protostellar disks, our model will provide more robust constraints on the sticking property of silicate grains.


INTRODUCTION
Understanding the formation of kilometer-sized planetesimals from (sub)micron-sized grains is a critical issue in planet formation. It is generally accepted that the grains grow into macroscopic aggregates through coagulation (Johansen et al. 2014;Drazkowska et al. 2022). However, whether the macroscopic aggregates can continue growing into kilometer-sized bodies by coagulation alone is much more uncertain because of their rapid radial drift (Weidenschilling 1977) and finite stickiness (Blum & Münch 1993). If the grains are sticky enough, they may form fluffy aggregates and then grow into planetesimals more rapidly than they drift Kataoka et al. 2013). Otherwise, other mechanisms that rapidly concentrate dust, such as the streaming and gravitational instabilities (e.g., Goldreich & Ward 1973;Youdin & Goodman 2005;Johansen et al. 2007;Youdin 2011), are necessary for planetesimal formation. Therefore, in order to understand how planetesimals form, one must know how sticky the grains are.
A major concern in planetesimal formation studies is that the stickiness of silicate grains is highly uncertain. Silicates are the ultimate building blocks of rocky bodies, and therefore knowing their stickiness is necessary for understanding how rocky planets like the Earth form. Early grain sticking models (Chokshi et al. 1993;Dominik & Tielens 1997) predicted that micronsized rocky grains can stick at collision velocities up to ∼ 0.1 m s −1 . However, laboratory experiments typically show that micron-sized silica grains can stick even at 1 m s −1 (Poppe et al. 2000;Blum & Wurm 2000). The source of the large discrepancy is yet to be confirmed, but recent studies (Kimura et al. 2015;Steinpilz et al. 2019) suggest that the early models underes-arXiv:2303.09148v2 [astro-ph.EP] 7 Jun 2023 timated the surface energy (a measure of the strength of the intermolecular forces) of silicates. The stickiness of dust aggregates also depends on the size of the constituent grains (Chokshi et al. 1993;Dominik & Tielens 1997). For instance, aggregates made of 0.1 µm-sized interstellar silicate grains may stick at velocities up to 50 m s −1 (Kimura et al. 2015). Because of this large uncertainty, one cannot yet conclude whether coagulation of silicate aggregates into kilometer-sized bodies is possible or not in protoplanetary disks, where the maximum grain collisional velocity is ≈ 20-50 m s −1 unless strong turbulence is present (e.g., Figure 4 of Johansen et al. 2014).
Given this importance, it is desirable to constrain the stickiness of silicate grains from observations of protoplanetary disks. However, it is challenging to do this because the disks around low-mass stars (∼ 1 M , 1 L ) are cold except in the innermost region or in the accretion outburst phase (e.g., Liu et al. 2021). In these disks, the radius of the snowline, where water ice sublimates, is only a few au or even smaller (e.g., Mori et al. 2021), and it is difficult to spatially resolve their silicatedust regions inside the snowlines. On the other hand, the snowline radii around massive protostars ( 10 M , 10 3 L ) are typically much large as several hundred au or more due to the intense stellar radiation and high accretion heating. Thus, silicate-dust regions in massive protostellar disks can be resolved by current telescopes, such as Atacama Large Millimeter/submillimeter Array (ALMA), although massive star-forming regions are typically distant as several kpc. In particular, the 1.14-mm polarimetric ALMA observations by Girart et al. (2018) found that the maximum grain size is as large as several hundred microns in the massive protostellar disk of GGD27-MM1. Given the high luminosity and the accretion rate of GGD27-MM1 (∼ 10 4 L and 10 −4 M yr −1 , Añez-López et al. 2020), its water snowline is likely to present outside of the accretion disk, suggesting that the observed large grain size is an excellent indicator of the coagulation of silicate dust.
Motivated by the detection of large-size silicate grains in GGD27-MM1, we develop the first theoretical model of dust coagulation in massive protostellar disks in this paper. Massive protostellar disks have not gotten significant attention as the site of dust coagulation, maybe because Earth-like habitable planets would not form around massive stars due to their intense radiation field and short lifetimes. However, massive protostellar disks are excellent observational targets to investigate the coagulation of silicate grains. Furthermore, ALMA observations have started to discover accretion disks around massive protostars in recent years (e.g., Ginsburg et al. 2018;Maud et al. 2019;Motogi et al. 2019;Tanaka et al. 2020;Johnston et al. 2020). We need to know the appropriate dust size (i.e., its opacity) to unravel their physical properties from observations, especially dust mass. Massive protostars have orders of magnitude higher luminosities (∼ 10 3 L ) and accretion rates (∼ 10 −3 -10 −4 M yr −1 ) than low-mass protostars (Hirota 2018). Thus, dust coagulation in massive protostellar disks is nontrivial: how much dust grains can grow, what process limits their coagulation, and how dust coagulation impact disk physical conditions. A theoretical model is required to understand observation results and constrain the sticking property of silicate grains.
The structure of this paper is as follows. We describe our disk model with dust coagulation in Section 2. We derive the analytic formula of the size limit of dust grains in Section 3, and show the results of numerical calculations in Section 4. We provide the discussion in 5, where we constrain the threshold fragmentation velocity of silicate dust by comparing our model with the observed grain size of GGD27-MM1. In Section 6, we present the summary of this paper.

MODEL
We present our model for massive protostellar disks with dust coagulation, which is schematically shown in Figure 1. We consider an accretion disk with a supply of gas and dust from the surrounding envelope. The disk's interior temperature is calculated by considering both stellar irradiation and accretion heating. For dust grains, we take into account collisional growth and fragmentation as well as radial drift and vertical settling. We detail each process in the following subsections.

Gas Disk
We assume that the mass flux from the envelope to the disk is approximately constant in time and axisymmetric. We then decompose the disk structure into a steady, axisymmetric part and time-dependent, non-axisymmetric disturbances. The latter may include density waves produced by the disk's self-gravity and turbulence. Below, we focus on the former structure and model the disk as being steady and radially one-dimensional. However, we also account for (timeaveraged) angular momentum transport induced by the disk's self-gravity.
The disk's rotation is assumed to be nearly Keplerian. In fact, the gas rotation is slightly slower than Keplerian owing to the outward pressure support. This small deviation is usually negligible regarding the gas disk structure. However, the sub-Keplerian motion must be taken into account when computing the radial motion of the grains (see Section 2.3). The disk is in a steady state with a supply of gas and dust from the envelope with the rates ofṀgas andṀ dust . The heating from the stellar irradiation and the accretion heating are considered. For the accretion heating, we take into account the vertical optical depth τ , which depends on the dust coagulation. We use the α viscosity to treat the angular momentum transfer of the gas disk. The α value increases where the disk is marginally unstable with the Toomre parameter of Q < 1.4. (right) For dust evolution, we consider dust collisional growth and fragmentation, as well as radial drift and vertical settling. The dust coagulation is suppressed if the collisional velocity is as high as or higher than the fragmentation velocity, i.e., ∆v v frag . The entire disk is the silicate-dust region inside the water snowline. The numerical calculation is executed from the disk outer radius r d to the inner radius of 10 au, or the dust sublimation radius with 1500 K.
The equation of continuity for an axisymmetric gas disk reads where Σ gas and v r,gas are the gas surface density and velocity, and r is the distance from the central protostar.
Assuming steady accretion, we obtaiṅ M gas ≡ 2πr|v r,gas |Σ gas = const., whereṀ gas is the gas accretion rate. We use the α viscosity to describe the gas angular momentum transport within the disk (Shakura & Sunyaev 1973). The accretion velocity in the steady viscous disk can then be written as where α is the dimensionless viscosity parameter and are the sound speed and Keplerian velocity, respectively, with M , G, k B , and m gas being the stellar mass, gravitational constant, Boltzmann constant, and mean molecular mass. The disk's vertically averaged temperature T is approximated by the midplane temperature. We use m gas = 3.9×10 −24 g, assuming a mean molecular weight of 2.34. A massive protostellar disk can be gravitationally unstable. In such a disk, the gravitational instability (GI) induces spiral density waves, whose torque transfer the disk's angular momentum. We therefore decompose the viscosity parameter α into two components: where α GI and α floor represent the contributions from the GI and any other angular-momentum transport mechanisms, respectively. Following Zhu et al. (2010), we model the former as where is the Toomre parameter for the GI (Toomre 1964) with Ω K = v K /r being the Keplerian angular velocity. Equation (7) mimics the efficient angular momentum transport by spiral density waves that develop when the disk is marginally unstable state of Q 1.4 (e.g., Boley et al. 2006;Zhu et al. 2010). Even in gravitationally stable disks (Q > 1.4), angular momentum transport by hydrodynamical and magnetohydrodynamical turbulence and by large-scale magnetic fields can still occur (see Lesur et al. (2022) for a recent review of potential angular momentum transport mechanisms in protoplanetary disks). Following Tanaka & Omukai (2014), we take α floor = 0.01.
To determine the disk temperature T , we assume radiative equilibrium and consider heating by stellar irradiation and disk accretion. We approximate T as where T irr and T acc are the temperatures in the limits where irradiation heating and accretion heating dominate, respectively. The former can be estimated as (Chiang & Goldreich 1997;Kusaka et al. 1970) where L is the stellar luminosity. Equation (10) assumes that the disk is optically thick to the starlight, which is valid unless small dust grains are heavily depleted. With vertically uniform viscosity, the accretiondominated temperature is (Hubeny 1990) where σ SB is the Stefan-Boltzmann constant, and τ /2 is the disk's vertical optical depth to its own thermal emission, measured from infinity to the midplane. An optically thicker disk is hotter because of inefficient cooling.

Optical Depth
Dust grains are the dominant source of the disk opacity. For grains with size distribution, the disk's vertical optical depth can be written as where σ abs (a) is the monochromatic absorption crosssection of a grain with the radius a, dN dust (a, r)/da is the number surface density of grains per unit grain radius, and a min and a max are the minimum and maximum grain radii, respectively. The grain size distribution is assumed to follow a power law Here, the normalization factor C satisfies where Σ dust is the dust surface density, and ρ int = 3.0 g cm −3 is the internal density of silicate grains. Equation (13) well approximates the high-mass end of the grain size distribution in coagulation-fragmentation equilibrium in turbulent disks (Birnstiel et al. 2012). We take a min = 0.1 µm because grains smaller than this size generally grow quickly thanks to Brownian motion (Birnstiel et al. 2011). The maximum grain size a max is determine from a dust coagulation model (Section 2.4).
For σ abs , we use a crude approximation (e.g., Ivezic et al. 1997;Ormel 2014;Fukuhara et al. 2021), where λ peak = 10(300 K/T ) µm is the peak wavelength of the Planck function with T . As long as T 100 K, Equation (15) approximates the real grain opacity to within a factor of a few (Dullemond et al. 2022). For convenience in presenting the results, we introduce the absorption opacity per gas mass as κ abs = τ /Σ gas .

Dust Surface Density
We employ the single-size approach of Sato et al. (2016) to compute the dust evolution in the disk. This approach assumes that grains of one characteristic size dominate the dust mass surface density at each orbit r. For the grain size distribution given by Equation (13), the mass-dominating grains are those with a ∼ a max 1 . In our model, we identify the mass-dominating grains with the largest grains.
In the single-size approach, the dust surface density Σ dust in an axisymmetric disk generally obeys where v r,dust is the radial velocity of the largest grains. Because the radial flow of the dust from the envelope to the star is assumed to be steady, Equation (17) yieldṡ whereṀ dust is the dust accretion rate. We takeṀ dust to beṀ gas times the interstellar dust-to-gas mass ratio For dN dust /da given by Equation (13) with amax a min , one has ap ≈ 0.52amax. of f dg,ISM = 0.01. We note that the ratio ofṀ gas anḋ M dust is constant in the entire disk (Equations 2 and 18), but the ratio of Σ gas and Σ dust is not, because the radial velocities of gas and dust are different due to the dust radial drift.
The dust radial velocity is determined by gas drag in the disk. We use the expression (Takeuchi & Lin 2002 where St is the largest grains' stopping time (see below) normalized by the inverse Keplerian frequency, i.e., St = Ω K t s , and is a dimensionless parameter characterizing the disk's sub-Keplerian motion (Adachi et al. 1976), with P being the gas pressure. Note that η is positive in a smooth gas disk with radially decreasing pressure. The parameter St, called the Stokes number, measures how strongly the grains' motion is coupled to the gas motion, with the limit St → 0 corresponding to the extreme case where the grains and gas move at the same velocity. On the right-hand side of Equation (19), the first term simply represents the radial motion induced by gas accretion. The second term represents the so-called radial inward drift in a sub-Keplerian rotating disk: grains feel the headwind of the background gas, lose angular momentum, and fall toward the central star (Whipple 1972;Adachi et al. 1976;Weidenschilling 1977). For simplicity, we fix dlnP/dlnr = −24/7, which exactly applies when T = T irr and Q = const. (Tsukamoto et al. 2017). The stopping time characterizes the timescale on which the grains' velocity relaxes into the terminal velocity under gas drag. If the largest grains are smaller than the gas molecules' mean free path, they obey Epstein's law where ρ gas is the gas density, and v th = 8/πc s is the molecules' mean thermal speed. Following Birnstiel et al. (2010), we approximate ρ gas with the midplane density in an isothermal disk, ρ gas = Σ gas /( √ 2πH gas ), where H gas = c s /Ω K is the gas scale height. We then have In our calculations, most of the largest grains are smaller than the molecules' mean free path of λ mfp = m gas /(σ mol ρ gas ), where σ mol = 2.0 × 10 −15 cm −2 is the collision cross-section of the gas molecules. Therefore, the use of Epstein's law is valid.

Maximum Grain Size
We take into account collisional fragmentation depending on the grains' collision velocity. The collisional growth/fragmentation and radial drift control the maximum grain size at each disk radius.
In the single-size approach, the mass of the largest grains, m max = (4π/3)ρ int a 3 max , obeys (Sato et al. 2016;Okuzumi et al. 2016) where t coll is the largest grains' mean collision time, and ξ frag is the fractional change of m max upon a single collision, accounting for collisional fragmentation. Because we consider steady state, Equation (23) reduces to an ordinary differential equation for a max , Following Okuzumi & Hirose (2012) and Okuzumi et al. (2016), we model ξ frag as where ∆v is the mean collision velocity and v frag is the threshold velocity above which a net mass loss of the largest grains occurs. Hereafter, we simply call v frag the fragmentation velocity. Equation (25) is based on dust aggregate collision simulations by Wada et al. (2009), which show that ξ frag is ≈ 1 at ∆v 0.2v frag and decreases approximately logarithmically with ∆v (see Figure 11 of Wada et al. 2009).
The collision time t coll is the time taken for dust grains of the same size to collide. With the number density of dust grains of n dust , the collision time is given as where is the number density of dust grains, and is the dust scale-height determined by the balance between vertical settling and turbulent diffusion (Dubrulle et al. 1995). In turbulent disks as considered in this study, turbulence is the main source of relative velocity (Johansen et al. 2014). Based on the analytic study by Ormel & Cuzzi (2007), we model the turbulence-induced relative velocity as ∆v = min (∆v 1 , ∆v 2 ) , where Re = αΣ gas σ mol /(2m gas ) is the turbulent Reynolds number. Equation (29) continuously connects the asymptotic expressions for ∆v in the limits of St Re −1/2 and Re −1/2 , i.e., ∆v 1 and ∆v 2 (see Equations (27) and (28)  It should be noted that we employ a single α value to describe both the angular momentum transfer efficiency (the so-called viscosity) and the velocity dispersion of gas turbulence. While the former drives gas accretion (Equation (3)), the latter induces dust diffusion (Equation (28)) and collision (Equation (29)). The two physically distinct quantities would be comparable in magnitude if nearly isotropic turbulence is the main driver of the accretion. However, they can have significantly different values in self-gravitating disks, where the accretion stress from the self-gravitating spirals dominates over the Reynolds stress from turbulence (Baehr & Zhu 2021). Observations of protoplanetary disks that exhibit a high level of dust settling also point to dust diffusivity much lower than the angular momentum transport efficiency (Pinte et al. 2016;Ribas et al. 2020;Villenave et al. 2022). In Section 5.3, we assess the effects of α acc = α turb on our results.

Boundary Conditions and Parameter Choices
We numerical solve Equation (24) from the disk outer edge at r = r d inward. The maximum grain size at r = r d is taken to be a max = a min = 0.1 µm.
We consider a wide range of physical parameters (Table 1). In the fiducial model, we consider the central massive protostar of M = 20M and L = 4×10 4 L with the disk of r dist = 150 au. For the gas accretion rate, we explore over two orders of magnitudeṀ gas = 10 −5 -10 −3 M yr −1 , selecting 10 −4 M yr −1 as the fiducial value. As explained in Section 1, there is a large degree of uncertainty in the fragmentation velocity of silicate grains. To investigate how the v frag value affects dust coagulation, we vary v frag between 1-100 m s −1 , with 10 m s −1 being the fiducial value.
Our ultimate goal is to constrain the sticking property of silicate grains by comparing our model and observations. To date, the massive protostellar disk of GGD27-MM1 is the only object where the maximum silicate grain size is observationally constrained (Girart et al. 2018). Thus, we also perform model calculations with the parameters of GGD27-MM1 (Añez-López et al. 2020), searching for the value of v frag that yields the observationally inferred grain size, analytically in Section 5.1 and numerically in Appendix A.

ANALYTIC ESTIMATES
As already demonstrated by previous studies (e.g., Birnstiel et al. 2009Birnstiel et al. , 2012Okuzumi et al. 2016), one can easily estimate the maximum size of dust grains when their growth is limited by either collisional fragmentation or radial drift. Our numerical calculations presented in Section 4 show that fragmentation is the dominant growth barrier (see Section 4) in massive protostellar disks. Before presenting the numerical results, we here derive analytic expressions for the fragmentationlimited grain size in self-gravitating disks, making use of the fact that the gas density and temperature in the disks can be expressed in terms of the Toomre parameter Q. Such expressions are useful for understanding how the maximum grain size in massive protostellar disks depends on various parameters.
We begin by using Equations (2) and (8) to rewrite the gas surface density and sound velocity as An advantage of these expressions is that they do not explicitly involve T and therefore do not explicitly depend on the disk's heating mechanisms. These expressions are useful for self-gravitating disks, where the accretion driven by the GI tends to regulate the value of Q and α(Q) to ≈ 1-1.4 and ∼ 0.1, respectively (Zhu et al. 2010, see also Section 2.1). When collisional fragmentation limits dust growth, the grains grow until the collision velocity ∆v reaches the fragmentation velocity v frag (Birnstiel et al. 2011(Birnstiel et al. , 2012. Therefore, the fragmentation-limited grain radius a frag is determined by the relation ∆v(a frag ) = v frag .
The corresponding Stokes numbers are Because St determines v r,dust , one can also derive estimate Σ dust =Ṁ dust /(2πr|v r,dust |). As shown below, the dust-to-gas surface density ratio Σ dust /Σ gas in massive protostellar disks tends to become radially uniform because v r,dust ≈ v r,gas . From Equations (3) and (19) with St 1, the relative difference in the dust and gas radial velocities can be described as v r,dust − v r,gas v r,gas = 2 3 dlnP dlnr For typical parameters for massive protostellar disks, Equation (34b) and (34c) show that St frag α floor < α, indicating that the radial velocity difference between the gas and dust is small. Dust vertical settling is also negligible (Equation (28)).

NUMERICAL RESULTS
Here we present the results of the numerical calculations for the models listed in Table 1. For all models presented here, the entire part of the disk has temperatures above 200 K and thus can be regarded as being inside the water snowline. Figure 2 shows the steady-state disk structure for the fiducial disk model. Figure 2(a)-(c) plot the maximum grain size a max , collision velocity ∆v, and the Stokes number St as a function of radial distance r. These plots indicate that the grain growth in this model is limited by collisional fragmentation. At the disk's outer edge lying at 150 au, the dust grains supplied from the envelope grow locally to ∼ 100 µm in size (Figure 2(a)). This local growth terminates as the collisional velocity ∆v reaches the fragmentation velocity v frag (Figure 2(b)). After that, the grains accrete inward, keeping the size determined by the balance between coagulation and fragmentation (∆v ≈ v frag ). The Stokes number of the inward accreting grains is ∼ 10 −4 , nearly independent of r (Figure 2(c)). This value of St is much smaller than α floor = 10 −2 , indicating that the grains' radial velocity is dominated by the gas accretion velocity, as mentioned in Section 3.

The Fiducial Model
The disk structure is regulated by the GI torque. Figure 2(d)-(f) present the gas and dust surface den-sities, Toomre parameter Q, and viscous parameter α in the fiducial model. Both the gas and dust surface densities increase toward smaller r, roughly following a single power of r −1.5 (Figure 2(d)). The radial dependence of Σ gas agrees with our analytic estimate, Equation (30), showing Σ gas ∝ Ω K ∝ r −3/2 . The entire disk is marginally unstable at Q ∼ 1, inducing a GI torque of α ∼ 0.02-0.2 (Figure 2(e) and (f)). Because α St, the gas and dust accrete toward the central star at nearly the same velocity (see Section 3), resulting in a radially constant dust-to-gas surface density ratio of Σ dust /Σ gas ≈ f dg,ISM = 0.01 (Figure 2(d)).
Because Q is radially nearly constant, the analytic estimates derived in Section 3 can be used to predict a max and St for the accreting grains. This is demonstrated in Figure 2(a) and (c), where the dashed lines show Equations (33) and (34) for Q = 1.2. We find that the analytic estimates well reproduce the numerical results at r 40 au. The agreement is less good at r 40 au because the actual value of Q depends weakly on r. Since we are primarily interested in the outer ∼ 100 au region, which is the main target of radio imaging observations, we use Q = 1.2 as the reference value for self-gravitating disks. Figure 2(g)-(i) show the temperature T , optical depth τ , and opacity κ abs for the fiducial model. Around r ∼ 100au, the disk heating is dominated by stellar irradiation, yielding T ≈ T irr (Figure 2(g)). The temperature appreciably exceeds T irr at r 100 au, where the accretion heating dominates owing to the deeper gravitational potential and higher optical depth (Figure 2(h)). Because the maximum grain size a max increases as the grains accrete toward the central star, the opacity decreases toward smaller r (Figure 2(i)).
To see how dust coagulation contributes to the disk structure shown above, we also show in Figure 2 the calculation results for a uniform, constant grain radius of 0.1 µm (see the black lines). One can see that dust coagulation leads to smaller κ abs and τ , resulting in lower T (see Figure 2(i), (h), (g), respectively). In contrast, the surface densities Σ gas and Σ dust are almost unchanged (Figure 2(d)). This is consistent with Equation (30), which shows that in a self-gravitating disk whose accretion is regulated by the GI torque (Q ∼ 1 and α ∼ 0.1), the surface density is determined independently of the disk temperature.

Dependence on the Fragmentation Velocity
We here investigate how dust coagulation and the disk structure depend on the fragmentation velocity v frag . Figure 3 shows the results for v frag = 1, 10 (fiducial), and 100 m s −1 .
All three cases show ∆v ≈ v frag (Figure 3(b)), suggesting that fragmentation mainly limits dust coagulation. The radial drift barrier only gives a minor contribution to the coagulation limit, although its effect is appreciable in the case of v frag = 100 m s −1 , which shows that ∆v is slightly lower than v frag . Figure 3(a) and (c) show that a higher v frag leads to a larger maximum grain size and a larger Stokes number, consistent with the analytic estimates given by Equations (33) and (34).
The gas surface density is insensitive to the fragmentation velocity (Figure 3(d)). This is because the disk accretion is self-regulated to the marginally unstable state with Q ∼ 1 and α ∼ 0.1 (Figure 3(e)-(f)), which only implicitly depends on T (see also Equation (30)). The relation Σ rmdust ≈ 0.01Σ gas seen in the fiducial model approximately holds as long as v frag < 100 m s −1 . For v frag = 100 m s −1 , Σ dust is slightly shallower than Σ gas due to the grains' non-negligible radial inward drift relative to the gas.
The disk temperature is lower with higher v frag (Figure 3(g)), since the optical depth and opacity are lower with larger grain size (Figure 3(h)-(i)). In the most optically-thin case of v frag = 100 m s −1 , stellar irradiation is the dominant heat source at r 30 au. As mentioned above, any change in temperature has little impact on the gas surface density, which is the characteristic of self-gravitating disks.

Dependence on the Accretion Rate
Next, we present the model calculations with various accretion rates, fixing the fragmentation velocity to v frag = 10 m s −1 . Figure 4 shows the results foṙ M gas = 10 −5 , 10 −4 (fiducial), and 10 −3 M yr −1 .
In all models, dust coagulation is regulated by the collisional fragmentation (Figure 4(b)). The maximum grain size is smaller for higher accretion rates (Figure 4(a)). The main reason for this trend is that a higheṙ M gas leads to higher α and T (see Figure 4(f)-(g)), both of which increase the turbulence-induced collision velocity ∝ α 1/2 c s ∝ (αT ) 1/2 . However, theṀ gas dependence of the maximum grain size is minor compared with the v frag dependence shown in Figure 3.
The disk with higherṀ gas has higher surface-densities Σ gas and Σ dust (Figure 4(g)). However, the dependence is rather weak, with a factor of 100 increase inṀ gas only resulting in a factor of several increases in Σ gas and Σ dust . Because a max increases withṀ gas , a higheṙ M gas yields a lower κ abs (Figure 4(i)). The increase in κ gas and Σ gas results in larger τ and hence higher T for higherṀ gas (Figure 4(g)-(h)).
In the case ofṀ gas = 10 −5 M yr −1 , the regions inside and outside r ≈ 30 au exhibit distinct radial structures. This is because, in this low-Ṁ gas disk, the inner part is gravitationally stable with Q > 1.4 (Figure 4(e)) and α ≈ α floor = 0.01 (Figure 4(f)).
In the lower and higherṀ models, the numerical results deviate more from the analytic estimates than in the fiducial model, because the Q values are greater or less than the reference value of 1.2 in those cases (Figure 4(e)). However, the deviations from the analytic estimates are only a factor of a few at the outer region of ∼ 100 au.  In the previous section, we have shown that collisional fragmentation limits the growth of silicate grains in massive protostellar disks. The maximum grain size of the silicate grains is primarily determined by their threshold fragmentation velocity v frag . The maximum size is only weakly dependent on the accretion rate and, more importantly, independent of the disk temperature owing to the self-regulating nature of the self-gravitating disks (see also Section 3). These suggest that measurements of the maximum grain size in massive protostellar disks can be used to constrain v frag for silicates, which is currently highly uncertain, as introduced in Section 1. In this subsection, we derive analytic formulas that are useful for this purpose. An application of the formulas is presented in Section 5.2.
Equations (33a)-(33c) in Section 3 provide analytic estimates for the fragmentation-limited grain size a max for given v frag . Now we invert them into formulas that return v frag for given a max . The result reads v frag = min (v frag,1 , v frag,2 ) , v frag,1 ≈ 35α 2/3 Q 5/6 a max 100 µm Putting Q = 1.2 and α(Q = 1.2) = 0.14 for typical selfgravitating disks and ρ int = 3.0 g cm −3 for silicates, we have v frag = min (v frag,1 , v frag,2 ) , (37c) Importantly, these estimates depend only weakly on the gas accretion rate and stellar mass, which are difficult to determine precisely from observations.

The Case of the GGD27-MM1 Disk and Its Implication for Rocky Planetesimal Formation
To date, the disk of GGD27-MM1 is the only massive protostellar disk for which a measurement of the maximum grain size from observations is available (Girart et al. 2018). We here apply Equation (37) to the observations to infer the fragmentation velocity for the silicates in this disk.
The 1.14 mm polarimetric observations of the GGD27-MM1 disk (Figure 3 of Girart et al. 2018) showed that the disk emission around r ≈ 100 au is polarized, with the polarization vector aligned with the disk's minor axis. Such polarization features can be explained by the self-scattering of the thermal emission by dust grains with a maximum radius of ∼ 100 µm (Kataoka et al. 2016;Yang et al. 2017). In principle, the grain size constraint derived from self-scattering models may not reflect the size of the grains lying at the disk midplane if Figure 5. Fragmentation-limited grain radii a frag (Equation (32)) for the GGD27-MM1 disk (see text for the adopted stellar parameters) with different values of the fragmentation threshold v frag . The vertical bar indicates the inferred range of the maximum grain radius amax at r = 100 au from the assumption that the polarized emission from this region is due to dust self-scattering. The inferred amax is consistent with v frag = 6.4-11 m s −1 (blue shaded area) and rules out v frag = 1 and 50 m s −1 (gray lines).
the disk is optically thick and if larger grains have settled to the midplane (Ueda et al. 2021). However, as discussed in Section 3, dust settling in massive protostellar disks should be negligible (St α), and thus the grains responsible for the polarized emission should also represent the grains at the midplane.
To give a more quantitative estimate of the maximum grain size a max , we focus on the polarized emission at r = 100 au in the GGD27-MM1 disk. According to Figure 3 of Girart et al. (2018), the emission from this position has a polarization degree of ∼ 1.5% (we here assume the source distance of 1.4 kpc; Añez-López et al. 2020). We then follow Girart et al. (2018) and convert the polarization degree into a max based on the self-scattering model of Kataoka et al. (2016, their Figure 3). We estimate a max ≈ 100-300 µm for the polarization degree of ∼ 1.5%. Our estimate here is still crude because the selfscattering model we use here was not originally designed to predict the polarization degree of the GGD27-MM1 disk. Moreover, the model of Kataoka et al. (2016) assumes perfectly spherical grains, but more recent studies show that the polarization degree also depends on the shape of the grains (Kirchschlager & Bertrang 2020;Tazaki & Dominik 2022;Lin et al. 2023). Future dedicated radiative transfer modeling of this particular system will allow a more accurate estimate for a max .
We now use Equation (37) to infer the fragmentation velocity of silicates in the GGD27-MM1 disk. We assume M = 20 M andṀ gas = 7 × 10 −5 M yr −1 for GGD27-MM1 (Añez-López et al. 2020). Equation (37) then yields v frag ≈ 6.4-11 m s −1 for the estimated maximum grain size of a amax ≈ 100-300 µm at r = 100 au. This is also shown in Figure 5, which plots the fragmentation-limited grain size a frag (Equation (33)) for various v frag values. Importantly, our estimate rules out the commonly quoted value of v frag = 1 m s −1 for silicates, along with the most sticky scenario of v frag = 50 m s −1 in the literature (see Section 1). It is worth noting that our work is consistent with the conclusion drawn by Liu et al. (2021), who similarly estimated a fragmentation velocity of ∼ 10 m s −1 for silicate dust based on observations of the accretion burst FU Ori Disk. For completeness, we show in Appendix A the radial structure of the GGD27-MM1 disk obtained from direct numerical calculations of the model equations, confirming that v frag ≈ 10 m s −1 best explains the polarized emission at r = 100 au.
The above exercise is only one example based on the observations of a single system at a single wavelength. Clearly, more case studies with other massive protostellar disks and at various wavelengths are needed to give a more robust constraint on v frag for silicates. With this caveat in mind, it is still interesting to discuss what the derived fragmentation velocity of v frag ≈ 10 m s −1 would imply for the formation of rocky planetesimals around low-mass protostars. As already mentioned in Section 1, the maximum value of grain collisional velocity in protoplanetary disks is ≈ 20-50 m s −1 , even if no strong turbulence is present (Johansen et al. 2014). With v frag ≈ 10 m s −1 , aggregates of silicate grains in protoplanetary disks are likely to experience catastrophic fragmentation at some point, meaning that rocky planetesimals would not form solely via the coagulation of silicates. The streaming and gravitational instabilities would be necessary to account for rocky planetesimal formation. Interestingly, the fragmentation velocity we derived is still considerably higher than the commonly assumed value of 1 m s −1 (Güttler et al. 2010). Because the streaming and gravitational instabilities favor large grains whose aerodynamical coupling to the surrounding gas is moderate, our estimate for v frag suggests that rocky planetesimal formation via these mechanisms may be more efficient than previously anticipated.

Caveats of Our Model
This work is only the first step toward a full understanding of dust evolution in massive protostellar disks and rests on several assumptions and simplifications. Here, we describe some caveats of our model.
Our model uses Epstein's drag law to compute the dust stopping time and thus the Stokes number. We have confirmed that the condition for Epstein's law, a max λ mfp , holds in all disk models presented in this study, except in the inner region of r < 20 au with the high fragmentation velocity of v frag = 100 m s −1 . Therefore, this assumption does not affect our main conclusions. In any case, it is straightforward to extend our model beyond the Epstein drag regime.
We have focused on the steady part of the disk structure by enforcing steady accretion. However, selfgravitating disks are often highly time-varying due to accretion bursts (e.g., Meyer et al. 2017). Vorobyov et al. (2022) simulated dust collisional evolution during accretion bursts in low-mass protostellar disks. In their simulations, the snowline moves away from a few au to several dozen of au during bursts due to increased stellar luminosity. The snowline excursion would induce a drop in the maximum grain size in the temporal silicatedust region if bare silicate grains are more fragile than ice grains. However, in massive protostellar disks, the snowline already lies at as far as several hundred au from the central star even without bursts. In these disks, accretion bursts would not affect the grain size in the r 100 au region.
We have neglected fragmentation of self-gravitating disks, which in fact could occur if Q < 0.6 (e.g., Takahashi et al. 2016). However, in all disk models presented in Section 4, the Toomre parameter is larger than 0.6 in the entire disks, indicating that disk fragmentation would be negligible.
We adopted the α value as the function of Toomre parameter Q (Equations (6) and (7)), and selected Q = 1.2 and α(Q = 1.2) = 0.14 as the reference values. Although this function is known to reproduce the qualitative of self-gravitating protostellar disks (e.g., Takahashi et al. 2013), the accuracy of the function has not been evaluated so far. For instance, if the GGD27-MM1 disk is moderately unstable (Q ≈ 1) but highly turbulent (α ≈ 1), Equation (36) predicts v frag ≈ 20 m s −1 instead of v frag ≈ 10 m s −1 . Because this fragmentation velocity is comparable to the low end of the range of the maximum grain collision velocity (20-50 m s −1 ), planetesimal formation by dust coagulation alone may be barely possible with this higher value of v frag .
As for the α parameter, we also note that it is non-trivial whether turbulence and angular-momentum transfer can be represented by a single value of α, especially in self-gravitating disks. For completeness, we re-derive the formula of the fragmentation velocity, distinguishing two α values: α turb represents turbulence intensity and α acc causes accretion through angular momentum transport, v frag = min v frag,1 α turb α acc 3/4 , v frag,2 α turb α acc 1/2 .
(38) This is equivalent to Equation (37) if α turb = α acc . In self-gravitating disks, non-turbulent (coherent) spiral arms could give a considerable contribution to the gravitational accretion stress. In this case, one would have α turb < α acc , and therefore v frag estimated by Equation (38) would be smaller than our reference value of ≈ 10 m s −1 , which would strengthen our conclusion that rocky planetesimal formation via silicate coagulation alone is unlikely (Section 5.2).
We evaluate the turbulent Reynolds number Re, assuming that molecular viscosity determines the smallest eddy scale. However, in partially ionized media, the turbulence cutoff scale may be much larger due to magnetohydrodynamic waves and radiative cooling (Xu et al. 2016;Silsbee et al. 2021). Despite the difficulty in determining this scale, we examine the impact of Re uncertainty. From Equations (29b) and (36b), we find v frag,1 ∝ Re 1/4 . For example, if Re is four orders of magnitude smaller than the value estimated with molecular viscosity, the fragmentation velocity inferred from the GGD27-MM1 observations would decrease to ≈ 1 m s −1 , further supporting our conclusion that v frag for rocky grains are not as high as 50 m s −1 .
The dust collisional velocity of Equation (29) is based on the Kolmogorov turbulence model (Ormel & Cuzzi 2007). However, earlier studies suggested that turbulence energy spectra in astronomical disks may deviate from the Kolmogorov law (e.g., Iroshnikov 1964;Kraichnan 1965). Gong et al. (2021) expanded upon this by developing collisional velocity formulas for arbitrary turbulence models represented by power-law spectra. They found that collisional velocities for the Iroshnikov-Kraichnan turbulence and the turbulence caused by magnetorotational instabilities are higher than those for the standard Kolmogorov turbulence. If the GGD27-MM1 disk has such a turbulence energy spectrum, we estimate that the fragmentation velocity might be as high as ∼ 50-150 m s −1 .
Finally, we discuss the validity and limitations of the single-size approach and the interpretation of radio polarimetric observations. In general, the single-size approach well approximates the evolution of the maxi-mum grain size and dust surface density when as long as the grain size distribution is top-heavy, i.e., the largest grains dominate the total dust surface density (Birnstiel et al. 2012;Sato et al. 2016). However, the single size approach itself does not derive the size distribution, so we have assumed a power-law size distribution to compute the disk's vertical optical depth (Section 2.2). We expect that this assumption would not affect the results presented in this study, because the temperature of selfgravitating disks is self-regulated such that Q ∼ 1 is maintained. The adopted size distribution is top-heavy (see footnote 1), so our model is internally consistent in this respect. However, depending on the cascade processes of collisional fragmentation (e.g., Kobayashi & Tanaka 2010), the grain size distribution may become bimodal rather than single-peaked (H. Kobayashi, priv. comm.). Our single-size approach would not be adequate for treating such distribution. The uncertainty in the grain size distribution, as well as in the grain shape, could affect the grain size constraint derived from millimeter polarimetric observations based on the dust self-scattering scenario (see also Section 5.2). Further investigation with detailed treatments of the collisional fragmentation and radio scattering processes will be required in future work.

SUMMARY
The threshold fragmentation velocity of silicate grains is a critical parameter to understanding the formation of rocky planetesimals, but it is currently highly uncertain. Silicate dust regions of 200 K are too small to observe in protoplanetary disks around low-mass protostars, while they are observable in disks around hotter massive protostars. In this study, we for the first time developed a theoretical model of dust evolution in massive protostellar disks, and constrained the fragmentation velocity of silicate dust.
We solved the gas disk structure and the dust coagulation self-consistently. For the gas disk structure, we took into account both stellar irradiation and accretion heating, together with efficient angular momentum transfer of self-gravitating disks. For the dust evolution, we con-sidered coagulation and fragmentation, as well as radial drift and vertical settling. We found that the maximum grain size is limited by the collisional fragmentation rather than by the radial drift in massive protostellar disks. We also derived the analytic formula of the fragmentation-limited grain size in self-gravitating disks (Equation (33)), which reproduces the numerical results well. The maximum grain size is larger for the higher fragmentation velocity, and only weakly depends on the heating mechanism and the accretion rate. This suggests that measurements of the maximum grain size in massive protostellar disks can be used to constrain the fragmentation velocity.
Based on the results of our analytic and numerical calculations, we derived a new simple formula to constrain the fragmentation velocity, i.e., Equation (37). Using this formula with the maximum size of silicate grains in the massive protostellar disk of GGD27-MM1 (Girart et al. 2018), we estimated the threshold fragmentation velocity of silicate grains as v frag ≈ 10 m s −1 . This obtained fragmentation velocity is lower than the maximum collisional velocity expected in protoplanetary disks around low-mass protostars (e.g., Johansen et al. 2014). This implies that rocky planetesimals form not by dust collisional growth but by other mechanisms, such as the gravitational instability of the dust layer.
We note that measurements of the maximum size of silicate grains in massive protostellar disks are still limited and indefinite. Further polarization observations at multiple wavelengths toward multiple objects are required to accurately constrain the maximum size of silicate dust, and thus the fragmentation velocity.  (Table 1) and v frag = 1,, 10, and 100 m s −1 (as indicated in panel a). The thin dashed lines represent the analytic solutions, assuming Q = 1.2. For reference, the range of the maximum grain radius suggested by the 1.14-mm polarization observation (Girart et al. 2018) is shown by the black vertical line in panel (a). The model with v frag = 10 m s −1 reproduces the observed grain size, as predicted by the analytic method (see Section 5.2). protostellar disk for which silicate grain coagulation has been confirmed (Girart et al. 2018).
In Section 5.2, we used our simple analytic formula of Equation (37) and found that the fragmentation velocity of v frag ≈ 10 m s −1 can reproduce the observed grain size, i.e., a max = 100-300 µm at r = 100 au. To confirm the evaluated v frag value, we here perform the detailed numerical model calculations with the parameters of the GGD27-MM1 disk (Table 1). Figure 6 shows the results of the GGD27-MM1 model calculations with v frag = 1, 10, and 100 m s −1 . The basic behaviors of all physical quantities are the same as in the fiducial and comparison models in Section 4.2 (see Figure 3). The maximum grain size is larger for higher fragmentation velocity. It can be seen that the model with v frag = 10 m s −1 falls within the range of the observed grain size (panel a), which confirms our analytic estimation.