© The Authors 2025

1 Introduction

The outer disk is a region of interest as multiple chemical species have been detected at different layers in the disk (Öberg et al. 2021). Icy volatiles transported on dust grains in proto- planetary disks (Henning & Semenov 2013; Krijt et al. 2023) constitute a significant portion of the ice budget in protostellar cores, accounting for approximately 26–60% of the total elemental oxygen along with H2O and CO2, whereas it is 14–27% for elemental C (Boogert et al. 2015). The James Webb Space Telescope (JWST) has provided a detailed and spatially resolved ice inventory, which has helped to trace the changes in ice composition over time and through different environmental conditions. A recent study by Sturm et al. (2023) shows that the prominent absorption features in CO₂ and CO ices are shifted to pure ice absorption features, indicative that ices are exposed to energetic processing at low temperatures.

Ices formed in the interstellar medium (ISM) carry a heavy isotope imprint in H, C, O, and N compared to the solar/ISM values that make up the solar gas disk and are crucial for life as we know it (see review by Broadley et al. 2022). This signature is seen in comets sampling pristine outer-disk ices and primitive meteorites. For example, the high D/H ratios in deuterated H₂O and deuterated H₂S (hydrogen sulfide) in comet 67P/Churyumov-Gerasimenko suggest a prestellar origin for some of its ices (Alexander et al. 2017). Other more complex CHO-, N-, and S-bearing molecules in comet 67P/Churyumov- Gerasimenko show correlated abundances with protostellar material (i.e., IRAS 16293-2422 B; Drozdovskaya et al. 2019). This suggests that volatiles in Sun-like systems may be preserved from the protostellar phase into planetesimals (Drozdovskaya et al. 2019). Primitive meteorites that formed closer to the Sun are subjected to more processing than comets and are estimated to contain 5–10% pristine ISM material, (Alexander et al. 2017) while comets contain about 10–30% (Mumma & Charnley 2011). This evidence suggests that significant amounts of ISM- like ice survived disk formation and evolution, becoming part of comets and, to a lesser extent, asteroids.

The formation and destruction reaction pathways of chemical species are driven by stellar radiation and thermal processes that further modify the ISM-like original composition of the disk. High energetic photons such as ultraviolet (UV) photons from stellar radiation can induce photodissociation reactions on ice (AB_ice + h_ν → A_ice + B_ice) and photodesorption reactions (AB_ice + hν → AB_gas) to trigger the formation of simple and complex organic molecules, which are incorporated into pebbles, planetesimals, and later planets. Volatiles and complex organic molecules are carried in particles that experience radial drift and vertical stirring due to disk instabilities, which ultimately influence the composition of planets.

Numerous previous models have examined the reprocessing of prestellar ices and have found that substantial reservoirs of prestellar ices are incorporated into the outer disk through infall accretion from the envelope (i.e., Visser et al. 2011; Yang et al. 2013; Hincelin et al. 2013; Drozdovskaya et al. 2014; Yoneda et al. 2016). More recently, Bergner & Ciesla (2021) used a constant α value to represent gas turbulence within the disk and made a model where individual pebbles were covered by pristine ice mantles. They found that ice destruction is typically ineffective, except for small grains (<100 µm). This statement supports the scenario that icy pebbles (>100 µm) originate in the outer disk and subsequently migrate inward, contributing to the formation of cometary bodies.

Turbulence in protoplanetary disks is crucial regarding angular momentum transport and dust dynamics, which can affect planet formation. Beyond 0.1 au, where the magnetorotational instability (MRI; Balbus & Hawley 1991) is suppressed due to low ionization rates (Turner et al. 2014; Dzyurkevich et al. 2013), turbulence can arise from hydrodynamic instabilities, for example the vertical shear instability (VSI; Nelson et al. 2013). The VSI is driven by vertical differential rotation in baroclinic disks (Stoll et al. 2017; Barker & Latter 2015), and the turbulence intensity increases with the aspect ratio of the disk on rapid cooling timescales (Manger et al. 2020, 2021). It causes anisotropic turbulence with Z and ϕ stress-to-pressure ratio typically reaching magnitudes of 10⁻², while the R and ϕ stress-to-pressure ratio range from 10⁻⁴ to 10⁻³ (Flock et al. 2017, 2020), and large- scale gas circulations, which can stir up dust particles and create vortices and zonal flows (Lyra & Umurhan 2019).

Global simulations that have incorporated Lagrangian particles of millimeter-sized dust demonstrated that these particle sizes can be lifted close to the gas scale height above the midplane, due to strong vertical turbulence (Stoll et al. 2017; Flock et al. 2017, 2020; Dullemond et al. 2022). High-resolution simulations indicate that strong shear between neighboring bands of these meridional circulations can generate vortices in the R – Z plane (Flores-Rivera et al. 2020; Melon Fuksman et al. 2024a,b) as well and in the Z – ϕ plane (Manger et al. 2020, 2021) that are related to the Kelvin-Helmholtz instability and could serve as dust trap locations, enhancing local dust concentrations. The corrugation modes of the VSI can be translated to rings with azimuthal structures in the disk that can potentially explain dust substructures observed by the Atacama Large Millimeter Array (ALMA) (Blanco et al. 2021).

The radial and vertical transport by the VSI within the disk is expected to impact the dust particle trajectory significantly, and consequently its ice composition. For this reason, this work is dedicated to investigating the gas and dust dynamics and providing predictions on the ice chemistry timescales in the outer disk region. Turbulence can lift and mix the dust particles at different vertical layers, exposing them to vastly different thermal and UV environments (Semenov & Wiebe 2011).

The outline of the paper is as follows. In Sect. 2 we describe the numerical methods for modeling accretion disks that leads to the generation of the VSI modes. In Sect. 3 we describe the physical properties used in the RADMC3D code. In Sect. 4 we describe the thermal and photochemical reactions rates used to calculate the CO ice destruction timescale. Sects. 5 and 6 we present the results and discussion, respectively. Finally, Sect. 7 we present our conclusions.

2 Numerical methods

We used the PLUTO code¹ (Mignone 2007) to conduct the simulations of the gas in the protoplanetary disk. Our HD simulations are described based on the conservation of mass and momentum $$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho v)=0,$$(1) $$\frac{\partial (\rho v)}{\partial t}+\nabla \cdot (\rho vv)+\nabla P=-\rho \nabla \text{Φ},$$(2)

where ρ is the mass density, v is the velocity vector, ρv is the momentum density vector, and P is the gas pressure. We select the isothermal equation of state with $P=c_{\text{s}}^2\rho$, where c_s is the isothermal sound speed. The simulations are in 2.5D (two dimensions, three components) using spherical coordinates (r, θ, ϕ) with axisymmetry in the azimuth. The grid cells are set up with a logarithmic increase in the radial domain and a uniform spacing in the meridional domain. The radial domain ranges from 20 to 75 au and the vertical domain covers approximately ±5 gas scale height.

The initial gas density profile in cylindrical coordinates (R, Z) derived by Nelson et al. (2013) is $$\rho (R,Z)={\rho }_0{\left(\frac{R}{R_0}\right)}^p\exp \left[\left(\frac{G{M}_{*}}{c_{\text{s}}^2}\left[\frac{1}{\sqrt{R^2+Z^2}}-\frac{1}{R}\right]\right)\right],$$(3) $$\text{Ω}(R,Z)={\text{Ω}}_{\text{k}}{\left[(p+q){\left(\frac{H_{\text{R}}}{R}\right)}^2+(1+q)-\frac{qR}{\sqrt{R^2+Z^2}}\right]}^{1/2},$$(4)

where G is the gravitational constant, M_* is the mass of the star, and p and q are the radial profile exponent of the temperature and density, respectively. The initial density at the midplane is defined as, ${\rho }_0=\frac{{\mathrm{\Sigma }}_0}{\sqrt{2\pi }{H}_0{R}_0}$ where Σ₀ is the initial surface density at R0 is the reference radius, H₀ is the reference scale height (see Table 1). The Keplerian frequency is ${\text{Ω}}_{\text{k}}=\sqrt{G{M}_{*}/R^3}$. The gas scale height of the gas is defined as $$H_{\text{R}}=H_0{\left(\frac{R}{R_0}\right)}^{(q+3)/2}.$$(5)

The reference scale height ratio is H₀ = 0.1 at R₀ (Table 1). The temperature profile is $T\propto c_{\text{s}}^2$ is the sound speed and is given based on Equation (5) so that c_s= H_RΩ_k. The scale length of the vertical modes is proportional to H₀. The vertical component of the velocity in cylindrical coordinates, v_Z = v_r cos(θ) – v_θ sin(θ), is used for the analysis of the VSI (Figures 1 and 2).

Fig. 1 Dust particles of 100 µm (gray) being stirred by disk turbulence high above the disk midplane, compared to the more settled 1 mm particles (black). The disk turbulence is driven by the VSI represented by the vertical gas velocities (vz) expressed in terms of the local sound speed. This simulation snapshot at 300 orbits represents a disk time evolution for 0.1 Myr at 50 AU.

2.1 Dust

