The Effect of Luminosity Outbursts on the Abundance of Pebbles and Their Ice Mantles in Protoplanetary Disks

Hydrodynamic equations describing mass, momentum and internal energy conservation for the gas component are:

where $\Sigma_{\rm g}$ is gas surface density, $e$ is internal energy per unit surface area, $\cal P$ is vertically integrated gas pressure from the ideal gas equation ${\cal P}=(\gamma-1)e$ with $\gamma=7/5$, $f$ is gas-dust friction force, $\mbox{\boldmath$v$}=v_{r}\hat{r}+v_{\phi}\hat{\phi}$ is gas velocity in the disk plane, $\nabla=\hat{\mbox{\boldmath$r$}}\partial/\partial r+\hat{\mbox{\boldmath$\phi$}}r^{-1}\partial/\partial\phi$ is the planar gradient. Gravitational acceleration in the disk plane $\mbox{\boldmath$g$}=g_{r}\hat{\mbox{\boldmath$r$}}+g_{\phi}\hat{\mbox{\boldmath$\phi$}}$ includes central protostar gravity and gas/dust self-gravity calculated by solving the Poisson integral (Vorobyov & Basu, 2010). Expressions for radiative cooling $\Lambda$ and radiative heating $\Gamma$ (the latter including stellar and background radiation) can be found in Vorobyov & Elbakyan (2018). Background radiation temperature is set to 15 K. In FEOSAD, turbulent viscosity is included via the viscous stress tensor $\Pi$. Viscosity parameterization is described in more detail in Section 2.2.

Gas and dust co-evolve and interact via gravitational and friction forces, the latter including dust back reaction on gas following the analytical method by Stoyanovskaya et al. (2018). Initially all dust in the cloud is in submicron form but grows as the protoplanetary disk forms and evolves. Small submicron dust is assumed to be coupled to the gas, while grown dust can dynamically decouple from gas. This paper describes the dust growth scheme presented in Molyarova et al. (2021). In FEOSAD, small and grown dust populations are represented by their surface densities $\Sigma_{\rm d,sm}$ and $\Sigma_{\rm d,gr}$, respectively. Each dust population has a power-law size distribution $N(a)=Ca^{-p}$ with normalization constant $C$ and fixed exponent $p=3.5$. For small dust, minimum size $a_{\rm min}=5\times 10^{-7}$ cm, maximum size $a_{\ast}=10^{-4}$ cm. For grown dust, $a_{\ast}$ is the minimum size and $a_{\rm max}$ is the maximum size that can change via dust coagulation and fragmentation, as well as advective transport through the disk.

Dynamics of small and grown dust grains is described by the following continuity and momentum equations (small dust is strictly dynamically coupled to gas):

where $u$ is grown dust velocity. $S(a_{\rm max})$ represents the small-to-grown dust conversion rate per unit disk surface area as dust grows, i.e., as maximum dust size $a_{\rm max}$ increases (in g s^-1 cm^-2) (see details in Vorobyov et al., 2020, equation (14)). The diffusion coefficient is calculated from kinematic viscosity $D=\nu/Sc$, with Schmidt number $Sc=1$. Friction force $f$ calculation is described in Stoyanovskaya et al. (2020); Vorobyov et al. (2023).

In the dust growth model, maximum grown dust size $a_{\rm max}$ is calculated at each timestep for each cell as the solution to

where dust growth rate $\cal{D}$ describes $a_{\rm max}$ change due to coagulation, and the second left-hand term describes $a_{\rm max}$ change from dust advection between cells. We calculate $\cal{D}$ as

where $\rho_{\rm d}$ is dust volume density in the midplane, $\rho_{\rm s}=3$ g cm^-3 is refractory dust core density, and $v_{\rm rel}$ is dust grain collision velocity. The main collision sources are assumed to be Brownian motion and turbulence determined by the Shakura-Sunyaev $\alpha$ parameter. This approach is described in more detail in Vorobyov et al. (2018).

