Monte Carlo Simulation of Dust Particles in a Protoplanetary Disk: Crystalline to Amorphous Silicate Ratio in Comets

Observationally inferred crystalline abundance in silicates in comets, which should have been formed in the outer region of a protoplanetary disk, is relatively high (~ 10-60%), although crystalline silicates would be formed by annealing of amorphous precursors in the disk inner region. In order to quantitatively address this puzzle, we have performed Monte Carlo simulation of advection/diffusion of silicate particles in a turbulent disk, in the setting based on pebble accretion model: pebbles consisting of many small amorphous silicates embedded in icy mantle are formed in the disk outer region, silicate particles are released at the snow line, crystalline silicate particles are produced at the annealing line, the silicate particles diffused beyond the snow line, and they eventually stick to drifting pebbles to come back to the snow line. In a simple case without the sticking and with a steady pebble flux, we show through the simulations and analytical arguments that crystalline components in silicate materials beyond the snow line is robustly and uniformly ~ 5%. On the other hand, in a more realistic case with the sticking and with a decaying pebble flux, the crystalline abundance is raised up to ~ 20-25%, depending on the ratio of decay and diffusion timescales. This abundance is consistent with the observations. In this investigation, we assume a simple steady accretion disk. The simulations coupled with the disk evolution is needed for more detailed comparison with observed data.