In this paper, we are interested in how the dust particles are stirred up by the VSI turbulent motions (Figure 1) and whether they are processed in different vertical layers, closer to high UV flux regions in the outer disk. From the 2D axisymmetric model of the density and rotational velocity described in Sect. 2, we introduced two Lagrangian particles with sizes of 100 µm and 1 mm in the simulation, and their vertical dust distribution is described as a Gaussian profile. Initially, the particles are distributed within a scale height close to the gas scale height to capture the initial dust-settling phase before grains settle more in the midplane. The particle dynamics is described using a Lagrangian dust module (Flock & Mignone 2021) in the PLUTO code using a spherical coordinates system. These dust grains are assumed to be homogeneous, compact, and spherical with a dust density of, ρ_dust =1.7 g cm⁻³. We do not include the effects of the dust-to-gas feedback. We also assume that dust particles are small enough compared to the gas mean free path length so they can be treated in the Epstein regime (Epstein 1924). The stopping time is important as the interaction that dust particles have with gas is different depending on the dust size. The stopping time dictates the time that it takes for a dust particle to be slowed down by the gas flow and it is defined as $$t_s=\frac{\text{St}}{{\text{Ω}}_{\text{k}}}=\frac{a{\rho }_{\text{dust }}}{\rho (R,Z)H_{\text{R}}{\text{Ω}}_{\text{k}}}=\frac{a{\rho }_{\text{dust }}}{\rho (R,Z)c_s}$$(6)

where a is the dust particle size, which is kept constant throughout the dynamical evolution of the disk. The Stokes number can be rewritten in terms of the disk physical parameters: $$\text{St}=\frac{a{\rho }_{\text{dust }}\sqrt{2\pi }}{{\mathrm{\Sigma }}_0}\exp \left[\frac{Z^2}{2H_{\text{R}}^2}\right]{\left(\frac{R}{R_0}\right)}^{1/2}$$(7)

The dynamical trajectory of our particles is influenced by their sizes and the turbulent gas motion. For the midplane, where Z/H_R=0, exp[…]=1, and at Z/H_R=1, exp[…]~1.65, therefore, for one scale height variation the Stokes number or particle size changes by less than a factor of 2. For example, the Stokes number for a 0.1 cm (1 mm) particle size at the midplane at 50 au is St=0.00426, and at Z = 5 au, the St=0.0070. Similarly, the Stokes number is an order of magnitude lower for the smaller size of 0.01 cm (100 µm). This is further motivated if we consider the turbulent fragmentation velocity assumed to be 1 m s⁻¹ and α = 10⁻³ (a typical level of turbulence consistent with VSI), for which the fragmentation-limited Stokes number results in 0.0019 with a corresponding particle size of 0.045 cm at the midplane (see Equation (33) in Birnstiel 2024). This implies that particles are small enough to avoid destructive collisions, which enables growth to larger sizes through mechanisms such as sticking and coagulation sticking collisions until the fragmentation barrier is reached. We inject a total of 10⁶ Lagrangian particles in the simulation uniformly distributed over two large grain sizes, and uniformly distributed in the r – θ plane. The initial velocity of the particles is adjusted to align the local Keplerian velocity.

Fig. 2 Time evolution of the global kinetic energy of the gas at 50 au normalized with the kinetic energy from pure Keplerian velocity.

2.2 Dust mass distribution of large and small dust particles

We follow the approach in Ruge et al. (2016) to calculate the mass distribution between the small and large dust grains as seen in Figure 3. For this, we first use the dust-size distribution (Mathis et al. 1977): $$\frac{dn(a)}{da}\propto a^{-3.5}.$$(8)

We do not model the small grains as Lagrangian particles in the VSI simulations, as these will remain well coupled to the gas. We use grain sizes between 0.1 < a_small < 30 µm discretized in 10 bins. The number 10 is arbitrary as it is only used to determine the dust-to-gas mass fraction in small size bins. These dust particles are well coupled with the gas motion and are distributed uniformly in the radial and meridional domain. The mass of small dust, M_small,_d, can be determined as $$M_{\text{small},\text{ d}}=\frac{4}{3}\pi {\rho }_{\text{dust}}{\displaystyle \sum _{i=1}^{10}{a}_{\text{small},i}^3}{N}_i\left(a_i\right),$$(9)

where N(a_i) is the absolute number of particles with radius a_i. The N(a_i) can be expressed as a proportion relative to a reference number of particles. We select the total number of particles, N_max , with radius a_max : $$N_i\left(a_i\right)=N_{\max }{\left(\frac{a_i}{a_{\max }}\right)}^{-2.5}.$$(10)

Similarly the amount of mass for the large particles is described as $$M_{\text{large, d}}=\frac{4}{3}\pi {\rho }_{\text{dust }}{\displaystyle \sum _{g=1}^2{a}_{\text{large },g}^3}{N}_g\left(a_g\right),$$(11)

where $$N_g\left(a_g\right)=N_{\max }{\left(\frac{a_g}{a_{\max }}\right)}^{-2.5}.$$(12)

The total mass of small and large grains is $$M_{\text{large,d }}+M_{\text{small,d }}=\frac{4}{3}\pi {\rho }_{\text{dust }}\frac{N_{\max }}{a_{\max }}\left({\displaystyle \sum _{i=1}^{10}{a}_i^{0.5}}+{\displaystyle \sum _{g=1}^2{a}_g^{0.5}}\right).$$(13)

In the case of large dust grains, each individual test particle contains a certain number of dust particles, Ñ_𝑔 (a_𝑔). Finally, the swarm particle mass for each large particle is given by $${\tilde{N}}_g\left(a_g\right)=\frac{N_g\left(a_g\right)}{{10}^6},$$(14)

where 10⁶ is the total number of large particles. For a given distribution we obtain the dust-to-gas mass ratio of M_small,d/M_gas = 0.0003, M_100µm,d/M_gas = 0.0015, and M_1mm,d/M_gas = 0.0036.

Fig. 3 Dust density distribution of our model at 200 orbits. The large grains are 100 µm (green) and 1 mm (orange). The small grains spanning from 0.1 to 30 µm is plotted as a black shade in the disk and are assumed to be well coupled with the gas.

3 Radiative-transfer as post-processing

We use hydrodynamical simulations to obtain the evolution of the particle component due to settling, drift, and stirring caused by VSI turbulence. For simplicity, these simulations use a local isothermal temperature profile. Using RADMC-3D (Dullemond et al. 2012), we compute the local UV field and the dust temperature for each dust species to calculate ice processing on dust grains as a post-process step. We do not attempt to fit our physical modeling structure with actual observations. We rather use the stellar template as a model for the stellar properties required by RADMC3D to reproduce the dust temperature and radiation field profiles.

Our fiducial case is described using the TW Hydra stellar template, a well-known low-luminosity (L_∗ = 0.2 L_⊙ and L_UV ~ 10⁻³ L_⊙) T Tauri star with a M_* = 0.8 M_⊙, T_eff = 4110 K, and R_∗ = 1.04 R_⊙ (Andrews et al. 2011). We also attempted to do a test comparison with another source with a higher UV field to explore the photodesorption and photodissociation rates of particles subjected to ice processing. Note that by changing the stellar properties one ought to adjust the disk physical properties accordingly; however, performing a parameter study of the disk properties in conjunction with the stellar properties is out of the scope of this study. This second source is the high-luminosity star (L_∗ = 40 L_⊙ and L_UV ∼ 10⁻¹ L_⊙), AB Aur (see Appendix A for more details about its stellar properties). Both sources have different UV intensity fields. The total integrated UV luminosity for the high-luminosity case is two orders of magnitude greater than that of the low-luminosity case.

In order to determine the dust opacity, our small dust grain sizes, 0.1–30 µm, are assumed to be well coupled with the gas, sharing the same spatial distribution as the gas. The large grains of sizes 100 µm and 1 mm are, although, more settled in the disk, both particle sizes experience much stirring by the upward and downward motions of the VSI, clearly seen in Figures 1 and 3. We use the DIANA standard dust opacities that consider a mixture of 60% amorphous laboratory silicates, 15% amorphous carbon and 25% porosity by volume, well-mixed on small scales (Woitke et al. 2016). We performed another run using the DSHARP opacities and the temperature profile is about a factor of 0.94 lower than the temperature profile with the adopted DIANA opacities (see more Appendix B). We do not include ice-coated grains in our models, that could show characteristic opacity features at ~3 µm (i.e., Ossenkopf & Henning 1994). Here we focus on the analysis of the ice processing timescale compared to the dynamical timescale that the particles spend above certain vertical spatial height in the disk (see Figure 4).

The local radiation field is calculated at each grid point in our disk model. The fluxes are computed at 71 wavelength bins spanning 91–200 nm. Similarly as described in Bergner & Ciesla (2021), the wavelengths are chosen such that, given the unattenuated stellar spectrum, the total wavelength-summed photodissociation differs by less than 5% from the wavelength- integrated photodissociation rate for the molecules considered here. We do not consider the scattering effects on the propagation of Ly-α photons, which has been shown to show some enhancement deeper disk layers with the Ly-α regime (Bethell & Bergin 2011). For this reason, our treatment for the photodesorption and photodissociation rates may be underestimated here. The resulting temperature and radiation field profiles are shown later in Sect. 5.4.

4 Chemical reaction rates

Our chemical modeling framework does not include multiple icy monolayers in the mantle on the grains. However, we assume the particles are icy grains with a CO and H₂O composition reminiscent of protostellar conditions. We focus on the chemical reaction rates in s⁻¹ that quantify the ice destruction timescales of dust grains through thermal- and photodesorption as well as photodissociation. The dust surface composition is given by the adopted binding energies representative of the first generation molecule (see Bergner & Ciesla 2021).

The thermal desorption reaction rate from the ice surface is $$k_{\text{th},\text{des}}=\nu \exp \left[\frac{-E_{\text{des}}}{k_B{T}_{\text{dust}}}\right],$$(15)