The $a_{\rm max}$ value is limited by the fragmentation barrier (Birnstiel et al., 2012); maximum dust size cannot exceed $a_{\rm frag}$ determined by dust properties. If $a_{\rm max}$ exceeds $a_{\rm frag}$ at any point, $\cal{D}$ is set to zero and $a_{\rm max}$ is set to $a_{\rm afrag}$. The maximum fragmentation-limited dust size is

where $v_{\rm frag}$ is the fragmentation threshold depending on ice mantle presence (see Section 2.4). Note the dust drift barrier is self-consistently included via coupled gas-dust dynamics with mutual friction.

The model implements a variable Shakura-Sunyaev $\alpha$-parameter used to calculate turbulent viscosity. Viscosity, significantly contributing to mass and angular momentum transport in protoplanetary disks, is primarily considered as caused by turbulence generated by magnetorotational instability (Balbus & Hawley, 1991; Turner et al., 2014). In FEOSAD, turbulent viscosity is included via the viscous stress tensor $\Pi$ (see equations  (2), (3)), with kinematic viscosity $\nu=\alpha c_{\rm s}H$, where $c_{\rm s}$ is sound speed and $H$ is the vertical scale height of the gas calculated assuming local hydrostatic equilibrium. To mimic accretion through a layered disk, we consider an effective adaptive $\alpha_{\rm eff}$ parameter (Eq. (12) from Kadam et al., 2022). The $\alpha_{\rm eff}$ parameter is calculated as a weighted average for combined active and dead magnetorotational instability layers:

The active layer sufficiently ionized by cosmic rays for magnetorotational instability development has surface density $\Sigma_{\rm MRI}=200$ g cm^-2 (100 g cm^-2 from each disk surface), while the dead zone has $\Sigma_{\rm dz}=\Sigma_{\rm g}-\Sigma_{\rm MRI}$. The viscosity parameter in the active layer is $\alpha_{\rm MRI}=10^{-3}$, while in the dead zone $\alpha_{\rm dz}=10^{-5}$ for temperatures below $1300$ K and $\alpha_{\rm dz}=10^{-1}$ otherwise, which corresponds to magnetorotational instability development from thermal ionization of alkali metals in inner disk regions. This $\alpha_{\rm eff}$ parameterization helps approximately account for turbulence suppression in the so-called dead zone and corresponding reduction in mass transport, as well as accretion and luminosity outbursts from episodic magnetorotational instability development in inner hot disk regions. More detailed descriptions can be found in Bae et al. (2014); Kadam et al. (2019, 2022).

Luminosity outbursts can be triggered by various mechanisms causing temporary increases in stellar accretion rates. A key mechanism is magnetorotational instability (MRI, Balbus & Hawley, 1991, 1998; Turner et al., 2014) development in inner disk regions. Corresponding $\alpha$ growth leads to accretion outbursts, or MRI outbursts (Zhu et al., 2010; Kadam et al., 2020; Vorobyov et al., 2020). This mechanism is self-consistently implemented in the presented model, leading to episodic growth of accretion luminosity. Each considered model experiences several dozen such outbursts during simulation, with amplitudes and durations similar to FUor outbursts (Audard et al., 2014; Connelley & Reipurth, 2018). These outbursts affect disk thermal structure, which in turn influences dust and pebble properties.

The Stokes number in protoplanetary disks typically describes dust dynamics. It is defined as:

where $\Omega_{\rm k}$ is Keplerian angular velocity, $\rho_{\rm g}=\Sigma_{\rm g}/\sqrt{\rm 2\pi}H$ is midplane gas volume density. Grains with $St\sim$1 are most strongly affected by gas and rapidly drift radially toward pressure maxima. Millimeter-to-centimeter sized dust grains called pebbles play a crucial role in planet formation scenarios involving pebble accretion (Ormel & Klahr, 2010; Lambrechts & Johansen, 2012; Ida et al., 2016; Lambrechts et al., 2017). Pebbles have relatively high $St>0.01$ and move relative to gas, unlike submicron grains that are coupled to gas and move together with it.