INTRODUCTION
Crystalline silicates, such as forsterite, are formed by re-condensation after sublimation at T 1400 K of amorphous precursors or annealing of them at T 1000 K (e.g., Wooden 2008). However, they are observed in outer cold regions of protoplanetary disks and the Stardust samples (e.g., Brownlee et al. 2006;Ciesla 2011, and references therein). Infrared observations of comets also show 10-60% of silicates are crystalline (e.g., Sitko et al. 2011;Harker et al. 2011, also see compiled data for several comets in Shinnaka et al. 2018). Comets were formed in outer region of a protoplanetary disk, such as Kuiper belt or Uranus-Neptune zone, and scattered to Oort Clouds (e.g., Dones et al. 2015). Some mechanism should have transported the crystalline silicate particles that experienced T 1000 K in the inner disk region to the outer region. The radial turbulent diffusion is one of available transporting mechanisms (e.g., Gail 2001;Bockelée-Morvan et al. 2002;Cuzzi et al. 2003), while centrifugal jets (Shu et al. 1996) or global disk flow (e.g., Keller & Gail 2004;Ciesla 2007;Desch 2007) could also be available. In this paper, we focus on the radial turbulent diffusion. Ciesla (2010Ciesla ( , 2011) developed a 3D Monte Carlo simulation code for turbulent diffusion of crystalline silicates in a protoplanetary disk to successfully show that the particles initially at ∼ 5 au can diffuse out beyond 20 au. However, because the initial distributions of amorphous and crystalline silicates were not clear, the results were not compared quantitatively with the observed abundance of crystalline components in silicates in the comets. Pavlyuchenkov & Dullemond (2007) numerically solved one-dimensional diffusion equation for crystalline materials formed at the formation front in the inner disk region. Arakawa et al. (2021) discussed the modulations from the Pavlyuchenkov & Dullemond (2007)'s results by silicate particle growth and by drift due to gas drag in the disks with and without a pressure bump created by magnetically driven disk winds. Both papers assumed that amorphous silicates are stationarily distributed with the uniform solid-to-gas ratio in the steady accretion disk, and Arakawa et al. arXiv:2201.05507v2 [astro-ph.EP] 21 Jan 2022 (2021) concluded that the crystalline abundance is very low unless a global pressure bump is introduced. However, their setting on silicate particles is not consistent with the new model of planetary accretion called "pebble accretion" (e.g., Ormel & Klahr 2010;Lambrechts & Johansen 2014).
In the pebble accretion model, pebbles with cm-m size are formed in disk outer regions, drift all through the disk, and some fraction of the drifting pebbles are accreted by planetary seeds (e.g., Lambrechts & Johansen 2014;Sato et al. 2016). The pebbles would consist of many small amorphous silicate particles embedded in icy mantle and the silicate particles are released by ice sublimation at the snow line (e.g., Saito & Sirono 2011;Morbidelli et al. 2015;. The reason to consider this pebble structure is that the threshold velocity for fragmentation of silicates would be generally one order smaller than that of ices (e.g., Blum & Wurm 2000;Zsom et al. 2011;Wada et al. 2011), although this conventional view is recently challenged (Kimura et al. 2015;Gundlach et al. 2018;Musiolik & Wurm 2019;Steinpilz et al. 2019). Assuming this silicate-ice structure of pebbles, Ida et al. (2021) derived the conditions for the gravitational instability of the pile-up of the released silicate particles by local 2D (r − z) Monte Carlo simulation, following the method by Ciesla (2011).
The motivation of this paper is to apply the above setting of injection of amorphous silicate particles at the snow line for predicting the crystalline silicate ratio in the total (silicate and amorphous) silicates in disk outer regions. We find that about 5% of the released amorphous particles drifting inward from the snow line come back to disk outer regions by outward diffusion, after experiencing high temperature T 1000 K to be transferred to crystalline particles. To calculate the global evolution of silicate dust particle distributions, we adopt global 3D (x − y − z) Monte Carlo simulation by Ciesla (2011). We do not apparently incorporate a coupled silicate dust growth such as the model by Misener et al. (2019), because the results here do not change, as long as Stokes number of the silicate particles is smaller than the α parameter of the turbulent mixing, while the effect of the silicate growth is discussed in Appendix. On the other hand, we include the effect of rapid inward drift of the silicate particles that coagulate with drifting pebbles, which was discussed by Misener et al. (2019). Because we can calculate surface densities both of crystalline and amorphous silicate particles as a function of orbital radius, a quantitative comparison between theoretical prediction and the observation of comets can be done.
In Section 2, we describe the protoplanetary disk model that we adopt and Monte Carlo methods for advection due to gas drag and turbulent diffusion of silicate particles in the disk. The sticking rate of the particles to drifting pebbles beyond the snow line is also evaluated. In Section 3, we describe the analytical steady solution of the radial advection-diffusion equation to explain the simulation results with and without silicates sticking to pebbles. In Section 4, we show the simulation results for the simple case neglecting the sticking (Section 4.1), the case with sticking to icy pebbles (Section 4.2), and the case with the attenuation of a pebble flux due to the depletion of solid materials in the reservoir in the outer disk (Section 4.3). While the sticking to the pebbles inhibits outward diffusion of the silicate particles, the decaying pebble flux tends to diminish the sticking effect. We show that the radial diffusion in the case with the decaying pebble flux can quantitatively explain the relatively high abundance in the comets. Section 5 is the conclusion and discussion.

Disk model
We adopt the quasi-steady composite disk model derived by  and Ida et al. (2019). For simplicity, the depletion due to photoevaporation is neglected and we assume that the disk evolution is regulated by turbulent viscous diffusion. In general, the viscous heating is dominant in the inner disk regions, while the irradiative heating dominates in the outer disk regions. We composite these two regions by the disk mid-plane temperature T = max(T vis , T irr ), where T vis and T irr are the temperature in the viscous-heating and irradiation dominated regions given respectively by and whereṀ g is the gas accretion rate through the disk and α is a parameter of turbulent viscosity, which are scaled by the fiducial parameters in our disk model,Ṁ g = 10 −7 M /yr and α = 10 −2 .
The temperature in the viscous-heating depends on the opacity. Equation (1) corresponds to relatively low opacity with relatively large silicate particles with the size ∼ 0.1 mm and the temperature is higher for µm size particles (Oka et al. 2011). As we will show in section 4.1, as long as both the snow and the annealing lines are in the viscousheating dominated regime, the predicted crystalline silicate abundance is not affected by the numerical factor of T ; the abundance depends only on the r-dependence of T . From Eqs (1) and (2) au. (3) The gas disk aspect ratios corresponding to Eqs. (1) and (2) are From the assumption of steady accretion,Ṁ g = 3πΣ g αH 2 g Ω, where Σ g is the disk gas surface density and Ω is the Keplerian frequency, the gas surface densities the viscous-heating and irradiation dominated regions are We use the gas disk aspect ratio h g = H g /r = max(h g,vis , h g,irr ) and the gas surface density Σ g = min(Σ g,vis , Σ g,irr ). We simply set the snow line by T = 170 K. From Eqs. (1) and (2), the snowline radius is given by r snow = max(r snow,vis , r snow,irr ), where The snow line is determined by viscous heating whenṀ g 5.0 × 10 −9 M /yr for α = 10 −2 . Similarly, we set the silicate line radius r sil as the location of silicate sublimation temperature (T = 1400 K), and the annealing line r annl by T = 1000 K. When the viscous heating dominates at the snow line, it also dominates at the silicate and annealing lines and they are given by r annl 0.14 r snow (10) If the irradiation dominates, the radial dependence of T is weaker and r annl 0.016 r snow . We use M * = M and L * = L . As a result, the disk parameters are only the disk accretion rateṀ g and the viscosity parameter α. The fiducial values of our disk model areṀ g = 10 −7 M /yr and α = 10 −2 . We will also show the results with α = 10 −3 . Another important parameter is the pebble to gas mass flux ratio, F p/g =Ṁ peb /Ṁ g (wherė M peb is the pebble mass flux). Because pebbles start effective drift when they grow enough that the drift timescale becomes shorter than the growth timescale, the pebble mass flux is regulated by the outward migration of the pebble formation front where the two timescales coincide (Lambrechts & Johansen 2014;.  estimated that F p/g ∼ 0.03(α/10 −2 ) −1 (t/10 6 yr) −8/21 . Because F p/g depends on pebble growth model and disk structure, we treat F p/g as a parameter. How to reflect the value of F p/g to Monte Carlo simulation is described in section 2.3.

Overview
Our models for gas disk and particle evolution in different disk regions are illustrated in Figure 1. Pebbles are formed in outer disk regions and drift inward all through the disk due to gas drag. When icy pebbles are formed, many small "amorphous" silicate dust particles (represented by dark blue dots) are embedded in the icy mantle (represented by light blue), because we assume that collisions between silicate and icy particles are as sticky as ice-ice collisions (the threshold collision velocity for fragmentation, v frag 10 m/s (e.g., Wada et al. 2013;Gundlach & Blum 2015), while silicate-silicate collisions are not sticky, v frag ∼ 1 m/s (e.g., Wada et al. 2013). The Stokes number τ s,peb is given by Eq. (37) ( 0.1). The typical size of pebbles is 10-100 cm (Eq. (36)).
When the pebbles drift inward to r < r snow , the ice mantle is sublimated and the embedded amorphous silicates are released. During the drift, crystalline and amorphous silicate particles that diffused out from inner regions (see below) also stick to the surface of drifting icy pebbles with the probability derived in section 2.4. They are also re-released at the snow line.
The Stokes number of silicate particles is assigned to be τ s,sil , which we use τ sil = 10 −5 , corresponding to the dust particle size of 10µm for α = 10 −2 . The radial drift velocity of the silicate particles is given by Eq. (22). It shows that the radial velocity of the silicate particles and the calculation results are independent of the detailed value of τ s,sil (equivalently, the silicate particle size), as long as τ s,sil < α (diffusion-dominated for the silicate particles) and τ s,peb > α (drift-dominated for pebbles). We adopt a relatively severe rebounding/fragmentation limit to neglect growth of silicate particles through mutual collisions. Even if the detailed silicate dust growth model is adopted, the above condition would be satisfied for α = 10 −3 -10 −2 that we use in this paper. Detailed discussions on τ s,sil is given in the Appendix.
Once the silicates arrive at r < r annl , they are identified as "crystalline" silicates (represented by orange dots). Even if they go back to the regions at r > r annl , they keep the state of crystals, because crystalline silicates are stable. When the silicates arrive at r < r sil , they should become gas molecules. We assign τ sil = 0 for them in this region. When they go back to the regions r > r sil , they re-condense as crystalline silicates. About 5% of the crystalline silicate particles diffuse out beyond the snow line and eventually diffuse inward toward the star. If the icy pebbles are still drifting there, the crystalline silicates attach to the pebbles to rapidly return to the snow line again. Once the steady state is established, the silicate loss rate to the star is balanced with their injection rate at the snow line, although individual silicate particles have different radial (sometimes, rather complicated) travel paths. These silicate particles tracking is summarized in Table 1.

Detailed Monte Carlo simulation methods
In order to follow the orbital evolution of silicate particles described in section 2.2.1, we employ the 3D Monte Carlo simulation method for advection and diffusion of silicate particles developed by Ciesla (2011). Ciesla (2011) showed location Stokes number structure r > rsnow τ s,peb (sticking) amorphous/crystalline τ s,sil (not sticking) r annl < r < rsnow τ s,sil amorphous/crystalline r sil < r < r annl τ s,sil → crystalline r < r sil 0 → gas Table 1. Changes of silicates. If the particles stick to icy pebbles, they drift/diffuse with the pebbles of τ s,peb given by Eq. (37), which is usually 0.1. If not, they drift/diffuse with τ s,sil = 10 −5 . At r sil < r < r annl , they change to crystalline silicate particles, while they change to gas molecules at r < r sil . At r > r annl , they do not change structure.
that the surface density evolution of silicate particles determined by the concentration equation, where Σ sil is the surface density of the silicate particles, D is their diffusivity, and v r is their radial advection velocity, is reproduced by the following 3D Monte Carlo simulations of many particles. The changes in the positions in Cartesian coordinates (x, y, z) of each particle after the timestep δt are given by where δt is given by the inverse of the local Keplerian frequency, δt = Ω −1 , and the first and second terms in the right hand sides represent the advection and diffusion. In the diffusion terms, R x , R y and R z are independent random numbers in a range of [-1,1] with a root mean square of 1/ √ 3. Because τ s,sil 1 and the back-reaction due to silicate particle concentration is neglected (assuming the silicate to gas density ρ sil /ρ g 1) in this paper, we identify the diffusivity of the particles, D, as the gas turbulent viscosity (Hyodo et al. 2019), The advection velocities are given by (Ciesla 2010(Ciesla , 2011, In the steady accretion gas disk (∂Ṁ g /∂r = ∂(3πΣ g ν)/∂r = 0), the second terms in Eqs. (17) and (18) vanish. Assuming the vertically isothermal disk with gas density ρ g ∝ exp(−z 2 /2H 2 g ), the second term in Eq.(19) is −αΩz. The radial drift velocity due to gas drag is given by Schoonenberg & Ormel 2017 where v K is the local Keplerian velocity, u ν is the disk gas (inward) accretion velocity given by u ν = 3ν/2r = (3αh 2 g /2)v K , and η is the degree of deviation of the gas rotation angular velocity Ω from Keplerian one, given by Hereafter we approximate 1/(1 + τ 2 s,sil ) as ∼ 1, because τ s,sil 1. Substituting the expression of u ν and Eq. (21) For simplicity, we take C η = 11/8 for T ∝ r −1/2 . The vertical drift velocity v drag,z is given by (Ciesla 2010) Substituting Eqs. (17) to (23) into Eqs. (13) to (15), we obtain We confirmed that these equations reproduce the evolution of a Gaussian-ring (Ciesla 2011, Fig. 1), as shown in Fig. 2 in our paper. We inject super-particles that represent a swarm of silicate particles and suffer specific drag force for realistic particles near the snow line at every timestep δt inj . The initial radii are randomly distributed in the range of [r snow − 0.5∆r 0 , r snow + 0.5∆r 0 ] and the vertical height is determined by a Gaussian distribution of the root mean square ∆z 0 . Following Ida et al. (2021), we adopt ∆r 0 = 0.1 H g and ∆z 0 = H p where H p is a scale-hight of pebbles given by (Dubrulle et al. 1995;Youdin & Lithwick 2007) Shorter (longer) δt inj corresponds to a larger (smaller) number of super-particles with relatively smaller (larger) individual masses in the simulations. We choose δt inj short enough to statistically calculate the silicate surface densities and large enough for numerical simulations. When silicate particles enter the region at r < 0.01r snow , we regard that the particles are lost to the central star and remove them from the simulations.

Surface density of particles
If we are concerned only with relative abundance between amorphous and crystalline silicates at each r, we do not need the information of the surface density of silicate particles. However, to calculate the collision probability of silicate particles with pebbles at r > r snow , we need to specify the surface density of the pebbles. Accordingly, we can also calculate the silicate particle surface density, Σ sil . The surface density of the silicate particles is calculated by the mass of individual silicate (super-)particles given from the pebble mass flux by the method of Ida et al. (2021), as shown below.
The silicate injection mass rate is given byṀ where f sil is the silicate mass fraction in drifting pebbles (f sil = 0.5 in the normal case) andṀ g is the disk gas accretion rate. The mass of an individual super-particle, m, is given bẏ The surface density of silicate particles at r is given with the number of particles ∆N r in the radial width of ∆r around r by where we used Σ g =Ṁ g /2πru ν =Ṁ g /3παH 2 g Ω.

Sticking probability to icy pebbles
While we neglect silicate particle growth due to silicate-silicate collisions, we include the effect of coagulation of silicate particles with pebbles, assuming that silicate-ice collisions are more sticky. By the sticking of diffused-out silicate particles to drifting pebbles at r > r snow , the silicate components are quickly returned to the disk inner region (Misener et al. 2019). We will show later that this effect is very important for the distributions of silicates beyond the snow line.
Assuming the perfect sticking of silicate particles to icy pebbles, the probability P for a silicate particle to collide with any of pebbles during ∆t, which we call "sticking probability," is calculated by where n peb is the number density of the pebbles, R peb is the physical radius of the pebbles, and v is the relative velocity between a silicate dust particle and a pebble. For α 10 −3 , it is dominated by turbulent mixing and is given by (Sato et al. 2016) The factor, H peb /H g , in Eq. (31) reflects that the scale height of the pebbles, H peb , is smaller than that of silicate dust particles, H d ∼ H g , and the collisions are possible only during the silicate particles are passing the region of | z |< H peb . The number density of pebbles is where Σ peb is the surface density of pebbles. Because the scale height of the pebbles in the above equation, H peb , cancels out with the factor in Equation (31), we obtain where we used Eq. (32).
We assume that the bulk density of icy pebbles is ρ bulk 1.0 g cm −3 . From Eq.(22) with replacing τ s,sil by τ s,peb and α τ s,peb , where we usedṀ g = 3πΣ g αH 2 g Ω. In the Stoke regime 1 , the radial drift dominates over collisional growth for pebbles, and the pebble size is approximated by an r-independent form , The corresponding Stoke number of pebbles is given by ( Substituting Eqs. (6) and (37) Therefore, the sticking probability (Eq. (34)) is given by

ANALYTICAL SOLUTION OF SILICATE SURFACE DENSITY
In order to understand the results of our Monte Carlo simulations presented in section 4, we here derive analytical steady solutions to the advection-diffusion equation of concentration given by Eq.(12). We also derive modifications due to the effect of the sticking to pebbles.

Simple case
Equation (12) has two possible solutions in a steady state (∂Σ sil /∂t = 0), whereṀ The former solution represents the inwardly accreting steady flow, M sil = 2πrΣ sil v r,sil .
1 In outer regions far beyond the snow line, the Epstein drag is applied and the pebble size is smaller. However, as we will show in section 4, silicate particles cannot diffuse out so much to the outer regions when the collisions with pebbles are efficient. They can significantly spread only after pebbles are significantly depleted. In that case, while the Epstein drag must be applied, the collision effect is not important. Thereby we only use the expression in the Stokes regime.
Equation (44) means Σ sil (r) ∝ Σ g (r). Assuming D = ν and v r,sil −(3/2)(ν/r) with τ s,sil α, Eq. (45) witḣ The latter (Eq. 42) is the solution of zero net flux, given by which means that the radial gradient is steeper for Σ sil than for Σ g . This dependence of r is equivalent to Eq.(23) with Sc = 1 (D = ν) in Pavlyuchenkov & Dullemond (2007). Setting the radial gradient of gas and dust surface densities Σ g ∝ r −p and Σ sil ∝ r −q , the former and latter solutions are In the viscous-heating and irradiation dominated regions, p = 3/5 (Eq. 6) and p = 15/14 (Eq. 7), respectively. We will show in section 4.1 that the silicate distributions inside and outside the snow line correspond to q = p and q = p + 3/2, respectively.

Case with sticking to icy pebbles
Once silicate particles stick to icy pebbles, they inwardly drift with relatively high velocity as a part of the pebbles. If the sticking is considered, Eq. (12) should additionally have an extra term at r > r snow . Because it is no more exactly solved, we approximately estimate the modification to Σ sil due to the sticking.
The sticking to icy pebbles is a sink of silicates beyond snow line. The sticking rate is given by ζ P Ω in Eq. (39). Thus, the decay rate of silicate outward diffusion flux due to the sticking is where ∆t = t − t 0 and t 0 is the time at which the diffusion flux departs from the snow line. Typical diffusion length from r snow at t is ∆r = < R 2 r > 6αH 2 g Ω∆t = 2αH 2 g Ω∆t.
The solution to Eq. (50) corresponds to whereṀ sil,out,0 is the silicate outward flux at the snow line. Solving this equation in terms of ∆r/r snow , ∆r r where we substituted ζ P given by Eq. (40), and h g is given by Eq.(4) since viscous heating dominates at the snow line wheṅ M g 1.5 × 10 −9 M /yr for α = 10 −3 . The effect of the decay inṀ sil,out can be expressed as modulation to q = p + 3/2 for the zero net flux solution in the simple case in section 3.1. The additional power-law index would be q = − ∂ ln(Ṁ sil,out /Ṁ sil,out,0 ) ∂ ln r − ln(Ṁ sil,out /Ṁ sil,out,0 ) ln((r snow + ∆r)/r snow ) = γ ln(1 + 0.23γ 1/2 ζ r ) . The modified power-law index is given by where we adopt γ 1. ForṀ g = 10 −7 M /y, α = 10 −2 , and F p/g = 0.1, we obtain ζ r = 1.1 and q = 4.4 at 3 au. In the viscous heating dominated regime, q = 6.5. For α = 10 −3 , we obtain ζ r = 0.83, q = 5.7, and q = 7.8 at 3 au. In the irradiation dominated regime (p = 15/14), q is higher by 0.47. Because ζ r is insensitive to the disk parameters (Eq. 54), the modified values of q are not sensitively affected by the disk parameters. In section 4.2, we will show that Eq. (56) reproduces the results of the Monte Carlo simulations.

Simple case
We first show the results in a simple case without the sticking to pebbles and with time-constant F p/g . Figure 3 show snapshots of amorphous (the blue dots) and crystalline (the orange dots) silicate particles and their surface density evolution in a run with fiducial parameters, α = 10 −2 ,Ṁ g = 10 −7 M yr −1 , and F p/g = 0.1. We always adopt M * = M , L * = L , and f sil = 0.5, in this paper. With these parameters, the disk is viscous-heating dominated at r < 7.1 au and irradiation dominated at r > 7.1 au. The disk gas surface density (subsection 2.1) is Σ g = 1.2 × 10 3 (r/1 au) −3/5 gcm −2 (r < 7.1 au) 3.0 × 10 3 (r/1 au) −15/14 gcm −2 (r > 7.1 au).  The snow and annealing lines are at r = 2.0 and 0.28 au. We use the injection time interval δt inj = 1.0 year in all the runs in this paper (in the runs with decaying pebble flux, δt inj is increased after the decay starts). The corresponding super-particle mass is m = 1.0 × 10 25 g. Figure 3a-d are the snapshots at t = 1.8 × 10 4 , 3.6 × 10 4 , 1.8 × 10 5 , and 9.0 × 10 5 years. The characteristic radial diffusion timescale is t diff ∼ r 2 /(6νR 2 ) = r 2 /2ν ∼ (2αh 2 g ) −1 Ω −1 (Eqs. 13 and 14). In the unit of t diff,snow , the diffusion timescale at r = r snow , the time in panels a-d correspond to t ∼ 0.5 t diff,snow , 1.0 t diff,snow , 5 t diff,snow , and 25 t diff,snow , respectively. Silicate particles released at the snow line diffuse both inward and outward with a net inward mean flow inside the snow line. Accordingly, the surface density distribution of silicate particles, which is calculated from the snapshot with Eq. (6), increases with time. We find that the surface density distribution at t = 25 t diff,snow (panel d) does not change any more even if simulations are continued, implying that the distribution at t = 25 t diff,snow is already in equilibrium in the range of r in the plot (r < 40 au).
At t ∼ 0.5 t diff,snow (panel a), about a half of silicate particles released at t ∼ 0 arrive at the annealing line and the crystalline silicate surface density increases to a half of that in the equilibrium state at r ∼ r annl . This timescale is explained as follows. The radial drift timescale for silicate particles by disk gas accretion flow is ∼ r/v r ∼ 2r 2 /3ν ∼ t diff . In the viscous heating dominant region, h g ∝ r −1/20 and accordingly t diff ∝ h −2 g Ω −1 ∝ r 7/5 ; the diffusion spends more time at larger r. Thereby, the total timescale for particles released at r ∼ r snow to migrate to the annealing line would be ∼ t diff,snow , which agrees with the simulation result within a factor of 2. At t ∼ 1.0 t diff,snow (panel (b)), roughly 5 % of crystalline silicate particles come back to the snow line and the crystalline silicate surface density increases to a half of that in the equilibrium state at r ∼ r snow . The timescale to diffuse out from r ∼ r annl to r ∼ r snow is also ∼ t diff,snow , which is consistent with the simulation result in panel b. Although, in fact, the diffusion timescale for inward transport would be shorter than outward transport by inward disk accretion, accretion timescale may be comparable to diffusion timescale for transport from the snow line and the annealing line and the simple estimation with diffusion timescale helps to roughly interpret the simulation results. At t ∼ 5 t diff,snow (panel c), the crystalline silicate surface density reaches the equilibrium value at r ∼ r snow , but it has not reached the equilibrium value in outer region, because t diff is longer for larger r. Figure 4a shows the silicate particle surface density at t = 1.0×10 6 years of the run in Figure 3, where the equilibrium state is established. The orange and blue lines are the surface densities of crystalline and amorphous silicate particles calculated by Eq. (30) from the particle distributions obtained by the Monte Carlo simulations, In this plot, the blue dashed line represents the assumed gas surface density given by Eq. (57). The other dashed lines represent the analytically predicted surface densities of silicate particles, as explained below. These analytical predictions reproduce the simulation results very well. The magenta and olive dashed lines indicate the analytical zero net flux solutions that depart from the values of the steady accretion solution, given by Eq. (46), at the annealing and snow lines, respectively. The radial power-law index of the zero net flux solution (Σ sil ∝ r −q ) is given by q = p + 3/2 (Eq. (49)) where p = 3/5 and q = 21/10 in the viscous-heating dominated regime and p = 15/14 and q = 18/7 in the irradiation dominated regime (Eq. 57). The sum of numerically obtained silicate surface densities of the crystalline (the orange line) and amorphous (the blue line) agrees with the analytical steady accretion solution inside the snow line. Taking account of the decrease due to the transformation from amorphous to crystalline silicates at r ∼ r annl , the simulation result shows that the amorphous and crystalline silicate surface densities individually follow the steady accretion solution inside the snow and annealing lines, respectively. On the other hand, they fit with the zero net flux solution outside the snow and annealing lines, respectively. In Figure 4b, we adopt α = 10 −3 , so that t diff is 10 times longer and Σ sil has not reached the equilibrium state beyond the snow line. Therefore, the equilibrium surface densities of crystalline silicates Σ cry,sil and the sum of crystalline and amorphous silicates Σ tot,sil ≡ Σ cry,sil + Σ amr,sil obtained by the Monte Carlo simulations are fitted as and where Σ sil,0 (r) ≡ f sil F p/g Σ g (r) = 0.05 f sil 0.5 At r > r snow , because Σ cry,sil Σ amr,sil , Σ amr,sil (r) Σ sil,0 (r) r r snow These surface density distributions are physically interpreted as follows. We release amorphous silicate particles at the snow line. While the outward diffusion length increases with the square root of t, inward drift length due to disk gas accretion increases with t. Eventually, they balance to establish an equilibrium zero net flux distribution. After that, the total mass flux of the particles injected at different times that are passing the inner edge becomes equal to the given injection mass flux at the snow line (Ṁ sil = f sil F p/gṀg ). Thereby, amorphous silicate surface density follows the zero net flux solution outside the snow line and the steady accretion solution inside it. At the annealing line, crystalline silicate particles are generated from the drifting amorphous silicate particles with mass fluxṀ sil , which is similar to the injection of silicate particles at the snow line due to ice sublimation. Accordingly, crystalline silicate surface density follows the zero net flux solution outside the annealing line and the steady accretion solution inside it.
Equations (58) and (59) show that for r > r snow , independent of r, Σ cry,sil (r) Σ tot,sil (r) When the snow line is within the viscous heating dominated region, r annl /r snow 0.14. Equation (62) shows Σ cry,sil /Σ tot,sil 0.05. We note that the value of (r annl /r snow ) depends only on radial gradient of T . It is independent of other detailed disk parameters, such as α andṀ g , and pebble parameters, such as f sil , F p/g , and τ s,sil as long as τ s,sil < α. This robust, predicted value of the crystalline abundance in silicates in this simple case, Σ cry,sil /Σ tot,sil 0.05, is smaller than the observationally inferred abundance, ∼ 0.3 in Sitko et al. (2011) and ∼ 0.1-0.6 in Shinnaka et al. (2018). It strongly suggests that the effects of sticking to drifting pebbles and time-dependent F p/g must be incorporated into the simulations to predict the crystalline abundance in comets.

4.2.
Case with sticking to icy pebbles Figure 5 shows the results with the effect of sticking to pebbles after the establishment of an equilibrium state. Comparing to the simple case without the sticking effect (section 4.1), the decay of the silicate surface densities with r is faster outside the snow line, by the effect of rapid inward drift with τ s,peb after the sticking. The radial gradient of the silicate surface densities outside the snow line are analytically predicted by Eq. (56) as q 6.4 for α = 10 −2 and q 7.7 for α = 10 −3 in the viscous-heating dominated regime. In the irradiation dominated regime, q is larger by 15/14 − 3/5 0.47. The numerical results agree with the analytical prediction. Silicate particles have a steeper radial gradient than that without the effect of the sticking to pebbles (the power-law index is larger by ∼ 4 − 6). However, because we assume that the sticking probability is the same for crystalline and amorphous silicates, their surface density distributions are smaller by the same factor as in the non-sticking case, and Σ cry,sil /Σ tot,sil is still given by Eq. (62) as ∼ 0.05.
So far, we have considered a time-independent pebble mass flux, equivalently, a time-independent sticking probability. However, the pebble flux should start decaying after the pebble formation front reaches the characteristic disk size and the solid materials in the disk are reduced. Because the sticking effect substantially modifies the silicate surface density outside the snow line, the time dependence of the pebble flux changes the Σ cry,sil /Σ tot,sil ratio, which we investigate below. Sato et al. (2016) and Ida et al. (2019) showed that icy pebble flux decays rapidly after the pebble formation front reaches the characteristic disk radius (r disk ) at t 2.0 × 10 5 (r disk /100 au) 3/2 years, in the fast limit of pebble growth with perfect accretion. According to this result, we perform simulations in which the pebble flux starts decaying from the steady distribution in section 4.2. We set the decaying pebble flux as

Case with decaying pebble flux
where t 1 = t − t 0 is the time from which the decay is switched on at t = t 0 and t peb is the decay timescale with a typical value of a few ×10 5 years in the case of r disk ∼ 30 au (Sato et al. 2016;Ida et al. 2019). We also perform a simulation with t peb = 0 as an extreme case, which means sudden truncation of pebble flux. To keep the mass of individual super-particles unchanged, we inject particles with a time interval, δt inj = δt inj,0 exp(t 1 /t peb ). As shown in section 4.1, the typical timescale for crystalline silicates to return to the snow line, after they are released at the snow line and pass the annealing line, is Most of the amorphous silicates beyond the snow line may directly diffuse out from the snow line. It is expected that crystalline silicates beyond the snow line were released earlier on average by ∼ t diff,snow than amorphous silicates there. Because the crystalline silicate distribution would reflect an earlier pebble flux, Σ cry,sil /Σ tot,sil for decaying F p/g would be higher than the ratio in the case of the time-independent F p/g (Eq. (62)). According to the above prediction, we perform simulations with various values of t peb relative to t diff,snow , from the equilibrium states with F p/g = 0.1. Figure 6 show the distributions of Σ cry,sil and Σ amr,sil at t 1 = 2.5 t diff,snow of (a) t peb = 0, (b) 1.8 × 10 4 years (∼ 0.5 t diff,snow ), (c) 3.6 × 10 4 years (∼ t diff,snow ), and (d) 7.2 × 10 4 years (∼ 2 t diff,snow ), for α = 10 −2 andṀ g = 10 −7 M yr −1 . We also show (e) the initial equilibrium state, which corresponds to t peb = ∞. The ratio Σ cry,sil /Σ amr,sil is almost independent of r at r > r snow in all the panels, although Σ cry,sil and Σ amr,sil individually decay according to the decay of the pebble mass flux. Because the pebbles capturing the silicate particles also decay, both the crystalline and amorphous silicates extend to outer regions with time. The asymptotic value of Σ cry,sil /Σ tot,sil is (a) 0.25, (b) 0.20, (c) 0.13, (d) 0.10, and (e) 0.05. It increases with the decrease in t peb /t diff,snow , because the pebble fluxes reflected by the crystalline and amorphous silicate distributions are more different for more rapid pebble flux decay with smaller t peb . Because the diffusion effect smooths out the contrast in the pebble fluxes, Σ cry,sil /Σ tot,sil is saturated at ∼ 0.25 even in the limit of t peb = 0.
The observationally inferred crystalline abundance of comets is ∼ 0.3 in Sitko et al. (2011) and ∼ 0.1-0.6 in Shinnaka et al. (2018) (section1). Figure 6 show that our simulations predict the abundance of crystalline silicates is ∼ 0.2-0.25 for t peb 0.5 t diff,snow , while it is 0.1 for t peb 2 t diff,snow . The relatively fast decay with t peb 0.5 t diff,snow produces the results that may be consistent with the observed abundance. Ida et al. (2019) showed that t peb is a few ×10 5 years for a relatively compact disk with r disk ∼ 30 au, which can be comparable to ∼ 0.5 t diff,snow for α ∼ 10 −3 (Eq. (64)). The consistency with the observed data is discussed more in the next section.

CONCLUSION AND DISCUSSION
Crystalline silicates found in comets should have been formed in the disk inner region. One of the possible transfer mechanisms of the crystalline silicates to the disk outer region where cometary cores were formed is radial diffusion of silicate dust particles due to disk gas turbulence (Ciesla 2011). However, Ciesla (2011) did not present a quantitative estimate of the crystalline abundance for the comparison with the observations of comets has not been done, because the supply mechanism of amorphous and crystalline silicate particles during planet formation was not specified. Pavlyuchenkov & Dullemond (2007) and Arakawa et al. (2021) discussed an equilibrium radial distribution of crystalline materials. Because they assumed that amorphous silicates have a stationary distribution with the uniform solid-to-gas ratio in the steady accretion disk, they found that crystalline abundance is generally very low well beyond the snow line. However, the radial distribution of amorphous silicates must be derived by planet formation model for a qualitative estimate of crystalline abundance.
In this paper, adopting the "pebble accretion" model, we have performed Monte Carlo simulations, taking into account the release of silicate dust particles from pebbles at the snow line at T ∼ 170 K, transformation of amorphous silicate particles to crystalline particles at the annealing line at T ∼ 1000 K, sticking of silicate particles onto drifting pebbles beyond the snow line, and attenuation of the pebble flux due to consumption of solid material reservoir in the disk outer region. For the gas disk, a steady accretion model with viscous heating and irradiation from the host star is assumed. With this setting, we can quantitatively calculate the abundance of crystalline materials in all the silicate dust particles beyond the snow line.
In the simple case without sticking silicate particles onto drifting pebbles and with a steady pebble accretion, we found through the Monte Carlo simulation and analytical argument that the surface density of crystalline silicates scaled by the total surface density of crystalline and amorphous silicates is given by (r annl /r snow ) 3/2 uniformly beyond the snow line, independent of disk parameters and silicate particle size, as long as α > τ s (where α is the viscosity parameter and St is Stokes number of the silicate particles). The validity of the condition, α > τ s , is discussed in Appendix. When the viscous heating dominates at r ∼ r snow , r annl /r snow 0.14 and the crystalline abundance is 5 %. Although this robust value is substantially lower than the observationally inferred values, we found that with a more realistic condition where the sticking to icy pebbles and the pebble flux decay are included, the crystalline silicate abundance rises up to 20-25%.
Our simulation shows that the crystalline silicate abundance outside the snow line is independent of distance from the central star, r. This result is different from the results by Pavlyuchenkov & Dullemond (2007) and Arakawa et al. (2021), which showed the crystalline abundance decrease as r becomes larger. This difference results from the initial distribution of amorphous silicates. Because amorphous silicates are released from sublimating pebbles at the snow line, the radial gradient of the surface density is the same for both crystalline and amorphous silicates, in the framework of the pebble accretion scenario. The observationally inferred abundance has dispersion among comets of 10-60% (e.g., Shinnaka et al. 2018), although observations would include some uncertainty. There may be other factors that change the abundance, such as the evolution of a disk.
We note that annealing can proceed even if T < 1000 K (e.g., Yamamoto & Tachibana 2018;Ciesla 2011, and references therein), while we assumed that amorphous silicates are annealed immediately at T = 1000 K for simplicity. Whether or not they are annealed is determined by how long they experienced individual temperature. It could increase the crystalline silicate abundance from that obtained with the simple treatment at T = 1000 K, as well as increase the dispersion among comets. Because we track the particle motions including random walks due to gas turbulence, we can evaluate the progress of annealing for individual particles, which is left for our next paper.
In this paper, we assume a simple steady accretion disk, in order to highlight the physics associated with the pebble accretion model. In an early disk expansion phase, the snow and annealing lines are farther from the host star and crystalline silicates can be dragged to outer regions by the expanding disk gas (Ciesla 2011). However, if the dust particle diffusion timescale is shorter than the disk evolution timescale, the distribution of crystalline silicates would become the same as that in a steady accretion disk. The α value in the local turbulent diffusion can be much smaller than the effective α value for global disk gas accretion (angular momentum transfer) (e.g., Armitage et al. 2013;Hasegawa et al. 2017). The simulations coupled with the disk evolution are left to future study. Ida et al. (2021) considered the disk model with two different effective α parameters, one for disk accretion (α acc ) and the other for turbulent diffusion (α D ), because it is possible that the disk accretion is regulated by disk wind angular momentum removal rather than turbulent-diffusion angular momentum transfer. They showed that small silicate dust particles released from sublimating icy pebbles undergo gravitational collapse into rocky planetesimals just inside the snow line, if α acc α D and pebble mass flux is relatively high. In the case of α acc α D , crystalline silicate particles hardly diffuse out to regions much beyond the snow line. Relatively active disk wind and relatively high pebble mass flux would be realized in an early disk evolution phase, implying that rocky planetesimal formation by this mechanism is preferred to occur in the early phase. It could correspond to the formation of parent bodies of iron/stony meteorites. A later phase, in which the disk wind decays (α acc ∼ α D ) and the pebble mass flux decreases, may correspond to the situation in Section 4.3. In this case, the rocky planetesimal formation does not occur with Ida et al. (2021)'s mechanism. If rocky protoplanets in the terrestrial planet region have not fully grown beyond the asteroid parent body size, accretion of small dust particles is inefficient (e.g., Guillot et al. 2014) and the crystalline silicate particles diffuse out to regions much beyond the snow line. If cometary cores are formed through icy planetesimal formation in this timing, the cometary cores can include crystalline particles up to 20-25%. More discussions on a consistent scenario for rocky and icy planetesimal formation are left for future work.
We thank Takafumi Ootsubo, Hideyo Kawakita, Aki Takigawa, and Shogo Tachibana for fruitful and helpful discussions about observations and cosmo-chemistry. We also thank the referee for the helpful comments. This work was supported by JSPS Kakenhi 21H04512. α. Silicate particles can grow to some degree and accordingly their Stokes number can increase, after the release from sublimating icy pebbles. However, the results here do not change, as long as diffusion is dominated over advection due to gas drag for silicate dust particles (τ s < α). Here we estimate the maximum Stokes number of silicate particles, when the fragmentation limit is applied. As explained below, τ s is likely to be safely smaller than α for a steady accretion disk with our fiducial parameters, α = 10 −2 anḋ M g = 10 −7 M /yr.
The relative velocity v rel is regulated by turbulent diffusion for α 10 −3 (e.g., Sato et al. 2016). In this case, v rel √ 3α τ s c s .
Because τ s increases with the particle physical radius R d , the particle growth is limited when v rel exceeds a threshold fragmentation velocity v frag . The maximum Stokes number for the fragmentation limit is given by v rel v frag as τ s,max 1 3α where we scale v frag by a typical value of 1m/s for silicate-silicate collisions (e.g., Wada et al. 2013). In the region near or inside the snow line that we are mainly concerned, viscous heating is dominant (see section 2.1). Substituting T in the viscous heating region (Eq. (1)) into the above equation, we obtain τ s,max 3.1 × 10 −5 α 10 −2 −4/5 v frag 1 m/s 2 Ṁ g 10 −7 M /yr −2/5 r 1 au 9/10 . (A3) For the fiducial parameters, α = 10 −2 andṀ g = 10 −7 M /yr, in our disk model, τ s,max 1 × 10 −4 inside the snow line. Even if α = 10 −3 is used, τ s,max is still smaller than α. Therefore, our assumption that we simply set τ s = 10 −5 for silicate particles is justified. We refer to this estimate in sections 1, 2.2.1, 4.1, and 5. IfṀ g = 10 −8 M /yr and α = 10 −3 are assumed, τ s,max 5 × 10 −4 (r/1 au) 9/10 . We note that this τ s,max is ∼ α at a few au and our assumption of diffusion-dominance is marginal.
The particle physical radius R d corresponding to τ s,max is as follows. Stokes number of a particle with the bulk density ρ bulk which follows Epstein's law is given by where ρ g , Σ g , and H g are the spatial and surface densities of the disk and the pressure scale height, respectively. Substituting τ s,max into this equation, the typical particle size with the fragmentation limit of v frag = 1m/s in the viscous heating dominated regime is estimated as In outer disk regions, where pebbles are formed, irradiation would be dominated. In this case, τ s,max ∝ T −1 ∝ r 3/7 and Σ g ∝ r −15/14 (section 2.1). Accordingly, R d ∝ r 3/7−15/14 = r −9/14 . It means that silicate particles formed in outer regions to be embedded in pebbles are much smaller than that at the snow line. After their release at the snow line, they can grow up to the size given in Eq. (A6). The size released at a few au is ∼ 0.05 mm for α 10 −2 and ∼ 3 mm for α 10 −3 .