where T_dust is the temperature of the small dust grain and E_des is the binding energy for CO which is taken to be 1150 K (Garrod 2013). We exclude to map the thermal desorption for H₂O as it only becomes important inside 20 au, even for the high-luminosity test case. The pre-exponential factor, ν, is the first-order characteristic vibrational frequency of species CO on the dust grain (see Hasegawa et al. 1992) and is defined as $$\nu =\sqrt{\frac{2N_s{k}_B{E}_{\text{des}}}{{\pi }^2{m}_{\text{CO}}}},$$(16)

where N_s = 10¹⁵ cm⁻² is the number of binding energy sites per surface area of the grain, k_B is the Boltzmann constant, and m_co = 4.65 × 10⁻²³g is the mass of CO. For CO, ν is on the order of just a few 10¹² s⁻¹.

We also calculate the photodesorption induced by UV photons: $$k_{\text{ph},\text{des}}=Y_{\text{pdes}}{\sigma }_{\text{des}}{F}_{\text{UV}}/2.$$(17)

Here ${\sigma }_{\text{des}}=\frac{\pi a^2}{N_b}=\frac{1}{N_s}={10}^{-15}{\text{cm}}^2$, is the so-called standard photodesorption geometric cross-section for the small grains (see Bergner & Ciesla 2021) and N_b = N_sπa² = 10⁶ when using N_s = 10¹⁵ cm⁻² as the number of binding sites per grain surface for a standard grain size of 0.1 µm (Hollenbach et al. 2009). The Y_pdes = 10⁻⁴ in units of molecules per grain per photon is the photodesorption yield assumed to be constant for CO and H₂O based on laboratory experiments (Visser et al. 2011; Öberg et al. 2009). A constant photodesorption yield is assumed because it is difficult to assign a molecule-specific yield as photon absorption happens under the ice layer (see discussion in Ruaud et al. 2016; Bertin et al. 2012; Muñoz Caro et al. 2010). The factor 1/2 reflects that the UV photons reach the dust grain from one direction as in Bergner & Ciesla (2021). The F_UV is the total UV flux integrated from 91 to 200 nm.

We include the photodissociation rate as $$k_{\text{ph},\text{dis}}={\displaystyle \sum _n{\overline{\sigma }}_{\text{pdiss},n}}{\overline{F}}_{\text{UV},n}/2,$$(18)

where ${\overline{\sigma }}_{\text{pdiss },n}$ is the average UV photodissociation cross-sections for the CO and H₂O molecule that is taken from the Leiden Observatory database² (Heays et al. 2017) and n is the wavelength bins sampled in our model. The cross-sections are scaled down to by a factor of 10 to approximate the ice phase rates since photodissociation products can be trapped in the surface and rapidly react to reform the parent molecule lowering the effective photodissociation cross-section compared to gas phase (i.e., Öberg 2016).

5 Results

5.1 Dynamical trajectory of particles

The particles analyzed are subject to turbulence induced by the VSI. The initial distribution of the large grains is close to the gas scale height to reproduce the initial dust-settling phase for young disks. To analyze their dynamical trajectory, around 5000 particles of both sizes were randomly selected, from a set initially localized within 40 au and 50 au radial distances from the star for all the vertical extents. We do not consider choosing particles close to the inner and outer boundaries to avoid any boundary condition effects. The trajectory of the particles reached a dynamical evolution time of 300 orbits, corresponding to a time evolution of 0.1 Myr at 50 au. The linear phase of the VSI starts off after t > 1 orbit from the initial radius at the inner boundary (Figure 1). Initially, all particles go under the process of dust settling toward the midplane until the VSI acts upon the particles at slightly later times, usually at the first ~30 or 40 orbits, depending on the initial radial location. During the linear phase of the VSI, the particles experience small kicks or radial and vertical displacements. After t = 40 orbits, the generation of the VSI modes continue to generate turbulence through the entire radial extension. At t =70–100 orbits, also called the non-linear phase, the VSI modes reach convergence and saturation in turbulence levels of α ≈ 10⁻³. The tracks of the particles at the non-linear VSI phase tell us that particles with particles size of 1 mm and 100 µm undergo significant vertical stirring when the particles are subject to turbulence in the disk (see Figure 4). We also tried to reproduce a map of the local kinetic energy and local 1 mm dust density over time and height at 50 AU (see Figure C.1). From the top panel in Figure C.1, we depict the common features in the local kinetic energy as vertical lanes, so-called finger modes, within the first few tens of orbits corresponding to the VSI linear phase. Then, the body modes appear during the merging of the finger modes and the local kinetic energy increases slowly until the VSI reaches the saturation phase at around 150 orbits. The local dust density follows the VSI modes. The dust scale height reaches equilibrium at 5 H_R (red-dashed line in Figure C.1 middle panel) after 250 orbits during the saturation phase. In this stage, the width of the dust density distribution is around 2.3 H_R.

Figure 4 shows the spatial track of three selected particle cases to illustrate different evolutions. Two 100 µm particles were chosen at different vertical extents (left and middle panels, respectively). A 1 mm particle is also shown being shielded close to the disk midplane (right panel). The gray dashed line corresponds to the Z/R = 0.2 layer (also called the UV layer, see UV radiation field map in Figure 8 for reference), defined to be the bottom-most layer of the UV field and supported by recent observations of the CN (1-0) emission layer that is found to be a good indicator of UV penetration in the MAPS IV/ALMA imaging and fitting data from Bergner et al. (2021). This Z/R = 0.2 layer is the reference layer we use to quantify the amount of particles subject to some degree of ice processing. A particle that crosses this UV layer at any dynamical evolution time is susceptible to some degree of UV radiation. Therefore, some of its pristine composition could be susceptible to a photodestruction process (more detailed analysis in Sect. 5.4). None of the 1 mm particles cross the Z/R = 0.2 layer in our drawn particle sample. The dynamical track frame for the 100 µm particle (left panel) shows that the particle crosses the UV layer during the non-linear phase of the VSI (seen in orange color in the colorbar from Figure 4), whereas a particle below the UV layer remains shielded from the UV radiation keeping, in principle, its pristine composition in the outer disk at least within 300 orbits. Since our analysis of the composition of the particle is limited, quantifying the number of times per orbit each particle crosses the UV layer gives a very clear understanding of the exposure time that the particles are subject to UV radiation and, therefore, more likely to photodestruction processes that affect their pristine composition (see more details in Sect. 5.4).

Fig. 4 Spatial distribution of the dynamical trajectory for three different particles. The top panels show the track of the particles and the colorbar represents the dynamical evolution in orbits. The dashed gray line marks the elevation of Z/R = 0.2, our definition of the bottom-most UV layer. The top left and top middle panels show the dynamical trajectory of two 100 µm size particles, one showing that it crosses the Z/R = 0.2 during the saturation phase of the VSI after ~100 orbits, and the other being shielded from UV radiation all the dynamical evolution of the disk. The top right panel shows the trajectory for a 1 mm particle being shielded in the disk all the time. Our results show that no 1 mm particle crosses Z/R = 0.2. At the beginning of the simulation, all the particles experience dust settling toward the midplane and inward radial drift. After ~40 orbits, VSI induces diffusion in the radial and vertical motion of the particles.

Fig. 5 Time evolution of the radial velocity in m s−1 for the same particle locations as in Figure 4. The red horizontal line in each panel represents the mean radial velocity of the particles after 40 orbits (from left to right): 33.6, 12.5, and −7.7 m s−1, respectively.

5.2 Quantification of particles crossing the UV layer

From the same drawn particle sample, which includes the 1 mm and 100 µm size particles, we quantify the number of particles crossing the Z/R = 0.2 layer. As described earlier, these particles were initially located within a radial range of 40 < R < 50 au. However, at 300 orbits, these particles would end in different radial locations that could go beyond the initial radial range condition. At 300 orbits, we found that the percentage of particles crossing the Z/R = 0.2 layer is about 21% (1120 particles from the random selection sample), and 79% of the particles are UV- shielded in the disk. We emphasize that none of the 1 mm size particles in the sample cross the Z/R = 0.2; therefore, all particles in the UV crossing category are 100 µm size particles. The 1 mm particles remain shielded in the disk for all the dynamical disk time considered here.

The occurrence rate at which particles cross the Z/R = 0.2 layer happens mainly during the non-linear phase of the VSI, after 150 orbits, when the instability reaches turbulence convergence (see Figure D.1 in Appendix D). This could suggest that even at later times, beyond 300 orbits, fewer particles can still be vertically lifted at Z/R = 0.2, but less likely as according to Figure C.1 the dust scale height reaches equilibrium with the gas. From this UV crossing classification, we cannot infer the degree of ice processing these particles are subject to. For this, we have to compare the total time that particles spent above the Z/R = 0.2 layer, which involves the calculation of the orbital period from their radial location in the disk, together with the occurrence rate at which they are above the Z/R = 0.2 layer (more details and analysis are described in Sect. 5.4).

5.3 Vertical and radial velocity of particles

To visualize and exploit the velocity profile of the particles that undergo settling, radial drift, and turbulence, Figure 5 shows the time evolution of the radial velocity of the same particles. The initial dust-settling phase is hard to visualize in Figure 4 because of the crowding of the same particle track points being very spatially close with respect to its initial location; however, in the zoom-in box below of each panel in Figure 4 shows that the 1 mm dust particles settle inward and something similarly happens to 100 µm size particles. After t > 40 orbits, when the VSI linearly grows in its turbulence strength, the particles experience radial and vertical diffusion over time. During 40 < t < 110 orbits, the average radial velocity of the particles are −3.9 (left), 5.0 (middle), and −7.7 m s⁻¹ (right). The negative sign in the velocity reflects the particle’s inward motion, while the positive velocity reflects the particle’s outward motion.