This work uses the pebble definition from Vorobyov et al. (2023), with modifications also described in Topchieva et al. (2024). FEOSAD considers two dust populations, small ($a<a_{\rm*}=1$ $\mu$m) and grown ($a>a_{\rm*}$) dust. We define pebbles as grown dust with Stokes number and grain radius exceeding threshold values $St>St_{\rm 0}$ and $a_{\rm peb,0}$. In disk regions where both conditions are met, pebbles have a size distribution that is part of the grown dust size distribution $N(a)=Ca^{-p}$, ranging from $a_{\rm peb,min}$ to $a_{\rm max}$. The minimum pebble size $a_{\rm peb,min}$ is:

Here $a_{\rm St_{\rm 0}}$ is the grain size with Stokes number $St_{\rm 0}$. In regions containing pebbles, it is defined as

We adopt $St_{\rm 0}=0.01$ and $a_{\rm peb,0}=0.05$ based on studies in Lambrechts & Johansen (2012); Lenz et al. (2019). The Stokes number identifies large grains moving independently from gas, but it also depends on local gas density. Our simulation includes an envelope surrounding the disk where $St$ reaches high values for micron-sized grains due to low gas density. The minimum pebble size condition $a_{\rm peb,0}$ excludes these small grains from consideration.

Pebble surface density $\Sigma_{\rm peb}$ in each computational cell is calculated assuming its size distribution continues the grown dust size distribution $N(a)=Ca^{-p}$, leading to:

Pebbles defined this way are part of grown dust within our two-population dust model constraints. Their dynamics is indistinguishable from overall grown dust dynamics described by equations (5) and (6). To distinguish these populations, we consider grown dust surface density minus pebble surface density unless stated otherwise. Ice surface densities on pebbles are calculated assuming ice mass fractions on pebbles equal those on grown dust.

Importantly, dust destruction at sudden $v_{\rm frag}$ drops is instantaneous in our modeling. This stems from the implementation of dust growth rate $\cal{D}$ and $v_{\rm frag}$ described in Section 2.1. When mantles evaporate, the fragmentation barrier $a_{\rm frag}$ changes abruptly, and all dust exceeding the new $a_{\rm frag}$ is recycled mainly into small dust plus some larger grains determined by the new $a_{\rm frag}$. This instantaneous destruction scenario is justified if grains are aggregates of small grains stuck together by ice mantles. Then mantle disappearance should instantly disrupt grains (see e.g. Schoonenberg et al., 2017). An alternative approach explicitly includes collisional destruction (erosion), as in Stammler & Birnstiel (2022). In that case destruction is not instantaneous but occurs on coagulation timescales of order $10^{3}$ years (Vorobyov et al., 2022), potentially longer than outburst durations. That scenario assumes ice mantles completely cover silicate grain cores. Both scenarios were considered in Houge et al. (2024), and the comparison to V883 Ori spectral indices favors prolonged dust destruction in that object. However, in self-gravitating protoplanetary disks with prominent spiral structure, complex snowline geometry should lead to mantle-stuck aggregates as shown in Molyarova et al. (2021). This motivates more detailed consideration of different grain disruption scenarios. Within FEOSAD this could be implemented via a bidisperse dust model including erosion as described in Akimkin et al. (2020), planned as future work extending this study.

The numerical simulation considers four volatile molecules: H₂O, CO₂, CH₄ and CO, which can exist in gas, on small dust and on grown dust. When we identify pebbles as part of grown dust in post-processing, we can separately consider ices on pebbles. Distribution of volatiles in the disk and between phases changes via three main processes. First, they advect with gas and dust. Second, during collisional dust evolution ice mass redistributes between dust fractions similarly to refractory core mass redistribution: if a certain mass fraction of small dust converts to grown dust, the same fraction of each ice type on small dust converts to ice on grown dust. We also consider phase transitions: adsorption onto dust grain surface, thermal desorption and photodesorption. The volatile evolution model is described in more detail in Molyarova et al. (2021).