Beyond t > 110 orbits, the turbulence strength of the VSI is greater, injecting momentum into the particles, reflected as diffusion in the radial and vertical directions (Figure 5). The maximum particle diffusion is clearly seen in Figure 5, when the highest radial velocity, whether inward or outward, happens after t = 100 orbits. The average outward radial velocity of the 1 mm particles are 164.7 (left), 91.8 (middle), and 62 m s⁻¹ (right). The average inward radial velocities are −108.5 (left), and −76.5 m s⁻¹ for the 100 µm (middle), and −66.6 m s⁻¹ for the 1 mm particle (right). The highest radial velocity for the 100 µm size particle (left) occurs where the gas density is already very low, at Z ~ 12 au disk elevation. On average, the gas radial velocity initially that carries the momentum outward is similar to the dust radial velocity, and the difference between them is very minimal (about 0.01 m s⁻¹).

Figures 6 and 7 top and middle panels shows the vertical and radial velocity profile of both the gas (blue) and the swarm of 1 mm dust particles (orange) per orbit, respectively. The dust velocities are the mean velocity per orbit among multiple particles of 1 mm and 100 µm sizes. Initially, within the first 40 orbits, the particle starts settling toward the midplane. After that, the VSI starts to act upon the particles, and the radial and vertical back-and-forth motion of the particles follows closely the gas motion. The maximum radial and vertical velocity occurs after 150 orbits when the VSI reaches the saturation phase. At this phase, the particles experience vertical circling as upward and downward motions and short radial displacements that continue relatively constant until 300 orbits. The velocity profiles between the dust and gas look similar, however, Figures 6 and 7 bottom panels show the mean difference between the radial dust and gas velocity in red and the mean difference between the vertical dust and gas velocity in gray. The most significant difference between remains on the 1 mm particles, which is expected as these particles would have a greater de-coupling with the gas compared to the 100 µm particles. The difference in the mean vertical velocity is also greater than the mean radial velocity in both particle sizes.

Fig. 6 Velocity profiles of the particles. Top and middle: time evolution of the mean vertical and radial velocity of a swarm of 1 mm particles around the midplane. Bottom: mean difference between the vertical dust and gas velocity (top panel) in gray, and the radial dust and gas velocity (middle panel) in light red.

Fig. 7 Same as in Figure 6 but for a swarm of 100 µm particles.

5.4 Ice particle processing by thermal and by photo-destruction processes in the outer disk

To quantify the level of ice processing in the particles, we model the local temperature profile and the UV radiation field using RADMC3D, as described in Sect. 3. Figure 8 shows the physical structure of the low-luminosity case, L_* = 0.2 L_⊙. The left panel shows the temperature profile of the small grains, where grains of all sizes remain relatively cold in the midplane. Overplotted are three Z/R elevations for reference. The integrated UV luminosity, from 91 to 200 nm, is L_UV = 10⁻³ L_⊙ (Figure 8, right panel). The Z/R = 0.2 level is located already in a UV-active region corresponding to log₁₀(F_UV) ~ 8.0 photons s⁻¹ cm⁻², where icy particles would already be influenced by UV radiation if crossing this layer. We adopt the Z/R = 0.2 layer to be the bottom-most UV layer. As described in Sect. 3, we additionally explore the physical conditions of a high UV luminosity source. The UV field for the high-luminosity case penetrates deeper (see Figure A.2 for reference), reaching higher optical depth regions closer to the disk midplane than the L_* = 0.2 L_⊙ case. Although the UV layer at the early stages of disk formation is highly uncertain, the combination of high disk accretion and particles being less settled, variations of the UV layer for younger sources could be extended to slightly higher disk elevation than Z/R = 0.2.

Figure 9 (left) shows the CO ice thermal desorption timescale for the L_* = 0.2 L_⊙. The thermal desorption timescale, t_th,des, is given by $t_{\text{th},\text{des}}=k_{\text{th},\text{des}}^{-1}$ (Eq. (15)) for so-called first- order desorption of a mono-layer of ice (i.e., Minissale et al. 2022). The CO ice destruction time for the L_* = 0.2 L_⊙ case (Figure 9 in left panel) demonstrates that the CO snowline is located at around 30 au; inside this radius, the CO ice is thermally desorbed from the dust particle within a year. At Z/R = 0.2 no particles get processed by thermal desorption. However, thermal desorption seems to be an important ice destruction reaction at elevations higher than Z/R > 0.3 at around 50 au. For the high-luminosity case, as the temperature is considerably higher, the CO snowline is pushed outward at a radius beyond 70 au, which implies that all small and large particles get instantaneously desorbed at a more extended disk radius. For H₂O ice, thermal desorption is only effective inside 20 au, even for the high-luminosity case, which is out of the radial range we consider here. Thermal desorption for CO ice and H₂O ice play a minor role in most outer disk regions, particularly at Z/R = 0.2.

Figure 9 (middle) shows the photodesorption destruction timescales for CO ice in Myr. The photodesorption timescale, t_ph,des, is similarly set by $t_{\text{ph},\text{des}}=k_{\text{ph},\text{des}}^{-1}$ (Eq. (17)), corresponding to the first-order desorption timescale considering one active monolayer. The middle panel in Fig. 8 shows the CO photodesorption timescale in Myr. Since photodesorption rates considered are unvarying across CO ice and H₂O ice, meaning fixed photodesorption yield and geometrical cross-section, the H₂O ice destruction timescale is identical to that of CO ice. The photodesorption time at Z/R = 0.2 is way too long t_0.2,phdes ≫ 10⁴ years for the L_* = 0.2 L_⊙ case, demonstrating that this photodestruction process is inefficient for the 100 µm particles that cross the Z/R = 0.2 layer. Compared with the high-luminosity case, photodesorption is very effective (t_40,phdes ~ 41.7 years), removing all ices from the dust particle at Z/R = 0.2. Above Z/R = 0.4, the photodesorption time is rapid (~36.5 days), and we expect all pristine CO ice and H₂O ice to be completely photo desorbed from the dust particle.

Figure 9 (right) shows the photodissociation of CO ice in Myr. The CO ice photodissociation timescale at Z/R = 0.2 is t_0.2,phdiss~ 1737.8 years, which is much faster, thus more efficient than photodesorption, to destroy CO ice. For the high stellar luminosity case, all 100 µm particles that cross the Z/R = 0.2 get photodissociated at a shorter timescale (t_40,phdiss~ 3.98 years) than t_40,phdes. In the case of H₂O ice, t_0.2,phdiss at Z/R = 0.2 is ~ 1584.7 years, which is slightly more efficient than that of CO ice (Figure 10). It is important to note that both photodesorption and photodissociation reactions are ineffective in the midplane, even in cases of high stellar luminosity. This is because the UV photons already get absorbed by dust particles at Z/R = 0.2.

As for the small grain sample (< 100 µm), the ice destruction timescale at Z/R ≥ 0.2 is very efficient given that the aerodynamic coupling with the gas will keep the small particles continuously lofted, prolonging settling toward the midplane, and cycled in the upper low-density regions by turbulence. For this reason, ice survival in small particles is expected to be considerably inefficient as they will tend to spend a more significant orbital period than the larger particles in the UV layer.

5.5 Dynamical timescale of particles at Z/R=0.2

In Figure 11, we summarize the processing timescale of all the 1120 particles (100 µm) that cross the Z/R = 0.2 layer after the initial dust-settling (>40 orbits) to capture the radial and vertical excursion of the particles due to VSI. The left plot in Figure 11 refers to the ice processing timescale of the particles above the Z/R = 0.2 layer. The y-axis represents the total orbital time that particles spent above Z/R = 0.2, and the x-axis represents the last crossing time in simulation orbit when the particle last crossed the Z/R = 0.2. The last orbital time parameter is helpful to quantify the total number of times the particles have crossed the Z/R = 0.2 layer (see Figure D.1). Note that the shape of the particle distribution follows Figure D.1, the occurrence rate of particles crossing Z/R ≥ 0.2. A handful of 100 µm particles have their last Z/R = 0.2 crossing happening more beyond 150 orbits and before 250 orbits. After 250 orbits, fewer particles have their last crossing at Z/R = 0.2. This is because settling and inward radial drift occur throughout the entire disk time evolution; thus, the trajectory of the particles will tend to be more settled toward the midplane as the disk continues to evolve in time. We expect only a few to no particles later crossing the Z/R = 0.2 is expected as after 250 orbits, the dust scale height dynamics reaches equilibrium with the gas. If the number of particles in the statistical analysis is increased, and the simulation is run beyond 300 orbits, the occurrence rate of particles having their last cross at Z/R = 0.2 can increase, filling in the space after 250 orbits. The colorbar represents the total distance traveled by the 100 µm size particles, which cover a range from Z ~ 4 au to Z ~ 8.5 au. Since some particles cross the midplane, this total distance is calculated based on the average of the positive Zs and the negative Zs around the midplane for each particle track; the mean total distance peaks around 5.4 au.

The following calculation steps were considered regarding the ice processing timescale of the particles at Z/R = 0.2, t_lifted,0.2. First, we calculated the orbital period of the 1120 particles at every location when the particles cross the Z/R = 0.2 layer. This step requires considering the radial distance from the star at every particle radial track point. Because every time the particle crosses the Z/R = 0.2 has a different radius, we consider the averaged distance traveled by the particle above the Z/R = 0.2 layer. Second, we recovered the amount of times the particles cross the Z/R = 0.2 layer (see Figure D.1 in Appendix). Similarly to Figure D.1, it shows that at most 1120 particles cross the Z/R = 0.2 layer after 150 orbits of disk evolution time. Finally, we calculated the ice processing time of the particles in years by multiplying the last orbital time of the particle (from the first step) by the amount of times the particles cross the Z/R = 0.2 layer (from the second step). In other words, if the time the particles spent above the Z/R = 0.2 layer is greater than the CO ice and H₂O ice destruction timescale, as calculated in Sect. 5.4, then the particle is processed. The level of processing of the particles is thus quantified if the particle processing timescales at Z/R = 0.2, t_lifted,0.2, is greater than any CO ice or H₂O ice destruction timescale, t_destroy, summarized on the following condition: $$t_{\text{lifted },0.2}\ge t_{\text{destroy }}=\text{ processed }$$(19)

It is important to mention that all CO ice and H₂O ice destruction processes have a different destruction timescale at Z/R = 0.2. That is why in Figure 11 (left) we overplotted the CO ice (gray line) and H₂O ice (blue line) photodissociation timescale for the low-luminosity case, t_0.2,phdiss, and the photodesorption timescale at Z/R = 0.2 for the high-luminosity case, t_40,phdes, same for both molecules. The other timescales were not included because 1) the CO ice and H₂O ice photodesorption timescale for the low-luminosity case is way too long; thus, none of the 100 µm particles get photodesorbed, and 2) the CO ice and H₂O ice photodissociation timescales for the high- luminosity case is way too fast; less than 4 years all particles have a higher degree of processing at Z/R ≥ 0.2. Given the dependency of the orbital period on stellar mass, we also include the orbital period for the high-mass star shown in Figure 11 (left) as transparent markers. Note that their orbital timescales determine the ice processing degree in the particle. Therefore, when we say the particle is processed, we limit the discussion to whether the particle is fully processed or is partially processed. It can be intuitive that the faster the photodestruction process is, the more ice will be desorbed from the grain surface. However, this determination will require following the time-dependent chemistry, coupled with other chemical reactions and more species, to accurately determine the amount of CO ice and H₂O ice that remains on the surface after going through a long track in the disk. This is extremely expensive to do for many particles; however, in Sect. 6.5, we attempted to do this for one particle and generalize the fate for the other particles with the same ice processing timescale history. We expect the particles to be fully processed above the Z/R = 0.4 layer.

In summary, in Figure 11 the bar plot on the right shows the number of CO ice (black) and H₂O ice (light blue) particles that are processed by any thermal or photochemical processes, and those that, even though they cross the Z/R = 0.2 layer, do not get processed, for both stellar luminosity cases. The total of particles drawn in the sample was about 5340, which included both 1 mm and 100 µm size particles. The amount of particles crossing the Z/R = 0.2 layer is 1120. Thus, the other 4220 particles (mostly 1 mm) remain UV-shielded in the disk. For the low-luminosity case, only 16.7 and 17.2% of those 1120, respectively, particles experience some degree of CO ice and H₂O ice processing caused by photodissociation. We emphasize that all the icy particles processed in the low-luminosity case are 100 µm size particles. The remaining of 4.3 and 3.8% of the 1120 icy particles do not get photoprocessed, which is added up to the non-processed category, already encompassing the 79% that were already UV-shielded in the disk. No large particles get thermally processed beyond ~30 au for the low-luminosity case. However, inside this radius, all particles get thermally processed for all vertical extensions. On the other hand, for the high-luminosity case, all particles (including those that do not cross the Z/R = 0.2) get thermally photodesorbed for all disk radius and vertical layers for the same physical disk structure as the low-luminosity case. The photodesorption becomes comparable to the photodissociation where 21% of the CO ice and H₂O ice particles crossing the Z/R = 0.2 layer become processed, while below the Z/R = 0.2 are the remaining 79% of particles that are UV-shielded.

Fig. 8 Dust temperature and UV radiation field for small grains of the disk properties described in Sect. 3. The spectral shape is representative of a stellar template with a L* = 0.2 L⊙ and LUV = 10−3 L⊙. The UV flux is integrated from 91 to 200 nm.

Fig. 9 CO ice destruction timescale caused by thermal desorption (left,1/kth,des), photodesorption (middle, 1/kph,des), and photodissociation (right, 1/kph,diss) in units of Myr. Overplotted are the 100 µm particles (in gray) and the 1 mm size particles (in black) for a snapshot at 300 orbits.

Fig. 10 H2O ice destruction timescale caused by photodissociation (in units of Myr).

Fig. 11 Particle processing timescales, tlifted,0.2. Left: tlifted,0.2 compared to the orbital timescales of the 100 µm particles that cross the Z/R = 0.2 layer. The x-axis represents the last crossing time in simulation orbit when the particle last crossed the Z/R = 0.2. The horizontal dashed and solid lines represent the CO ice (gray) and H2O ice (light blue) photoprocessing timescale at Z/R = 0.2 for the low-luminosity case and the photodesorption timescale at Z/R = 0.2 for the high-luminosity test case. Particles above these photodestruction timescales are subject to ice processing. The other photodestruction timescales at Z/R = 0.2 lie outside the tlifted,0.2 range shown in this plot. The colorbar represents the total distance traveled by the 100 µm particles, where the mean total distance peaks around 5.4 au. Right: Bar plot showing the distribution of particles processed by thermal or photochemical processes vs. those that remain unprocessed despite crossing the Z/R = 0.2 layer, for both stellar luminosity cases. The sample size comprises 1 mm and 100 µm particles. In the low-luminosity scenario, 16.6 and 17.2% of the crossing particles (all 100 µm particles) undergo CO ice and H2O processing via photodissociation, respectively. Only 4.3 and 3.8% do not experience photoprocessing despite having crossed the Z/R = 0.2 layer, and are added to the non-processed category, which already includes 79% of the UV-shielded particles. For the high-luminosity case, thermal photodesorption affects all particles across all disk radii and vertical layers, whereas for photodesorption and photodissociation processes it is slightly higher than for the low-luminosity case.

6 Discussion

6.1 Model limitations

Our hydrodynamical simulations consider the effects of the time evolution turbulence generated during the formation and evolution of a protoplanetary accretion disk. We include the interaction that the gas exerts on the dust particles, but we do not consider the dust-to-gas feedback that significantly impacts the gas flow patterns and can modify local pressure gradients due to dust evolution (Schäfer & Johansen 2022). Thus, we leave that to future work. The equation of state of the gas dynamics is described in an isothermal configuration, meaning that the local temperature remains constant in the vertical direction. The gas cooling timescale in the outer disk is expected to be much less than one orbital period, making the isothermal configuration a good approximation to describe the gas dynamics in the outer disk. However, these two caveats on the treatment of the hydrodynamics are considered to be less realistic in regions very close to the midplane where the VSI is expected to be less active in the midplane primarily due to stabilizing effects of buoyancy, chemical composition (i.e., high molecular emis-sivity), dust opacity (i.e., high optical depths trap radiation that reduces the cooling efficiency), and specific thermal conditions (i.e., temperature from the central star) (Melon Fuksman et al. 2024a,b). Our analysis and discussion of the particle dynamics are descriptive of those particles that have not yet coagulated or significantly fragmented with other particles in the disk.

Our ice processing analysis starts by quantifying the radial and vertical excursion of the particles when exposed to time-dependent turbulence in the disk while maintaining computational efficiency good enough to run thousands of analyses of the trajectory of each particle. Our chemical analysis in CO ice and H₂O ice is in steady-state based on stellar templates for young Class II disks. As the strength level of the UV field in Class 0 and I sources is highly uncertain as well as the vertical height at which UV photons penetrate in the disk, we see that from Figure 8, the Z/R = 0.2 corresponds to the bottom-most layer UV photons could penetrate into. This is also supported by observational evidence of CN emission occupying Z/R = 0.2 elevations in some disks (more discussion in Sect. 6.2). During the Class I stage, disks are expected to be puffier, albeit the outflows clear some of it to enable UV-driven chemistry (Flores-Rivera et al. 2021, 2023). However, the disk can still be relatively embedded, meaning that the UV stellar photons might not reach the Z/R = 0.2 level as easily as it is during the later Class II stage. For this reason, we emphasize our ice processing percentages as upper limit values. We do not include the chemical analysis from X-rays and cosmic rays. Regarding the cosmic rays, Cleeves et al. (2013) found that magnetized winds in T-Tauri stars can significantly reduce the low-energy cosmic ray flux by ten orders of magnitude than the value used for the ISM. Stellar X-ray photons can play an important role by ionizing deeper in the surface layers, particularly when cosmic ray fluxes are reduced, which can increase the ice destruction timescales closer to the midplane.

6.2 Observational constraints of turbulence and photochemistry in outer disks