The model assumes that initially, small dust grains are covered by ice mantles, grown dust has not yet formed, and gas-phase volatiles are absent. Initial relative molecular abundances are based on ice observations in low-mass protostellar cores Öberg et al. (2011) with ratios H₂O : CO₂ : CO : CH₄ $=100:29:29:5$. Initial ice mass fraction relative to refractory material is $\approx$8.5%. Protoplanetary disks can have different abundances of volatiles. Since volatile mass is not included in dust dynamics, initial composition choice will not qualitatively affect the results. It may influence ice ratios but not snowline positions or dust properties.

Adopted initial abundances characterize prestellar core compositions representing the stage immediately preceding protostar and protoplanetary disk formation, they are often used for astrochemical modeling (Eistrup et al., 2016, 2018). Carbon elemental abundance in this chemical composition is reduced relative to solar, implying the presence of solid carbonaceous dust and refractory organics in the disk not included in current modeling. Accounting for these components and their possible transition to gas phase requires considering chemical processes beyond desorption and adsorption, which is beyond the scope of this work, but could potentially affect gas composition in the inner disk (Lee et al., 2010; Wei et al., 2019).

To estimate the effect of ice mantles on dust evolution we use fragmentation threshold $v_{\rm frag}$ as a parameter depending on ice mantle presence on grown dust. At collision velocities below $v_{\rm frag}$ grown grains stick together, above it they fragment. We assume that ice mantles generally increase $v_{\rm frag}$ regardless of composition (Molyarova et al., 2021, see), based on laboratory studies Wada et al. (2009); Gundlach & Blum (2015). Ice-covered grains have $v_{\rm frag}=5$ m s^-1, bare grains have $v_{\rm frag}=0.5$ m s^-1. Grown dust in a disk region is considered ice-covered if there is enough ice on grown dust to cover all present grown grains by at least one ice monolayer (thickness order of molecular size). Monolayer thickness is estimated using characteristic water molecule size $3\times 10^{-8}$ cm. After calculating ice surface densities at each timestep, $v_{\rm frag}$ values are updated.

We consider two models with different initial cloud masses $\rm M_{\rm core}=$ 0.66 $M_{\odot}$ and 1.0 $M_{\odot}$. Model M1 serves as baseline, while M2 investigates how model mass affects pebble property conclusions. Disk-to-star mass ratios are relatively high and similar (0.55 and 0.60), making both disks gravitationally unstable, but total system mass also affects gas and dust evolution. Central star mass in the model is self-consistently determined by accretion flow from the disk. Over 0.5 Myr the disk retains significant mass, maintaining gravitational instability and thus inward mass transport needed to activate magnetorotational instability (Bae et al., 2014).

At 0.5 Myr the disk remains massive, which is reproduced in other numerical protoplanetary disk formation models (e.g. Deng et al., 2020). Observational determination of disk masses mainly relies on indirect methods using CO and grown dust emission, and potentially underestimates masses by 1-2 orders of magnitude due to optical depth effects and model uncertainties (Dunham et al., 2014; Miotello et al., 2014; Molyarova et al., 2017). More direct methods like gas kinematics yield values comparable to our simulations (e.g. Lodato et al., 2023). Main model parameters are listed in Table LABEL:tab:model. The simulations used a $N_{r}\times N_{\phi}=390\times 256$ grid. Disk formation occurs at 53 kyr in model M1 and 79 kyr in model M2.

Luminosity outbursts accompanying accretion rate increases heat the disk causing thermal ice desorption, shifting all volatile snowlines outward (Cieza et al., 2016a; Vorobyov et al., 2022; Houge & Krijt, 2023). This inevitably affects ice mantle composition on dust grains, including pebbles. Fig. 2 shows temporal evolution of total mass of various ices on pebbles. At certain times ice masses drop sharply, corresponding to luminosity outbursts in Fig. 1.