The vertical extent to which turbulence operates in young disks with ages <1 Myr is quite challenging to determine due to the high molecular extinction in highly envelope-embedded disks. The eDisk/ALMA large program has found that some of the sub-mm dusty disks have not yet settled in the midplane (Ohashi et al. 2023). For example, the IRAS 04302 source has a dust scale height of 6 au at around 100 au disk radius (Lin et al. 2023). The integrated intensity maps of these young disks show, in general, that Class 0/I disks are more geometrically thick, potentially shaped by strong turbulence. This is the opposite scenario when comparing with Class II disks where the sub-mm dust scale height can be more than a factor of ten lower than the values inferred in Class 0/I disks (Villenave et al. 2020) with turbulence levels as low as 5–10% of the local sound speed (Flaherty et al. 2018; Teague et al. 2018), except for IM Lup which shows evidence of strong turbulence at levels of 0.18-0.30 the local sound speed (Flaherty et al. 2024). However, model predictions are optimistic that young Class II disks (~0.3 Myr) can experience levels of turbulence of 10⁻³–10⁻⁴ caused by the VSI and be observationally constrained in the CO upper layers of disks (Barraza-Alfaro et al. 2021). This prediction still holds when a planet carves a gap in a VSI-active disk. Further high angular (less than 10 au) and spectral (~0.05 km s⁻¹ with ALMA Band 6) is required to constrain the aforementioned turbulence level (Barraza-Alfaro et al. 2024). The James Webb Space Telescope (JWST) can also constrain the spatial vertical extent of dust. Duchêne et al. (2024) estimated the thickness of the 890 µm distribution to be around 14 au (a thickness consistent with our 1 mm and 100 µm dust excursion heights lifted by the VSI).

The young sources in the eDisk program have systematically analyzed the chemical morphology of some gas molecules to trace various components of the protostellar system (i.e., Ohashi et al. 2023). No photochemical imprints have been detected from this sample, therefore, it is uncertain to determine the location of a UV layer. UV photochemistry has been widely studied in Class II T Tauri systems, where the molecular envelope has been mostly dispersed. In situ photochemistry depends on the UV radiation field in which surface layers in the outer disk can still be strongly influenced. However, the vertical extent to which photochemistry is still active is somewhat unclear. The CN/HCN ratio has been proposed to trace UV fields in planetary nebulae and in photodissociation regions (Cox et al. 1992; Fuente et al. 1995). Motivated by the cross-section dependence of both molecules, recent chemical models have proposed that the most important formation pathway for CN involves the vibrational excited H₂ that is primarily driven by UV photons, <105 nm (Visser et al. 2018; Cazzoletti et al. 2018). This framework successfully reproduces the CN ring feature detected in TW Hydra (Teague et al. 2016). Moreover, an excitation analysis of CN emission in TW Hydra (Teague et al. 2017) determined that CN is emitted from low-density, more UV-exposed elevated regions of the disk, consistent with a UV-driven formation pathway. Later, Bergner et al. (2021) conducted a toy model of the local UV field supported by observational evidence and showed that in the region where CN emission is detected, the formation pathway via the vibrational excited H₂ (91–105 nm) would occupy Z/R = 0.2–0.4 elevations whereas, beyond ~100 au disk radius, the HCN photodissociation over CN (HCN + hν → CN + H at 113–150 nm) is more likely to reproduce the observations. Another work by Benz et al. (2016) explored the role of the UV field based CH⁺/OH⁺ column density ratios to determine that the ratios, analyzed in twelve young sources observed with the Herschel Space Telescope, are reproduced for a FUV flux that is 2–400 times greater than that of the interstellar radiation field (ISRF).

The JWST has observed ices through spectral absorption features in the envelopes of very young sources. Brunken et al. (2024) interestingly reported absorption features that appear as two peaks in the splitting bending mode of a couple of Class 0 sources. These absorption features were decomposed into different Gaussian fits augmented with laboratory data to match the spectra. They found that a high fraction of pure ¹²CO₂ ice (at 80 K) in the ¹³CO₂ ice feature indicating ice segregation, where pure ¹²CO₂ separates from other mixed ices due to thermal processing. The concentration of pure ¹²CO₂ ice caused by this process was correlated for sources with high stellar luminosity (as high temperature means a higher segregation process), making the ¹³CO₂ ice a useful tracer of protostellar thermal processing and composition. Furthermore, a work by Sturm et al. (2023) found that spatially resolved CO ice has been detected at unexpected altitudes within the disk, reaching up to the uppermost layers where CO gas emission is present. The CO ice absorption feature was found to be asymmetric and broader than expected for pure CO ice, which suggests the possible existence of an additional component where CO ice is mixed with either CO₂ or H₂O ices, motivating future work to focus on the modeling of a more realistic ice structure (i.e., ice layers) with mixed icy components (see more in Bergner et al. 2024). In the future, it could be possible for the JWST to detect an ice tracer of photochemistry in young disks.

6.3 Comparison with other Lagrangian models in VSI-active disks

Flock et al. (2017) used a 3D VSI configuration and treated dust particles in the Epstein regime. In their setup, the Lagrangian dust particles are initially located in the midplane, and as the disk evolves in time, the particles are vertically lifted by the VSI. Our results show that if particles are located initially close to the gas scale height, 0.1 mm particles can be lifted to slightly higher elevations. In our case, the dust scale varies with time, following the oscillations of the VSI, and continues the periodic behavior until the oscillation height reaches equilibrium after ~250 orbits (Figure C.1) concentrated into a narrow layer of width ~2.3 H_R that oscillates around the midplane. In Figure 13 of Flock et al. (2017), the mix of the particles is instead more scattered around the midplane, which we attribute to a 3D effect.

Similarly, Stoll & Kley (2016) and Dullemond et al. (2022) also use a method to treat particles in the Epstein regime by solving the equation of motion of the particles, including the force of gravity and the friction with the gas. They do not follow the tracking of individual particles over time. However, both works showed that the VSI is very effective in isothermal and radiative configurations in stirring up dust particles even with a Stokes number of 0.1 to significant heights equal to the gas pressure scale height, which implies that substantial vertical dust distribution is expected to be seen in observations.

6.4 Composition of planetary bodies in the comet-forming zone

The comet-forming zone in the early solar nebula is defined to be in the outer parts of the disk (>5 au), beyond the snowline of volatile materials (Whipple 1972; Bockelée-Morvan et al. 2004). The smaller grains (<100 µm) typically undergo significant levels of ice chemical processing above Z/R > 0.2. As a result, these slow-growing particles experienced a combination of different degrees of ice processing that will eventually be imprinted in the planetary body as isotopic ratio variations (Bergner & Ciesla 2021). We stress that 100 µm grains are still circulated by turbulence to a lower extent in the UV region, where ice destruction can still play a role by removing some pristine ice from the long-standing orbital period that the dust grains spent above the UV layer. In contrast, larger 1 mm grains are very well shielded in the disk and transported inwards toward the star. The drifting of the particles in the midplane is expected to contain a mix of processed particles inherited from the small grains as they settle toward the midplane. The millimeter-size grains with a icy mantle composition that have not yet accreted into a planetary-forming body will eventually cross their respective snowline, losing all their pristine composition.

van Kooten et al. (2024) reported the discovery of organic-rich “dark clasts” that represent pebble-sized outer disk materials related to the accretion region of comets. Based on the analysis of their chemical composition, the inferred ice component of this pristine material contains high δ¹⁵N and low δD values compared to pristine ISM ices, which are best explained by vertical mixing of small dust grains (≤100 µm) in the outer disk. In detail, these icy grains access the ¹⁵N-rich surface layers of the disk that experience N2 photodissociation corresponding to significant N isotope fractionation (Chakraborty et al. 2014). At the same time, they equilibrate with D-poor H₂ gas in the disk (Cleeves et al. 2014), resulting in ¹⁵N-rich and D-poor ice layers coated onto ISM-derived organic matter. This suggests that even the more pristine material from the comet-forming region shows signs of ice processing in its composition. Similarly, small grains are expected to rapidly lose their ice after prolonged exposure in the upper layers. Therefore, it is unlikely that the small particles at these radii are the leading carriers of pristine ISM-like for ice inheritance. Furthermore, they found that the dark clasts have very ¹⁶O-poor bulk compositions that are similar to cometary Interplanetary Dust Particles (IDPs; Starkey et al. 2014; Nguyen et al. 2022), suggesting a large addition of outer disk derived ice relative to chondrites, which aligns with the idea of the CO self-shielding isotope effect in the outer disk (Lyons & Young 2005). We expect that icy bodies, such as comets, inherited mainly their reservoirs of pristine composition from the inward drift and accreting icy pebbles in the midplane. An important future consideration will include dust growth treatment in the simulations to assess the amount of pristine dust aggregates incorporated into the planetary body via streaming instability, SI (Lorek et al. 2018).

The D/H ratio in water has traditionally served as a tracer for the origin of water on Earth (see review by Ceccarelli et al. 2014). The D/H ratio of water referenced by the vienna standard mean ocean water (VSMOW) is 1.56×10⁻⁴ (i.e., Hallis et al. 2015), where asteroids have the closest D/H ratio relative to Earth followed by comets (Altwegg et al. 2015). Variations of the D/H value in comets depend on how much the water is preserved when turbulent transport brings the icy pebbles to the disk surface and photodissociates by UV radiation. Such a photodestruction process in water ice can significantly modify the D/H value in the solar nebula. The influence that water ice pho-todissociation can have on the D/H ratio would also depend on the presence of other relatively abundant molecules (i.e., atomic oxygen) that help reform water ice back in the disk. Furuya et al. (2013) found that water ice incorporated into the disk from the molecular cloud is destroyed by UV at the disk surface (at a column density close to 10²² cm⁻² around 30 au) while the atomic oxygen entering the deeper disk condenses on the grain surface and converted back to water ice through OH-ice as a reactant, which can increase the water-ice abundance beyond 40 au. In their model with vertical mixing, the water-ice abundance decreases in the midplane as it gets destroyed in the surface layers. Consequently, water-ice abundance and other deuterated ices (i.e., HDO-ice) decrease over time, and the D/H column density ratio decreases by up to one order of magnitude within 10⁶ years between 30 and 50 au. The photoproducts from water-ice, OH-ice, and H-ice can recombine again in the deeper layers with themselves or with H₂CO-ice to form water ice. The reformation of water ice is efficient as OH-ice should be formed mainly on grain surfaces to compensate for the destruction. In the context of planetesimal composition, the D/H ratio of cometary water has been proven to have originated in the solar nebula given the high D/H ratio detection by JWST (Slavicinska et al. 2024). However, the D/H ratio decreases due to photoprocessing as the disk evolves.