Some outbursts affect all ices, while others only reduce CH₄ and CO masses. These molecules are more volatile with lower desorption energies, their snowlines are farther from the star and do not always coincide with pebble regions (Topchieva et al., 2024). Their masses on pebbles are orders of magnitude lower than those of H₂O and CO₂. Thus even relatively weak luminosity outbursts can strongly affect their abundance in the disk if the snowlines lie in regions where stellar irradiation dominates heating rather than viscous heating.

Fig. 3 shows the snowlines of the four molecules versus time. At certain times snowlines shift outward abruptly, corresponding to luminosity outbursts. For most weak outbursts the water snowline barely moves, because it lies in disk regions where turbulent viscous heating dominates over stellar irradiation. In contrast, CH₄, CO₂ and CO snowlines located at radii where stellar irradiation dominates heating show more pronounced shifts. This highlights the importance of analyzing radial boundaries where stellar irradiation heating dominates other heat sources for understanding the behavior of these molecules during outbursts Vorobyov et al. (2018).

At later times multiple water snowlines form in the disk. At $\approx 1-2$ au a region appears where water exists in ice phase. In model M1 this occurs at $\approx 250$ kyr, in model M2 it happens after $\approx 450$ kyr. The formation of water ice in the inner region is associated with a $\approx 1-2$ au dust ring accumulating grown dust (see Topchieva et al., 2024). This ring forms in the dead zone with the lowest $\alpha$-parameter and high gas/dust surface density. Consequently, heating is weakened there, causing relative temperature drops sufficient for water freeze-out.

Although the water snowline does not shift much and water ice mass on pebbles varies less than other molecules, water remains one of the most important molecules in protoplanetary disks. All disk pebbles are covered with water ice (Topchieva et al., 2024). This, combined with fragmentation velocity dependence on ice mantle presence (Section 2.4), makes water the key pebble surface molecule deserving detailed consideration.

Let us examine outburst effects on pebbles and their composition in more detail using characteristic events in both models. In M1 this is the 244.5 kyr event, in M2 the $\sim 340$ kyr event.

The first considered luminosity outburst profile and corresponding masses of refractory components are shown in right panels of Fig. 4. Luminosity peaks at 244.5 kyr. The outburst lasts about two hundred years with $\approx 70$ $L_{\odot}$ magnitude. Pebble mass begins decreasing simultaneously with luminosity rise; by 245 kyr, i.e. after the outburst, it drops by roughly half, then recovers to pre-outburst values in $\sim 4$ kyr. As pebble mass falls, small dust mass rises while grown dust mass changes only slightly. During the outburst pebbles are converted to small and grown dust. Grown dust also fragments to small dust in regions with no pebbles. Amount of grown dust increases from pebble destruction but it also fragments, so total dust mass changes weakly. In several thousand years around the outburst total dust mass shows slight decreases due to radial drift and accretion to the star.

This outburst is typical for early evolution in the M1 model. Similar changes in pebble mass occur during other outbursts, as seen in Fig. 4 lower left panel. During such outbursts mass of pebbles in the disk decreases by a factor of $\approx 2$.

The luminosity outburst heats the disk, but temperature increases vary between disk regions depending on local heating mechanism contributions to thermal balance. This is illustrated in Fig. 4 bottom panel showing radial temperature profiles before, during and after the outburst. Pre- and post-outburst profiles are similar but both differ strongly from the outburst profile. Beyond $\approx 7-8$ au midplane temperature $T_{\rm m.p.}$ approaches radiative temperature $T_{\rm irr}$, which is the temperature calculated assuming radiation is the sole heating source. In these regions external irradiation dominates heating, and luminosity outbursts have greatest impact, roughly doubling temperatures. Inside 7 au temperature exceeds radiative temperature by a factor of a few. Here thermal structure is determined by viscous heating. Within $\approx 4$ au MRI develops, increasing $\alpha$, enhancing viscous heating and raising temperatures. Between $4-7$ au $\alpha$ does not change during the outburst and radiative heating contribution is small, so temperatures stays approximately the same. The outburst significantly heats all disk regions except $4-7$ au. This effect occurs for outbursts with amplitudes around a hundred solar luminosities.

Pebble mass reduction directly results from luminosity outbursts. Additional accretion luminosity heating causes thermal desorption of ice mantles across large disk areas. This is evident from significant snowline shifts in Fig. 3. Near the water snowline on the inner side (closer to star), grains lose ice mantles to evaporation. Water desorption also occurs in spiral structures where temperatures are significantly elevated due to local gas heating. The zone where water does not desorb from dust includes not only regions outside the snowline but also spiral structures permeating the disk, which we examine in more detail below.

Grains in these regions are heated up and destroyed into monomers. This is also seen from corresponding increases in both grown and small dust masses from pebble fragmentation (we remind the reader that grown dust in our model is everything above one micron). When outbursts end, volatiles return to grain surfaces and dust growth resumes. Post-outburst volatile freeze-out times can be estimated as inverse adsorption rates onto dust (see equation (27) in Topchieva et al., 2024). Near water and CO₂ snowlines freeze-out times are at most days. In outer regions with CO and CH₄ snowlines, lower temperatures and dust surface densities can lead to freeze-out times up to decades.

Characteristic times for grains to return to pre-outburst states and become pebbles again are around thousands of years. They are determined by combined timescales of volatile freeze-out and dust growth. Water freeze-out times in these relatively high-density regions are very short compared to outburst durations, not exceeding several years as shown by astrochemical modeling of outburst chemical impacts Wiebe et al. (2019). Vorobyov et al. (2013) showed that CO and CO₂ freeze-out times in surrounding envelopes are thousands of years due to lower densities and temperatures favoring slower freeze-out. Here post-outburst pebble formation times around a thousand years are set by dust growth rates. Vorobyov et al. (2022) showed that post-outburst dust growth times are shorter closer to the star, with typical values $\sim 10^{3}$ years.

Fig. 5 shows radial distributions of solid disk components before, during and after the considered outburst. The comparison between the left and the middle panel demonstrates thermal desorption effects from additional stellar accretion luminosity heating. During the 245 kyr outburst CO₂ fraction drops sharply over a wide radial range, while CH₄ and CO disappear completely. Their snowlines lie beyond $>7$ au where accretion luminosity heating is important. Meanwhile the water snowline position does not change at all, remaining at 6 au. It lies at radii where outburst heating is weak.

As Fig. 5 shows, water ice radial distribution changes little with the outburst. Yet total disk water ice mass clearly decreases during outbursts, as seen in Fig. 2. This apparent contradiction results from non-axisymmetric structure of the disk. Fig. 6 shows 2D distributions of H₂O ice fraction on pebbles before and during the outburst. Water does sublimate between $\approx 7-11$ au but not at all azimuths, leaving a spiral-shaped region where pebbles remain icy. The H₂O snowline lies at radii where viscosity and radiation compete for thermal structure dominance, leading to complex snowline geometry evolution. In regions where grains lose ice mantles, they immediately fall apart into monomers and their sizes decrease. These regions contain about half of the total mass of pebbles in the disk. Pebbles are destroyed during the outburst as their ice mantles evaporate. This disruption is not visible in radially averaged distributions (Fig. 5) because regions with sublimated ice are averaged together with the regions having high H₂O ice fractions and pebble mass. Ice sublimation and pebble destruction during the outburst strongly depend on disk 2D structure.

Luminosity outbursts like those above also occur in model M2. Like the M1 outburst, these are accompanied by pebble mass reductions caused by disruption of grains that lost their ice mantles into monomers and subsequent pebble recovery as water ice refreezes. Radial temperature profiles before, during and after the 339 kyr outburst in Fig. 7 show temperature changes similar to Fig. 4 bottom panel.