The presence of complex organic molecules in protoplanetary disks is crucial for understanding the inventory of prebiotic compounds and their potential survival through the formation and accretion processes that contribute to cometary bodies (Walsh et al. 2014). Our results indicate that UV photodissoci-ation influences the destruction efficiency of simpler molecules such as CO and water; similarly, we expect the same for complex organics. Bergner & Ciesla (2021) found that for complex organic molecules rather than HCN, which is the most efficiently photodissociated molecule at Ly-α wavelength, have less Ly-α photodissociation rate suggesting that they are more likely to survive the disk environment.

6.5 Toy model: Surface chemistry evolution for a single particle

We follow the surface chemistry of a 100 µm particle that crosses the Z/R = 0.2 layer to explore the impact of the particle’s trajectory in the disk on the CO_ice and H₂O_ice abundance evolution at the surface of the particle. We use the NAUTILUS gas-grain astrochemical code (Hersant et al. 2009; Ruaud et al. 2016; Gavino et al. 2021) that considers the rate equation approach based on Hasegawa et al. (1992). We used the kinetic database for astrochemistry network (kida.uva.2014³; Wakelam et al. 2015). The code includes thermal desorption, cosmic-ray-induced desorption, photodesorption, and chemical desorption. However, the surface reactions, as well as adsorption, are turned off. There are two reasons to justify this prescription. First, because we want to investigate the impact of the UV field only on the evolution of the CO_ice and H₂O_ice, we removed the possible surface reactions that can affect their abundances. Secondly, removing adsorption processes allows us to avoid potential unrealistic re-adsorption events of gas-phase CO and H₂0 when the particles return to regions with lower UV flux. These two prescriptions combined ensure that the CO_ice and H₂O_ice are only affected by the local physical conditions such that we can obtain quantitative information on the UV flux value from which the CO_ice and H₂O_ice photodissociation rate becomes significant.

The NAUTILUS code considers the physical environment of the particle evolving with time. The code reads the evolving physical parameters (gas density, UV flux, gas temperature, and dust temperature) the particle encountered during its trajectory in the disk. The molecular initial abundances are extracted after an evolution of 10⁶ years in a dark molecular cloud (10 K, n_gas = lO⁴, A_v=10). The initial abundance used for this molecular cloud simulation is atomic and is the same as those used by Gavino etal. (2023).

Figure 12 shows that both the CO_ice and H₂O_ice abundance remain constant as the particle is UV-shielded in the disk during the whole settling and linear-VSI phase. At ~46.4 Kyr (at 130 orbits at 45.8 au), the particle first enters the UV-active region (log₁₀(F_UV) = 7.0 photons s⁻¹ cm⁻²) in the saturation-turbulence phase of the non-linear modes of the VSI. However, the abundances only significantly starts to drop at ~50.9 Kyr (at 143 orbits at 46.5 au) when the particle reaches a UV value of 9.1 photons s⁻¹ cm⁻². In agreement with the particle timescale calculation in Sect. 5.4, the t_lifted,0.2 =6.1 Kyr whereas the UV photodissociation timescale for the L_* = 0.2 L_⊙ case is t_destroy = 1.3 Kyr which means that the particle is subject to ice processing altering the initial composition of the particle.

As mentioned above, some assumptions when reproducing the abundance profiles include the absence of chemical surface reactions. In particular, this excludes the major surface reactions involving CO and water. We performed an additional chemical simulation in which the particle’s surface is chemically active, meaning surface reactions without adsorption. We find, for instance, that the activated reaction OH + CO → CO₂ + H (Ruaud et al. 2015) as well as the reaction CO₂ + hν → CO + O, can significantly change the CO_ice evolution profile. In the case of H₂O_ice, it remains adsorbed in grain surface. However, going beyond in the analysis is out of the scope of this study. Adsorption was also excluded in the main simulation (Figure 12). To see the impact of adsorption on the evolution of CO_ice, we performed an additional chemical simulation where adsorption is active. In this case, the CO_ice abundance re-increases when the particle returns to the UV-shielded midplane. But since there is a depletion of the CO gas in the midplane, in principle there should be no CO adsorption in this region (crosshatched region in Figure 12). This means that a fully processed particle would reach the disk midplane and potentially get mixed with unprocessed icy pebbles.

Fig. 12 Relative abundance of the 100 µm particles in the first panel in Figure 4. The COice and H2Oice abundance is a function of orbital time of the particle (in years). During the initial dust-settling phase, the COice and H2Oice abundance remains constant as the particle remains UV-shielded in the disk until it starts experiencing the vertical circling motion during the VSI non-linear regime at a few tens of hundreds of orbits. The big drop in both abundance profiles happens at ~50.9 Kyr (at 143 orbits at 46.5 au), when the particle reaches the UV-active layer (log10(FUV) = 9.11 photons s−1 cm−2) in the saturation-turbulence VSI phase. Here, the particle loses a significant amount of COice and H2Oice from its surface. These icy volatiles are expected to be completely absent from the grain surface when it is back to the disk midplane (crosshatched region).

7 Conclusions

Our modeling framework investigates dust particle dynamics and ice processing within a turbulent gas scenario caused by the VSI in the disk. The dust particles are described as Lagrangian particles of 1 mm and 100 µm in size, experiencing self-consistent gas drag in the disk. The spatial tracks of different particles were examined at different vertical extents. The bottom-most UV layer at Z/R = 0.2 represents the region where the particles spend enough time to be influenced by UV and thermal ice processing despite being exposed to lower levels of UV radiation, which can still lead to a loss of pristine ice composition in the dust particles. The main points of our conclusions are summarized as follows:

1. In the scenario where the particles are initially at a certain height above the midplane, the VSI causes significant stirring to them within the disk, with particles of 100 µm in size frequently crossing the bottom-most UV layer, thereby exposing them to UV radiation that leads to an effective photodissociation. Over time, the percentage of 100 µm size particles crossing the UV layer increases, highlighting the impact of VSI on particle dynamics and differential exposure to UV radiation. This suggests that intermediate-sized grains in the upper layers of disks can experience strong replenishment by turbulence in all directions and are expected to be subjected to ice processing due to prolonged exposure to the UV layer. In contrast, larger particles (1 mm) remain shielded below this layer throughout the disk’s dynamical evolution, and serve as the primary carriers of pristine ice composition within the disk.

2. The degree of ice processing in particles within a protoplan-etary disk is proportional to the occurrence rate of particles crossing the UV layer. For low-luminosity stars (L_* = 0.2 L_⊙ and L_UV = 10⁻³ L_⊙), only 17% of particles that cross the UV layer experience effective photodissociation. When UV-processed particles return to the midplane, they mix with pristine icy pebbles, leading to a mixed composition.

3. The VSI significantly influences the radial and vertical momentum of dust particles in the disk. As VSI grows over time, it causes particles to experience increasing radial and vertical diffusion. Beyond ~100 orbits, turbulence significantly increases the particle’s momentum until it reaches a maximum when in equilibrium with the gas (at ~ 150 orbits). The resulting dynamical track of the particles is in circular motion, causing the particles to frequently cross the UV layer, prolonging their exposure to UV radiation.

We highlight how different particle sizes and their initial positions within the disk influence their dynamic behavior and susceptibility to UV-induced ice processing. Our findings suggest that stellar luminosity plays a crucial role in determining the extent of photochemical processing, impacting the evolution and composition of particles within the disk. A future project would follow the chemical analysis for organic and refractory elements. Our findings align with a previous study from Bergner & Ciesla (2021), which supports the formation of icy bodies, such as comets, via the accretion of drifting icy pebbles in the outer disk and the inclusion of highly processed material that make up their chemical composition. It further reinforces the role of icy pebbles as the leading volatile carriers, providing critical building blocks to planets where they may contribute to the conditions necessary for the emergence of life.

Acknowledgements

We thank the anonymous referee for providing supportive feedback. L.F-.R. and M.L. acknowledge this work is funded by the European Research Council (ERC Starting Grant 101041466-EXODOSS). S.G. acknowledges support from the Independent Research Fund Denmark (grant No. 0135-00123B). A.J. acknowledges funding from the Danish National Research Foundation (DNRF Chair Grant DNRF159) and the Carlsberg Foundation (Semper Ardens: Advance grant FIRSTATMO). M.F. has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement No. 757957). The Tycho supercomputer hosted at the SCIENCE HPC center at the University of Copenhagen and the astro-node9 cluster at the MPIA were used to support this work. L. Flores-Rivera thanks Troels Haugbølle for the excellent guidance toward using the Tycho cluster and Sebastian Marino for the useful help regarding some conceptual technicalities when using RADMC-3D. Thanks to Jennifer Bergner and Riohei Nakatani for the helpful discussion about some theoretical chemistry terminology. We thank the matplotlib team (Hunter 2007) for the good quality tool in order to better visualize the data.

Appendix A Stellar spectra template comparison

In addition to TW Hydra, our fiducial sample, we compared it with another case, AB Aur, to explore only the importance of a higher stellar UV field and temperature in the ice processing timescales. Figure A.1 shows both stellar templates for comparison. The stellar template of AB Aur, a Herbig Ae star, has different stellar properties compared to TW Hydra: M_* = 2.4 M_⊙, T_eff = 9770 K, and R_* = 2.37 R_⊙ (Woitke et al. 2016; Boc-caletti et al. 2020; Currie et al. 2022). We recognize that by changing the stellar parameters, the disk properties (i.e., disk temperature and density profiles) must be adjusted; however, doing a parameter study to explore how the chosen stellar properties influence the disk properties is out of the scope of this study.

Figure A.2 shows the UV radiation field of AB Aur, which we use to calculate the CO ice and water ice destruction timescale caused by photodesorption and photodissociation processes.

Fig. A.1 TW Hydra and AB Aur stellar spectra templates normalized at 1 pc used as an input parameter in RADMC3D.

Fig. A.2 TW Hydra and AB Aur stellar spectra templates normalized at 1 pc used as an input parameter in RADMC3D.

Appendix B Dust opacities

We compare the radial temperature profile in the midplane between the DSHARP and the DIANA opacities (Figure B.1). The DSHARP opacity follow previous work from Pollack et al. (1994) regarding protoplanetary disk composition that includes a mixture of water ice (20% by mass, Warren (1984)), astronomical silicates (Draine 2003), Troilite (Henning & Stognienko 1996), and refractory organics (Henning & Stognienko 1996) without porosity. The DIANA standard dust opacities consider a mixture of 60% amorphous laboratory silicates, 15% amorphous carbon and 25% porosity by volume, well-mixed on small scales (Woitke et al. 2016).

Fig. B.1 Radial temperature profile in the midplane comparison between DSHARP and DIANA opacities.

Appendix C Time evolution of the local kinetic energy and dust over height

We furthermore explore the time evolution of the local kinetic energy normalized with respect to the Keplerian velocity at the midplane, and local dust density at 50 au. In order to map the local kinetic energy and dust density in a more uniform spatial distribution, we transform the velocity vectors to cylindrical and then use a cubic method to interpolate the data on a new grid space with a Z direction ranging from −15 to 15 au at 50 au. As expected, the fast growth of the finger modes in the disk happens within the first few tens of orbits (see Figure C.1 top panel). The dark lane at the midplane shows v_z is nearly zero. Second energy saturation happens after ~150 orbits. The middle panel in Figure C.1 shows the time evolution of the 1 mm dust density distribution at 50 au. During the second energy saturation, we see particles already going upward and downward and in a circular motion. Notice that during the saturation phase, the particles reach a maximum in height and after 220 orbits y reaches equilibrium with the gas. The bottom panel shows the overlap between the dust density distribution (middle panel) and the local kinetic energy (top panel), both at 50 au.

Fig. C.1 Time evolution of the local kinetic energy (top) and local 1 mm dust distribution (middle) at 50 au. The red dashed line at Z = 5 HR represent the height when particles reach equilibrium with the gas. The bottom panel shows the middle panel overplotted on the top panel.

Appendix D Occurrence rate of particles above Z/R = 0.2

In other to calculate how long do particles spend above the Z/R = 0.2, we gathered how many times the particles cross this layer. Figure D.1 shows the amount of times (as a factor) the 100µm particles spent at Z/R ≥ 0.2.

Fig. D.1 Occurrence rate of particles above the Z/R = 0.2 layer. This occurrence rate represents the amount of times particles spent above Z/R = 0.2 layer. The x-axis shows time in orbits in which particles made their last UV cross.

References

All Tables

All Figures

	Fig. 1 Dust particles of 100 µm (gray) being stirred by disk turbulence high above the disk midplane, compared to the more settled 1 mm particles (black). The disk turbulence is driven by the VSI represented by the vertical gas velocities (vz) expressed in terms of the local sound speed. This simulation snapshot at 300 orbits represents a disk time evolution for 0.1 Myr at 50 AU.
In the text

	Fig. 2 Time evolution of the global kinetic energy of the gas at 50 au normalized with the kinetic energy from pure Keplerian velocity.
In the text

	Fig. 3 Dust density distribution of our model at 200 orbits. The large grains are 100 µm (green) and 1 mm (orange). The small grains spanning from 0.1 to 30 µm is plotted as a black shade in the disk and are assumed to be well coupled with the gas.
In the text

Fig. 4 Spatial distribution of the dynamical trajectory for three different particles. The top panels show the track of the particles and the colorbar represents the dynamical evolution in orbits. The dashed gray line marks the elevation of Z/R = 0.2, our definition of the bottom-most UV layer. The top left and top middle panels show the dynamical trajectory of two 100 µm size particles, one showing that it crosses the Z/R = 0.2 during the saturation phase of the VSI after ~100 orbits, and the other being shielded from UV radiation all the dynamical evolution of the disk. The top right panel shows the trajectory for a 1 mm particle being shielded in the disk all the time. Our results show that no 1 mm particle crosses Z/R = 0.2. At the beginning of the simulation, all the particles experience dust settling toward the midplane and inward radial drift. After ~40 orbits, VSI induces diffusion in the radial and vertical motion of the particles.

In the text

	Fig. 5 Time evolution of the radial velocity in m s−1 for the same particle locations as in Figure 4. The red horizontal line in each panel represents the mean radial velocity of the particles after 40 orbits (from left to right): 33.6, 12.5, and −7.7 m s−1, respectively.
In the text

	Fig. 6 Velocity profiles of the particles. Top and middle: time evolution of the mean vertical and radial velocity of a swarm of 1 mm particles around the midplane. Bottom: mean difference between the vertical dust and gas velocity (top panel) in gray, and the radial dust and gas velocity (middle panel) in light red.
In the text

	Fig. 7 Same as in Figure 6 but for a swarm of 100 µm particles.
In the text

	Fig. 8 Dust temperature and UV radiation field for small grains of the disk properties described in Sect. 3. The spectral shape is representative of a stellar template with a L* = 0.2 L⊙ and LUV = 10−3 L⊙. The UV flux is integrated from 91 to 200 nm.
In the text

	Fig. 9 CO ice destruction timescale caused by thermal desorption (left,1/kth,des), photodesorption (middle, 1/kph,des), and photodissociation (right, 1/kph,diss) in units of Myr. Overplotted are the 100 µm particles (in gray) and the 1 mm size particles (in black) for a snapshot at 300 orbits.
In the text

	Fig. 10 H2O ice destruction timescale caused by photodissociation (in units of Myr).
In the text

Fig. 11 Particle processing timescales, tlifted,0.2. Left: tlifted,0.2 compared to the orbital timescales of the 100 µm particles that cross the Z/R = 0.2 layer. The x-axis represents the last crossing time in simulation orbit when the particle last crossed the Z/R = 0.2. The horizontal dashed and solid lines represent the CO ice (gray) and H2O ice (light blue) photoprocessing timescale at Z/R = 0.2 for the low-luminosity case and the photodesorption timescale at Z/R = 0.2 for the high-luminosity test case. Particles above these photodestruction timescales are subject to ice processing. The other photodestruction timescales at Z/R = 0.2 lie outside the tlifted,0.2 range shown in this plot. The colorbar represents the total distance traveled by the 100 µm particles, where the mean total distance peaks around 5.4 au. Right: Bar plot showing the distribution of particles processed by thermal or photochemical processes vs. those that remain unprocessed despite crossing the Z/R = 0.2 layer, for both stellar luminosity cases. The sample size comprises 1 mm and 100 µm particles. In the low-luminosity scenario, 16.6 and 17.2% of the crossing particles (all 100 µm particles) undergo CO ice and H2O processing via photodissociation, respectively. Only 4.3 and 3.8% do not experience photoprocessing despite having crossed the Z/R = 0.2 layer, and are added to the non-processed category, which already includes 79% of the UV-shielded particles. For the high-luminosity case, thermal photodesorption affects all particles across all disk radii and vertical layers, whereas for photodesorption and photodissociation processes it is slightly higher than for the low-luminosity case.

In the text

Fig. 12 Relative abundance of the 100 µm particles in the first panel in Figure 4. The COice and H2Oice abundance is a function of orbital time of the particle (in years). During the initial dust-settling phase, the COice and H2Oice abundance remains constant as the particle remains UV-shielded in the disk until it starts experiencing the vertical circling motion during the VSI non-linear regime at a few tens of hundreds of orbits. The big drop in both abundance profiles happens at ~50.9 Kyr (at 143 orbits at 46.5 au), when the particle reaches the UV-active layer (log10(FUV) = 9.11 photons s−1 cm−2) in the saturation-turbulence VSI phase. Here, the particle loses a significant amount of COice and H2Oice from its surface. These icy volatiles are expected to be completely absent from the grain surface when it is back to the disk midplane (crosshatched region).

In the text

	Fig. A.1 TW Hydra and AB Aur stellar spectra templates normalized at 1 pc used as an input parameter in RADMC3D.
In the text

	Fig. A.2 TW Hydra and AB Aur stellar spectra templates normalized at 1 pc used as an input parameter in RADMC3D.
In the text

	Fig. B.1 Radial temperature profile in the midplane comparison between DSHARP and DIANA opacities.
In the text

	Fig. C.1 Time evolution of the local kinetic energy (top) and local 1 mm dust distribution (middle) at 50 au. The red dashed line at Z = 5 HR represent the height when particles reach equilibrium with the gas. The bottom panel shows the middle panel overplotted on the top panel.
In the text

	Fig. D.1 Occurrence rate of particles above the Z/R = 0.2 layer. This occurrence rate represents the amount of times particles spent above Z/R = 0.2 layer. The x-axis shows time in orbits in which particles made their last UV cross.
In the text