© The Authors 2026

1 Introduction

Observations of young stellar objects have revealed that many systems exhibit strong variability in their luminosity, often attributed to episodic accretion events (Herbig 1977; Audard et al. 2014). Of particular interest are the so-called FU Orionis-type outbursts (e.g., Hartmann & Kenyon 1996), characterized by sudden increases in brightness by several magnitudes that can last for decades to centuries (for a review, see Fischer et al. 2023).

Theoretical models have proposed various mechanisms to explain these outbursts, with a prominent candidate being a thermal instability (TI) operating at the interface between the inner edge of the dead zone and the ionized active interior (Dullemond & Monnier 2010). The latter region is thought to be sufficiently hot and ionized to sustain turbulence via the magnetorotational instability (MRI; Balbus & Hawley 1991; Hawley et al. 1995), with the dead zone being colder and significantly less turbulent due to insufficient ionization quenching the MRI (Gammie 1996; Bai & Stone 2013b), albeit with some weak level of hydrodynamic turbulence (for a review, see Lesur et al. 2023).

This contrast in turbulent angular momentum transport at the dead zone inner edge (DZIE) leads to steady mass accumulation there (Flock et al. 2017), until the TI is triggered when turbulent dissipation (i.e., viscous heating) exceeds radiative cooling even for the trace amounts of turbulence in the dead zone. In turn, this results in a runaway heating episode, which subsequently ionizes the gas, triggering the MRI and inducing a burst of viscous accretion onto the central star (Armitage et al. 2001; Zhu et al. 2009a; Cecil & Flock 2024). An example of a burst cycle via this process is shown in Fig. 1 using an “S-curve” representation of the thermal equilibrium states of the disk (e.g., Faulkner et al. 1983).

The evolution of an accretion burst via this variant of the TI has been the topic of several numerical studies. Early works by Armitage et al. (2001) and Zhu et al. (2009a) used onedimensional (1D) viscous disk models with simplified thermal physics to demonstrate the viability of this mechanism, while Chambers (2024) used more sophisticated 1D cooling prescriptions that recovered similar results. More recent studies have employed twodimensional (2D) radiation hydrodynamics simulations to capture the complex interplay between heating, cooling, and angular momentum transport at the DZIE (Wünsch et al. 2006; Cecil & Flock 2024), revealing the presence of reflares within a burst cycle and the formation of multiple rings of high gas density within the dead zone (∼1 au) as a byproduct of the outburst. Such features were also found by Chambers (2024).

This substructure is particularly interesting from the perspective of planet formation, with the pressure maxima associated with the rings acting as traps for inward-drifting dust particles (Weidenschilling 1977a; Takeuchi & Lin 2002). The accumulation of dust in these regions could facilitate the formation of planetesimals via the streaming instability (Youdin & Goodman 2005; Johansen & Youdin 2007; Baronett et al. 2024; Aly & Paardekooper 2025) and their subsequent growth via either pebble accretion (e.g., Lambrechts & Johansen 2012, 2014) or planetesimal accretion (e.g., Pollack et al. 1996; Weidenschilling 2000). If such a formation pathway is viable, it could explain the presence of thousands of super-Earths and mini-Neptunes found within ∼1 au of their host stars¹ (Fressin et al. 2013; Petigura et al. 2013, 2018) via essentially in situ formation.

With the dust playing a central role in the thermodynamics of accretion outbursts by being the primary carriers of opacity in the dead zone, it is crucial to treat dust dynamics self-consistently when modeling this process. However, previous works on the subject have assumed a fixed dust size distribution of grains perfectly coupled to the gas, with the dust-to-gas ratio being an input parameter rather than a dynamically evolving quantity through radial drift, accumulation at pressure bumps, and dust growth.

In this work, we aim to address this gap by modeling dust as a fully dynamic component that evolves via coagulation and fragmentation, and interacts with the gas both dynamically via drag forces but also thermally via size- and temperature-dependent dust opacities. This will allow us to assess the impact of dust evolution during and after an accretion outburst, as well as the potential for planet formation in the resulting dust rings.

Below, we will refer to the quiescent state before an outburst as the “pre-burst” phase, the outburst itself as a “burst cycle” (or simply “burst”), and the disk state after the burst has subsided as the “post-burst” phase. During the long quiescent phase, this post-burst state evolves viscously until the next pre-burst state is reached, where the conditions for a burst are once again satisfied.

The burst cycle itself can consist of multiple “flares,” with the first being the most intense one, and “reflares” being subsequent, less intense events within the same burst cycle. We will also use the terms “dead zone inner edge” (DZIE) and “inner rim” interchangeably.

In Sect. 2, we describe the theoretical framework used in our work. We then make semi-analytical predictions on the state of the disk before and after a burst in Sect. 3. The results of our radiation hydrodynamical simulations are presented in Sect. 4, followed by simplified long-term models of dust evolution that examine the prospects for planetesimal formation in Sect. 5. We discuss our results in Sect. 6, and conclude in Sect. 7.

Fig. 1 “S-curve” of the gas surface density versus temperature at R = 0.3 au for our fiducial model, discussed in detail in Sect. 4.1. 1: runaway heating is triggered at the DZIE. 2: the burst front passes through. 3: dust sublimation caps temperature during burst. 4: dust recondensation; 5: dust cooling. Gray bands highlight the stable branches corresponding to viscous evolution in the quiescent (bottom) and burst (top) phases. Different colors denote reflares during the same burst cycle.

2 Physics and numerics

In this section, we describe the physical and numerical framework used in our study. We begin by listing the relevant hydrodynamical equations and describing the different processes involved. We then describe the numerical setup used in our hydrodynamical simulations in Sect. 4.

2.1 Multifluid radiation hydrodynamics

We assumed a disk of gas with surface density Σ_g, velocity field u_g, internal energy density e, mean molecular weight µ = 2.353, and adiabatic index γ = 7/5, as well as a pressureless dust component with surface density Σ_d and velocity u_d, orbiting a star of mass M_⋆ = 1 M_⊙ and luminosity L_⋆ = 1.78 L_⊙ (Hayashi 1981). The vertically integrated Navier–Stokes equations then read $$\frac{{\text{dΣ}}_{\text{g}}}{\text{d}t}=-{\text{Σg}}_{}\nabla \cdot u_{\text{g}},$$(1a) $$\frac{\text{d}{u}_{\text{g}}}{\text{d}t}=-\frac{1}{{\text{Σg}}_{}}\nabla P-\nabla {\text{Φ}}_{\star }+\nabla \cdot \overline{\sigma }-{\displaystyle \sum }_i{ϵ}_i\frac{u_{\text{d},i}-u_{\text{g}}}{{\text{St}}_i}{\text{ΩK}}_{},$$(1b) $$\frac{\text{d}e}{\text{d}t}=-\gamma e\nabla \cdot u_{\text{g}}+Q_{\text{visc}}+Q_{\text{irr}}+Q_{\text{cool}}+Q_{\text{rad}},$$(1c) $$\frac{{\text{dΣ}}_{\text{d},i}}{\text{d}t}=-{\text{Σ}}_{\text{d},i}\nabla \cdot u_{\text{d},i}-\nabla \cdot F_{\text{diff},i},F_{\text{diff},i}=-{\text{Σg}}_{}\frac{\nu }{1+{\text{St}}_i^2}\nabla ϵ$$(1d) $$\frac{\text{d}{u}_{\text{d},i}}{\text{d}t}=-\nabla {\text{Φ}}_{\star }-\frac{u_{\text{d},i}-u_{\text{g}}}{{\text{St}}_i}{\text{ΩK}}_{}.$$(1e)

In the above equations, d/dt refers to the material derivative, P = (γ − 1)e is the gas pressure given by the perfect gas law, Φ_⋆ = −GM_⋆/R is the gravitational potential of the star at distance R, $\overline{\sigma }$ is the viscous stress tensor, ${\Omega }_{\text{K}}=\sqrt{G{M}_{*}/R^3}$ is the Keplerian angular frequency, and ν is the kinematic viscosity of the gas. The gravitational constant is denoted by G. Through the above we can define the isothermal sound speed $c_{\text{s}}=\sqrt{P/{\text{Σg}}_{}}$, the temperature $T=\mu c_s^2/R$, the pressure scale height H = c_s/Ω_K, and the aspect ratio h = H/R, with ℛ being the ideal gas constant.

Regarding the dust component for a species with index i, we defined the dust-to-gas ratio ϵ_i = Σ_d,i/Σ_g and the Stokes number (i.e., dimensionless stopping time) of a species with radius a_i and material density $\tilde{\rho }=2.08{\text{\hspace{0.17em}g/cm}}^3$ in the Epstein regime as $${\text{St}}_i=\frac{\pi }{2}\frac{a_i\tilde{\rho }}{{\text{Σg}}_{}}.$$(2)

The term F_diff,i then represents dust diffusion due to turbulence (Morfill & Voelk 1984; Youdin & Lithwick 2007).

The terms Q_visc, Q_irr, Q_surf, and Q_rad in Eq. (1c) represent viscous heating, stellar irradiation heating, radiative cooling through surface losses, and radiative diffusion through the disk midplane, respectively. Similar to Ziampras et al. (2025a), we defined these terms as follows: $$Q_{\text{visc}}=\frac{1}{2\nu {\text{Σg}}_{}}\text{Tr}({\overline{\sigma }}^2)\approx \frac{9}{4}\nu {\text{Σg}}_{}{\text{ΩK}}_2^{},\nu =\alpha \sqrt{\gamma }{c}_{\text{s}}H,$$(3a) $$Q_{\text{irr}}=2\frac{L_{\star }}{4\pi R^2}(1-\epsilon )\frac{\theta }{{\tau }_{\text{eff}}},$$(3b) $$Q_{\text{surf}}=-2\frac{{\sigma }_{\text{SB}}{T}^4}{{\tau }_{\text{eff}}},{\tau }_{\text{eff}}=\frac{3{\tau }_{\text{R}}}{8}+\frac{\sqrt{3}}{4}+\frac{1}{4{\tau }_{\text{P}}},{\tau }_{\text{R,P}}=\frac{{\kappa }_{\text{R,P}}{\text{Σg}}_{}}{2},$$(3c) $$Q_{\text{rad}}=\sqrt{2\pi }H\nabla \cdot \left(\lambda \frac{4{\sigma }_{\text{SB}}}{{\kappa }_{\text{R}}{\rho }_{\text{mid}}}\nabla T^4\right),{\rho }_{\text{mid}}=\frac{1}{\sqrt{2\pi }}\frac{{\text{Σg}}_{}}{H}.$$(3d)

In the above, we utilize the α-viscosity model of Shakura & Sunyaev (1973), an irradiation model based on Menou & Goodman (2004) with a disk albedo ε = 1/2 and flaring angle θ, an effective optical depth τ_eff following Hubeny (1990), and in-plane radiative transfer in the one-temperature flux-limited diffusion approximation (Levermore & Pomraning 1981; Rometsch et al. 2024). The Rosseland and Planck mean opacities are denoted by κ_R and κ_P, respectively; ρ_mid is the midplane gas volume density; λ is the flux limiter from Kley (1989), and σ_SB is the Stefan– Boltzmann constant. Our prescriptions specific to the dead zone inner edge for α, θ, κ_R, and κ_P are described in Sects. 2.2 and 2.3.

We note here that in-plane radiative diffusion via Q_rad is only included in a subset of our simulations due to its steep computational cost during the post-burst phase. Its effects and relevance are explored in detail in Sect. 4.2.

2.2 Dust evolution and opacity model

In addition to the set of equations in Eq. (1), we also solve for the evolution of the dust size distribution n(a) using the TriPoD method (Pfeil et al. 2024, hereafter P+24). This approach simplifies the full coagulation equations (Smoluchowski 1916) by instead tracking the evolution of representative “small” and “large” grain populations with surface densities Σ_small and Σ_large via coagulation and/or fragmentation, along with the maximum grain size a_max (following P+24, we fix a_min = 0.1 µm). This method can be seen as a more accurate extension of two-pop-py (Birnstiel et al. 2017), with a wider range of applicability. For the full documentation, implementation, and testing of the TriPoD method we refer the reader to P+24.

In the TriPoD approach, Σ_large and Σ_small define a power-law dust size distribution, truncated at a maximum grain size a_max. Given the two populations are divided at an intermediate grain size $a_{\mathrm{int}}=\sqrt{a_{\text{min}}{a}_{\text{max}}}$, they define the power-law exponent of the distribution q_dust via $$q_{\text{dust}}=\frac{\log ({\text{Σ}}_{\text{large}}/{\text{Σ}}_{\text{small}})}{\log (a_{\text{max}}/a_{\text{int}})}-4,$$(4)

such that ${\text{Σ}}_{\text{d},i}\left(a\right)\propto a^{q\text{dust+4}}$. The cutoff size a_max is treated as an additional “fluid” that is advected with the large population’s velocity field, that is, tracing the advection of the largest grains. Figure 2 illustrates the reconstructed dust size distribution using Σ_small, Σ_large, and a_max.

Since dust is modeled as a set of fluids with their own velocities in our PLUTO implementation, we further account for momentum conservation during coagulation and fragmentation by exchanging momenta between Σ_small and Σ_large during the TriPoD step. This is not the case in the original TriPoD implementation of P+24. Simple tests showed that this does not affect the results in any meaningful way, so a recalibration of the method was not necessary, but we include this detail for completeness.

Regarding our opacity model, we can write the total Rosseland or Planck opacity κ in Eqs. (3c) and (3d) as the sum of dust and gas opacities κ_d and κ_g, expressed in cm² per gram of gas: $$\kappa =f_{\text{subl}}ϵ{\kappa }_{\text{d}}+{\kappa }_{\text{g}},ϵ=\frac{{\text{Σ}}_{\text{small}}+{\text{Σ}}_{\text{large}}}{{\text{Σg}}_{}},$$(5)

with f_subl as a prefactor that accounts for dust sublimation at high temperatures (see Eq. (7) in Sect. 2.3). For our dust opacities, we used the a_max-, q_dust-, and T-dependent tables computed with the growpacity² module (Ziampras & Birnstiel 2026) and fixed κ_R,gas = κ_P,gas = 10⁻³ cm²/g_gas following Cecil & Flock (2024) (hereafter CF24). A summary of the growpacity method is provided in Appendix C.

Fig. 2 Sample dust size distribution, represented with a truncated power-law with qdust = −3.5 and reconstructed with the TriPoD method of Pfeil et al. (2024).

2.3 Model adaptation to the dead zone inner edge

Interior to the dead zone inner edge, the MRI is thought to operate at full strength, driving intense levels of turbulence (Flock et al. 2017; Iwasaki et al. 2024). Following CF24, we modeled the MRI activity as a temperature-dependent α, under the assumption that the gas is sufficiently ionized above a temperature T_MRI such that the MRI can operate: $$\alpha (T)={\alpha }_{\text{DZ}}+\frac{1}{2}({\alpha }_{\text{MRI}}-{\alpha }_{\text{DZ}})\left[1+\tanh \left(\frac{T-T_{\text{MRI}}}{\text{Δ}{T}_{\text{MRI}}}\right)\right],$$(6)

with α_MRI = 0.1, T_MRI = 900 K, and ∆T_MRI = 25 K. Similarly, we model the sublimation of dust grains beyond a critical temperature $T_{\text{subl}}=2000\text{K}{\rho }_{{\text{mid,g/cm}}^3}^{0.0195}$ (Isella & Natta 2005) with a sublimation fraction given similar to CF24: $$f_{\text{subl}}={10}^{-10}+\frac{1}{2}(1-{10}^{-10})\left[1-\tanh \left(\frac{T-T_{\text{subl}}}{\text{Δ}{T}_{\text{subl}}}\right)\right],$$(7)

with ∆T_subl = 100 K.

We note that, in our model, the sublimation fraction above does not directly affect the dust densities Σ_d,i, and sublimated dust vapor is allowed to instantly recondense as temperatures drop below T_subl. This simplification does not harm the accuracy of our model, however, as the temperatures needed to sublimate the dust practically always satisfy T_subl > T_MRI, translating to extreme values of α that lead to rapid fragmentation of grains. As a consequence to that, the typical grain sizes during this state are very low (∼µm), and the “dust vapor” has a Stokes number so small that it is practically coupled to the background gas flow, effectively evolving as a “gas” fluid. As the recondensed dust grows back to larger grains via coagulation, the dust population smoothly decouples from the gas once again.

Finally, to model the direct irradiation of the dust-free, hot inner disk interior to the dead zone edge, we modify the grazing angle in Eq. (3b) from its usual profile of θ₀ ≈ 2χh/7 with χ = 4 (Chiang & Goldreich 1997) to reach a value of 1 (i.e., direct illumination) in the dust-free region below R_rim = 0.1 au. This value for the dead zone inner edge was chosen matching the results of CF24, and the transition was modeled similarly to Eqs. (6) and (7) above as $$\theta (R)={\theta }_0+\frac{1}{2}(1-{\theta }_0)\left[1-\tanh \left(\frac{R-R_{\text{rim}}}{0.01\text{au}}\right)\right].$$(8)

This results in a significantly hotter region interior to R_rim that always satisfies α = α_MRI through Eq. (6) and ensures that the disk has a cavity within R_rim (see, e.g., Fig. 4). We note that more complex and realistic models of the temperature profile near the inner rim exist (e.g., Ueda et al. 2017), taking into account the full vertical structure of the disk, but this approach is sufficient for the purposes of this work as we are focused on the disk midplane.

2.4 Numerics

For our hydrodynamical simulations in Sect. 4 we use the finite-volume PLUTO v.4.4 code (Mignone et al. 2007) with the HLLC Riemann solver (Toro et al. 1994), the flux limiter by Van Leer (1974), and second-order spatial and temporal accuracy. The modules for FLD, dust–gas interaction, and dust growth via TriPoD are documented in Ziampras et al. (2020), Ziampras et al. (2025c), and P+24, respectively.

We use an axisymmetric cylindrical polar grid that extends radially between 0.05–10 au, with 1024 cells between 0.05–2 au and 128 beyond 2 au, for a total of 1152 cells. Both grid segments stretch outwards with logarithmic spacing. We use a strict outflow boundary condition at the inner boundary for all quantities, and apply a damping zone in the region beyond 8.5 au towards the initial conditions using the formula in de Val-Borro et al. (2006) and a damping timescale of 0.1 local orbits.

We typically initialize the disk with a surface density profile that is either nearly or already unstable to an accretion outburst, and a power-law temperature profile that corresponds to a passively irradiated disk with $T_0\approx 110R_{\text{au}}^{-3/7}$ that nevertheless quickly converges to the equilibrium set by Q_visc + Q_irr = Q_cool in regions where viscous heating or direct irradiation matter. To facilitate a comparison among our models with different α_DZ, we impose that the surface density in the outer disk is always ${\text{Σ}}_{\text{disk,0}}=200R_{\text{au}}^{-15/14}{\text{g/cm}}^2$. The gas velocity is then initialized near viscous and hydrodynamical equilibrium with u_R = −1.5ν/R and $u_{\varphi }=R{\Omega }_{\text{K}}\sqrt{1-3h^2\text{/}2.}$.

For the dust, we initialize both dust densities to 0.005 Σ₀ for a total dust-to-gas ratio ϵ = 0.01, the maximum grain size at a uniform a_max = 1 mm, and the dust velocities to u_R,i = 0 and u_ϕ,i = RΩ_K. As with the gas, these choices are largely inconsequential as dust drift and evolution will quickly adjust the dust size distribution to its equilibrium profile before a burst cycle, but we list them for completeness.

Similar to CF24, we treat our model with α_DZ = 10⁻³ as our fiducial setup and run additional models with α_DZ ∈ [3 × 10⁻⁴, 10⁻⁴] to investigate the role of diffusion and dust evolution during and shortly after the burst phase. To highlight the impact of dynamic dust evolution with TriPoD we additionally run a model where we disable dust growth/fragmentation but instead prescribe a radial profile of a_max and q_dust informed by our fiducial model. Finally, we note that due to very high computational costs during the quiescent phase we opted to not include in-plane radiative diffusion (Q_rad, Eq. (3d)) in most of our models, but we nevertheless carry out a run where Q_rad is included during the burst phase to showcase its effect and relevance during and after the burst.

For our results with α_DZ = 10⁻³ in Sect. 4.1, we simply initialize the gas surface density with Σ₀ = Σ_disk,0. This will result in an artificially more massive burst, but viscous timescales are short enough that we can discard this first burst cycle and evolve until the next burst self-consistently triggers, which we then analyze. For our models with α_DZ ≤ 3 × 10⁻⁴ these viscous timescales are prohibitively long, and we therefore initialize the gas surface density with the pre-burst state predicted using the method described in Sect. 3 and tested in Appendix A.

3 Constraints on the pre- and post-burst states

Here, we will use the physical framework described in Sect. 2.1 to derive approximate solutions to the pre- and post-burst states near the inner edge of the dead zone. We will also provide illustrative examples of these solutions.

3.1 Pre-burst state

As shown by the aforementioned studies (e.g., Armitage et al. 2001; Zhu et al. 2009a; Cecil & Flock 2024), an accretion burst is triggered once sufficient material has accumulated at the dead zone inner edge such that the viscous heating term Q_visc ∝ Σ_g can overcome the radiative cooling term $Q_{\text{cool}}\propto {\text{Σg}}_{-1}^{}$ (for τ ≫ 1, which is easily satisfied in the optically thick inner disk). This triggers a runaway heating phase where Q_visc ∝ α(T) further increases, until the temperature reaches T ≈ T_subl and is regulated by a “thermostat effect”: a drop in temperature would allow dust to recondense and thus the opacity to increase, promoting heating, whereas a further increase to temperature is not feasible due to the lack of a dust opacity and a very low gas opacity (see, however Cecil et al. 2026, for the role of realistic gas opacities in the disk atmosphere). During this “hot” state, the gas is evolving at viscous rates with α ≈ α_MRI, with a substantial fraction of the gas mass within 1 au accreted onto the central star by the end of this phase. The remaining gas is rarefied and optically thinner, allowing radiative cooling to once again overpower viscous heating and lowering temperatures back to the quiescent, “cold” state.

With this picture in mind, we can already place constraints on the shape of the pre- and post-burst states in terms of the gas density and temperature. In particular, during the burst phase, the “burst front” traveling outwards cannot exceed the radius beyond which Q_visc ≤ Q_cool even for α ∼ α_MRI. Assuming that the “hot” state is characterized by Q_visc ≫ Q_irr and T ≈ T_MRI, Eq. (6) yields α(T_MRI) = α_MRI/2, and we have $$\begin{array}{lll}{Q}_{\text{visc}}\left(T_{\text{MRI}}\right)\hfill & =\hfill & Q_{\text{cool}}\left(T_{\text{MRI}}\right)\hfill \\ \hfill & \Rightarrow \hfill & \frac{9}{4}\frac{{\alpha }_{\text{MRI}}}{2}\sqrt{\gamma }{c}_{\text{s}}^2{\text{Σ}}_{\text{g,min}}{\Omega }_{\text{K}}=2\frac{{\sigma }_{\text{SB}}{T}_{M\Re I}^4}{\frac{3}{8}\frac{1}{2}{\kappa }_{\text{R}}\left(T_{\text{MRI}}\right){\text{Σ}}_{\text{g,min}}}\hfill \\ \hfill & \Rightarrow \hfill & {\text{Σ}}_{\text{g,min}}=\sqrt{\frac{1}{{\alpha }_{\text{MRI}}{\Omega }_{\text{K}}}\frac{256}{27}\frac{\mu {\sigma }_{\text{SB}}}{R\sqrt{\gamma }}\frac{T_{\text{MRI}}^3}{{\kappa }_{\text{R}}\left(T_{\text{MRI}}\right)}}.\hfill \end{array}$$(9)

In other words, the density must be at least Σ_g,min in order to sustain the outward-traveling heating front, and the burst flare can no longer propagate outwards when this condition is no longer satisfied. As a result, excluding the always-MRI-active region interior to R_rim, the disk surface density profile in the immediate post-burst state can be approximated by $${\text{Σ}}_{\text{g,post}}\approx \min \left({\text{Σ}}_{\text{g,min}},{\text{Σ}}_{\text{g,disk}}\right),$$(10)

where Σ_g,disk is the unaffected density profile of the disk outside of the burst region. We note that, for constant α_MRI, T_MRI, and a temperature-dependent opacity model, Eq. (9) simplifies to $${\text{Σ}}_{\text{g,min}}\propto {\text{ΩK}}_{-1/2}^{}\propto R^{3/4},$$(11)

which provides a very good match to the 3D axisymmetric models of CF24 (Σ_g,min ∝ R^0.81 therein). As we show in Fig. 3, this profile of Σ_g,min ∝ R^0.81 is exactly reproduced in our models due to the dependence of κ on T through f_subl in Eq. (7).

The corresponding temperature profile in the post-burst state can then be computed by solving Eq. (1c) in thermal equilibrium, or $$Q_{\text{visc}}(T^{\prime })+Q_{\text{irr}}(T^{\prime })=Q_{\text{cool}}(T^{\prime })$$(12)

with respect to T′, using Σ_g,post from Eq. (10). This can be done numerically via a root-finding algorithm such as Newton– Raphson or bisection.

Equation (10) offers some insight on the radial extent of the disk that can be “ignited” during an accretion outburst via this mechanism, with more massive disks being prone to wider regions where the surface density profile can be affected by a burst event. In particular, for our choice of α_MRI, T_MRI, and κ_R(q_dust = −3.75, a_max = 200 µm, T = T_MRI) and for M_⋆ = 1 M_⊙, we find that ${\text{Σ}}_{\text{g,min}}~200R_{\text{au}}^{0.75}{\text{g/cm}}^2$. This would suggest that a Minimum Mass Solar Nebula (MMSN) disk with ${\text{Σ}}_{\text{g,disk}}\approx 1700R_{\text{au}}^{-3/2}{\text{g/cm}}^2$ (Weidenschilling 1977b) would feature an inverted density slope interior to 4 au, whereas the burst-prone region would shrink to within 1 au for a less massive disk with ${\text{Σ}}_{\text{g,disk}}\approx 200R_{\text{au}}^{-1}{\text{g/cm}}^2$.

Fig. 3 Pre- (left) and post-burst (right) states of the gas surface density (top) and temperature (bottom) for different disk configurations computed using the method described in Sect. 3. Here, we assumed constant dust opacities of κd = 700 cm2/gdust (including fsubl from Eq. (7)). The mass accretion rate through the outer boundary is approximately ${10}^{-9}{\text{M}}_{\odot }\text{/yr}\times {\alpha }_{\text{DZ}}{\text{/10}}^{-3}\times {\text{Σ}}_{\text{g}0}/200{\text{\hspace{0.17em}g}}^{/}$. Dashed-dotted lines in the top panels have a slope ∝ R0.81—in good agreement with Eq. (11) – and help guide the eye, while in the bottom panels they represent the irradiation temperature.

3.2 Post-burst state

Post-burst, the disk will spend a comparatively longer period in a quiescent state until the conditions for a burst can be fulfilled again. This phase more or less corresponds to the time it takes for the mass in the inner regions affected by the burst to be replenished via viscous accretion from the outer disk, a process happening on timeframes proportional to the viscous timescale with α = α_DZ. This state is not very exciting from a gas dynamics point of view, as it is largely determined by viscous evolution until enough material has been delivered to the inner disk such that a new burst cycle can begin.

The quiescent phase will last as long as the following condition is satisfied everywhere $$\frac{\partial Q_{\text{visc}}}{\partial T}<\frac{\partial \left|Q_{\text{cool}}\right|}{\partial T}.$$(13)

In other words, the disk must cool efficiently enough to prevent runaway heating after a small temperature increase, otherwise a new burst cycle will begin, erasing all “progress” in terms of mass accumulation in the inner disk. Unfortunately, modeling this would require following the hydro- and thermodynamical evolution of the quiescent phase, which can be quite long for low values of α_DZ. Nevertheless, assuming that the density evolves entirely viscously, that is (Lynden-Bell & Pringle 1974) $$\frac{\partial {\text{Σg}}_{}}{\partial t}=\frac{3}{R}\frac{\partial }{\partial R}\left[\sqrt{R}\frac{\partial }{\partial R}\left(\nu {\text{Σg}}_{}\sqrt{R}\right)\right],$$(14)

and that the temperature will respond to the density evolution via a balance among viscous heating, irradiation, and thermal cooling on timescales much shorter than the viscous timescale, we can estimate both the duration of the quiescent phase and the next pre-burst state by iterating between a viscous evolution solver that updates Σ_g and a root-finding algorithm that updates T via Eq. (12), checking if Eq. (13) is satisfied after every pair of updates³. In Appendix A we demonstrate that this method produces results that agree very well to a fully radiation-hydrodynamical model from our simulation set in this work.

In all cases, the criterion in Eq. (13) is first violated at the DZIE (∼0.1 au), which is where the next burst cycle will begin.

However, a clear explanation behind the timescale separation between bursts is still lacking from a theoretical point of view and is beyond the scope of this work. Nevertheless, after computing the pre- and post-burst states for three different disk configurations in Fig. 3, we find that the quiescent phase lasts for ∼3, 20, and 100 kyr for α_DZ = 10⁻³, 3 × 10⁻⁴, and 10⁻⁴, respectively, suggesting an empirical scaling of α^−3/2 for the quiescent phase duration. This could be related to the fact that the viscous timescale scales as t_visc = R²/ν ∝ 1/α, and that low-α models require more mass to accumulate before a burst can occur, to compensate for the weaker viscous heating (Q_visc ∝ α). From the latter point, the minimum surface density for a burst to occur then scales as Σ ∝ α^−1/2 (see also Eq. (9), albeit with α_MRI being used in that context as it refers to the post-burst state). Combining these scalings yields a quiescent phase duration ∝ α^−3/2, in line with our results.

Fig. 4 Radius–time heatmaps of the gas surface density (top), temperature (middle), and maximum dust grain size (bottom) as functions of radius for our fiducial model with αDZ = 10−3. The burst phase is highlighted with insets and marked with white boxes in the main panels, showcasing the emergence of reflares and substructures in the form of rings during the outburst. The rings diffuse away due to viscosity during the long quiescent phase following a burst. The behavior shown here repeats periodically, with the next burst cycle beginning at ∼6.8 kyr of evolution (see also Fig. 5).

3.3 Example pre- and post-burst states

We can now compute the pre- and post-burst states for different values of α_DZ and profiles of Σ_g,disk to showcase a few different disk configurations while accounting for their stability to accretion bursts near the inner edge of their dead zone. Some indicative results are shown in Fig. 3 for constant dust opacities of κ_d = 700 cm²/g_dust, and different combinations of α_DZ, Σ_g,disk, and ϵ. From the figure we can see how disks with lower α_DZ or ϵ tend to accumulate more mass in the inner regions before bursting, whereas more massive disks have more extended regions affected by the burst. Both trends are in agreement with the physical picture described above.

4 Results: Accretion outbursts

In this section we present the results from our radiation hydrodynamical models of the dead zone inner edge. We begin by describing in detail the results from our fiducial setup in terms of the behavior of the gas and dust during an accretion outburst as well as in the post-burst phase. We then analyze the impact of different thermodynamical prescriptions on the burst characteristics and dust evolution. Finally, we present models with different values of α_DZ to showcase the role of viscous diffusion in the burst and post-burst evolution.

4.1 Fiducial model

As viscous evolution is rather quick for α_DZ = 10⁻³, to the point that we can integrate over multiple burst cycles, we choose this value for our fiducial model. This choice has the additional advantage of more closely mimicking the setup in CF24, making a direct comparison easier. For this model, we initialized our disk with power-laws in density with $${\text{Σ}}_{{\text{g}}_0}=200R_{\text{au}}^{-15/14}{\text{g/cm}}^2,{\text{Σ}}_{{\text{small}}_0}=0.3{ϵ}_0{\text{Σ}}_{{\text{g}}_0},{\text{Σ}}_{{\text{large}}_0}=0.7{ϵ}_0{\text{Σ}}_{{\text{g}}_0},$$(15)

with ϵ₀ = 0.01 being the initial dust-to-gas ratio. The temperature is initialized with a simple power-law assuming Q_surf = Q_irr in the dead zone: $$T_0=157R_{\text{au}}^{-3/7}\text{\hspace{0.22em}K}\Rightarrow h_0=0.025R_{\text{au}}^{2/7}.$$(16)

Finally, the gas and dust velocities are set to (nearly) Keplerian for this configuration accounting for viscous drift, with $$u_{R\text{g}}=u_{R\text{d}}=-\frac{3}{2}\frac{\nu }{R},u_{\varphi \text{g}}=R{\text{ΣK}}_{}\sqrt{1-3/2h^2},u_{\varphi \text{d}}=R{\text{ΣK}}_{}.$$(17)

This setup would approximate a constant accretion rate $\dot{M}=-2\pi {\text{Σg}}_{}R{u}_{R_{\text{g}}}=3\pi v{\text{Σg}}_{}\approx {10}^{-9}{\text{M}}_{\odot }\text{/yr}$ in viscous equilibrium in the dead zone, but of course does not account for the direct illumination of the inner rim (see Eq. (8)), the dependence of α on temperature, or dust growth (Birnstiel et al. 2012). As a result of the former two, it is immediately unstable to the thermal instability near the α transition region at R_rim, and results in an accretion burst, albeit with unrealistic initial conditions. For this reason, similar to CF24, we discard the evolution of this first outburst and the viscous evolution that follows it and only analyze data starting from the next burst thereafter, where the dust–gas mixture has had the chance to adapt to the viscously evolving disk, the dust has reached a fragmentation/coagulation equilibrium, and the temperature has adjusted to the thermal equilibrium given by all terms in Eq. (1c).

Figure 4 shows an overview of the results of our fiducial model in terms of the gas density Σ_g, temperature T, and maximum dust grain size a_max as functions of radius and time. The accretion burst phase is highlighted with insets, revealing the formation of substructures in the form of rings and the presence of reflares similar to those documented in CF24. The figure also illustrates the viscous evolution phase of the disk between burst cycles, over which the substructures formed during the burst phase dissolve due to viscous diffusion. Other general features of the inner disk edge are also visible, such as the hot, MRI-active, low-density cavity interior to R_rim ≈ 0.1 au and the relative durations of the burst and quiescent phases, lasting ∼85 yr and ∼3200 yr, respectively. In Appendix E, we show an extended version of Fig. 4 covering the full duration of our simulation, to illustrate the periodicity of the burst cycles.

Assuming that all the mass exiting the domain radially inwards is accreted onto the star, we can compute the accretion rate onto the star Ṁ_acc as $${\stackrel{\cdot }{M}}_{\text{acc}}=-2\pi R{\left[{\text{Σg}}_{}{u}_{R\text{g}}+{\displaystyle \sum }_i{\text{Σ}}_{\text{d},i}{u}_{R\text{d},i}\right]|}_{R=R_{\text{in}}},$$(18)

and the corresponding accretion luminosity L_acc assuming a stellar radius of R_⋆ = 2.5 R_⊙ (T_⋆ = 4200 K) as $$L_{\text{acc}}=\frac{G{M}_{\star }{\stackrel{\cdot }{M}}_{\text{acc}}}{R_{\star }}.$$(19)

We then plot Ṁ_acc, L_acc, and the total stellar luminosity L_tot = L_⋆ + L_acc as functions of time in Fig. 5. We find that both Ṁ_acc and L_acc increase by about three orders of magnitude during the burst phase, reaching peak values of ∼10⁻⁷ M_⊙/yr and ∼1 L_⊙, respectively, before dropping back to quiescent values of ∼10⁻¹⁰ M_⊙/yr and ∼10⁻³ L_⊙. These values are in good agreement with CF24, albeit shifted downwards by a factor of ∼10 due to the lower reference surface density in our model. The star itself reaches a peak luminosity of ≈3.5 L_⊙ during the burst phase, with about half of its total luminosity coming from accretion. We note that modulations in the accretion rate due to the reflares can also be seen in the inset of Fig. 5.

For a more quantitative comparison of the pre- and post-burst states as well as the burst phase itself, we show radial profiles of various gas- and dust-related quantities at different times in Fig. 6⁴. The disk is practically featureless in the pre-burst state (blue), with the exception of a pile-up of dust near the inner edge of the dead zone at ≈0.12 au due to the presence of a pressure maximum induced by the viscosity transition there (see also Flock et al. 2017). It is nevertheless worth noting that the gas surface density interior to ∼1 au is far from a typical power-law profile, instead characterized by a combination of the inherited post-burst profile from the previous cycle and viscous accretion from the outer disk during the quiescent phase. The quiescent phase is terminated with the reignition of the MRI at ≈0.12 au, hinted at by a small kink in the temperature profile (top right).

During the burst phase (orange in Fig. 6), a heating front propagates out to ≈0.9 au, increasing temperatures above T_MRI and even T_subl within that region. This results in both significant dust sublimation due to high temperatures as well as intense fragmentation due to the dramatic increase in α due to the MRI (with α_MRI = 0.1). As a result, the dust-to-gas ratio, maximum grain size, and total opacity plummet within the burst region. As the front recedes due to cooling, dust is allowed to recondense and regrow via coagulation. This process, highlighted in the gray region, is evident by a drop in temperature to ≲500 K and the increase in q_dust to values near −3.2, indicating a non-equilibrium state dominated by coagulation (P+24). A few more reflares occur during the same burst cycle, each propagating outwards less than the previous one (with the second reflare stopping at ≈0.45 au as shown in the figure), until all available mass is accreted and the burst phase is terminated.

The post-burst phase (green in Fig. 6) is characterized by an inverted surface density profile within ∼1 au, as predicted by Eq. (10), albeit rich in radial substructure in both gas and dust, and a cooler temperature profile compared to the pre-burst state. The latter is a consequence of the lower surface densities leading to both less efficient viscous heating (∝ Σ_g) and more efficient radiative cooling $\left(\propto {\text{Σg}}_{-1}^{}\right)$.

Interestingly, the dust size distribution has already recovered to an equilibrium state dominated by fragmentation with q_dust ≈ −3.75, suggesting that dust growth occurs on very short timescales, possibly comparable to the burst duration itself. Nevertheless, as we will show in Sect. 4.2, dust evolution during the burst phase is actually crucial in determining both the radial extent of the burst region as well as the post-burst state. This is evident by the spike in opacity at the very edge of the heat front during the burst phase (at R ≈ 0.45 au in panel d of Fig. 6), which is a direct consequence of dust fragmentation increasing the abundance of small grains with higher opacities. The dust coagulating on timescales comparable to the burst duration is also hinted at by the gradual increase in a_max within the burst region in the bottom panel of Fig. 4 (see, e.g., the region between 0.3–0.8 au at t ≈ 3.4–3.42 kyr).

Finally, we discuss the “S-curve” of the disk at R = 0.3 au in Fig. 1 to illustrate the coevolution of the gas surface density and temperature during a burst cycle. The cycle consists of two “stable” branches (highlighted in gray bands in the figure) corresponding to the quiescent, cold state between the end of a burst cycle until the beginning of the next (bottom), and the hot state during which the MRI is active and the disk is rapidly accreting onto the star (top). Starting at the lower branch, and following the blue curve, the disk evolves as follows (see also the numbered points in Fig. 1):

1. The pre-burst state reaches the critical Σ_g for runaway heating at the DZIE (∼0.1 au).

2. The activation of the MRI and the resulting increase in α lead to viscous transport away from the triggering region, partly transporting mass outwards, and sequentially igniting the MRI in regions that satisfy Σ ≳ Σ_g,min (see Eq. (9)), akin to an outward-propagating burst front. This front passes through the highlighted radius (here R = 0.3 au) as it advances outwards, resulting in a momentary spike in Σ_g at this radius, while T continues to rise rapidly.

3. Temperatures reach T_MRI and quickly exceed T_subl, preventing further heating due to the thermostat effect. The disk enters a hot state of viscous accretion with α = α_MRI.

4. Σ_g and T have dropped enough for dust to recondense and allow the disk to cool efficiently.

5. The burst front recedes past this radius while the disk continues to cool.

6. Additional reflares may trigger at smaller radii (see remaining colors), repeating the process until Σ_g finally reaches its post-burst value and the disk returns to the quiescent state. While the precise condition for the triggering of such reflares is not entirely clear, it is very likely linked to a balance between viscous and cooling timescales during the evolution of the burst (see also Wünsch et al. 2006).

Overall, the behavior of the burst in terms of gas properties is very similar to that shown in CF24, including the presence of multiple reflares, with the addition of a decoupled dust size distribution and its evolution being detailed in this work. The latter shows quite interesting features, with dynamic dust evolution playing a key role in the dynamics of the burst cycle via the modification of the overall opacity due to fragmentation/sublimation, the coagulation of dust on a timescale comparable to that of the burst event itself, and the accumulation of grains on the pressure maxima induced by the burst flares. In the next paragraphs, we will stress the importance of dust coagulation in capturing the correct burst evolution, as well as the effects of in-plane radiative diffusion.

Fig. 5 Top: Accretion rate onto the star (through the inner radial boundary) and accretion luminosity as functions of time for our fiducial model. The inset zooms in on the narrow burst region highlighted in orange. The black line corresponds to exponential growth with an e-folding time of ∼2.2 kyr. Bottom: Total stellar luminosity as a function of time.

Fig. 6 Radial profiles of various gas- and dust-related quantities during the pre-burst (blue), burst (orange), and post-burst (green) phases for our fiducial model: gas surface density (panel a), temperature (b), dust-to-gas ratio (c), Rosseland mean opacity (d), maximum grain size (e), and dust size distribution exponent (f). The gray-shaded region in the burst phase highlights the area where dust is recondensing and recoagulating after the passage of the heating front. The associated movie is available online.

Fig. 7 Post-burst gas surface density (top) and temperature (bottom) profiles for models with different thermodynamical prescriptions: our fiducial model with dust evolution but no in-plane radiative diffusion (blue), a model without dust evolution (orange), and a model with in-plane radiative diffusion (green). Vertical lines mark the radial extent of the burst region in each model. A dashed line denotes the irradiation temperature profile.

4.2 Comparison among thermodynamical prescriptions

In the previous section, we showed the evolution of a burst cycle for our fiducial model, where we included dust fragmentation/coagulation but ignored in-plane radiation transport due to its extreme cost in the post-burst state. Here, we will compare that model against one where dust evolution is omitted, and one where we include radiative diffusion through the disk midplane via Eq. (3d) during a burst cycle.

For the model without dust evolution, we prescribe q_dust = −3.75 and a radial profile of a_max informed by the results of our fiducial model: $$\begin{array}{l}{a}_{\text{in}}=260\text{μm}{\left(\frac{R}{1.26\text{au}}\right)}^{-0.59},a_{\text{out}}=260\text{μm}{\left(\frac{R}{1.26\text{au}}\right)}^{0.38},\hfill \\ a_{\text{max}}^{\prime }(R)=a_{\text{in}}{a}_{\text{out}}/{\left(a_{\text{in}}^{\zeta }+a_{\text{out}}^{\zeta }\right)}^{1/\zeta },\zeta =3.38,\hfill \\ a_{\text{max}}(R)=a_{\text{min}}+\frac{1}{2}(a_{\text{max}}^{\prime }-a_{\text{min}})\left[1+\tanh \left(\frac{R_{\text{au}}-0.11}{0.002}\right)\right].\hfill \end{array}$$(20)

In the above, a_in and a_out are power-laws fitted to the inner and outer disk regions, respectively, joined smoothly around R = 1.26 au. This profile follows closely the equilibrium profile of a_max in the pre-burst state of our fiducial model (see blue line in bottom panel of Fig. 6). The dust fluids otherwise inherit their properties via a_min and a_max in Eq. (20) as if TriPoD were active (i.e., effective grain sizes, Stokes numbers, and velocities, see P+24).

Both models are initialized very close to the pre-burst state of our fiducial model, at t = 3.37 kyr, and evolved through a full burst cycle. The resulting post-burst states in gas surface density and temperature are shown in Fig. 7, with the burst regions marked with vertical lines. Both models feature several reflare events and similar post-burst temperature profiles.

The model without dust evolution (orange curves) results in a weaker burst, propagating out to ≈0.6 au only, compared to ≈0.9 au in the fiducial case (blue), and lasting ≈60 yr compared to ≈85 for the fiducial model. This is a direct consequence of the lower opacity in the burst region compared to the fiducial model, due to the lack of fragmentation increasing the abundance of small grains which would otherwise dominate the opacity and promote heating. As a secondary effect, the post-burst state has roughly 1.5× higher surface densities within ∼1 au compared to the fiducial model, as less mass has been accreted onto the star during the weaker burst. This is also reflected in the ∼1.5× lower accretion rate during the burst (not shown).

As for the model with in-plane radiative diffusion (green curves), the burst region extends to ≈0.8 au, slightly less than in the fiducial case. This is likely due to the radial diffusion of heat away from the burst front, resulting in slightly lower peak temperatures. The latter effect can be seen in the post-burst temperature profile (lower panel of Fig. 7), which features very slightly smoother temperature peaks compared to the other two models. A by-product of this is an overall slower evolution of the burst, resulting in one fewer reflare event. Finally, the temperature profile within the always-MRI-active region interior to R_rim is also different in the model with radiative diffusion.

All in all, we find that the effects of in-plane radiative diffusion are secondary to those of dust evolution during the burst phase, with the former slightly modifying the specifics of the burst evolution but otherwise resulting in a similar post-burst state. In Appendix B, we further show that in-plane radiative diffusion has no effect on the dynamics during the quiescent, viscous post-burst phase. Dust evolution, on the other hand, plays a key role in determining both the radial extent of the burst as well as the total mass accreted onto the star during a burst event, by a factor of ∼1.5 in both aspects in our models. For these reasons, we omit in-plane radiative diffusion for the rest of our models in this work, but always include dynamic dust evolution with TriPoD.

4.3 Different levels of turbulence in the dead zone

Having analyzed the behavior of our fiducial model in detail and demonstrated the importance of dust evolution, we now shift our focus to models with different values of α_DZ to investigate the role of viscous diffusion during and after a burst event. Using the method described in Sect. 3, we can first evaluate whether a model with a given α_DZ would contain an unstable configuration in the first place, and find that models with α_DZ ≥ 8 × 10⁻⁵ are indeed prone to accretion bursts for our choice of disk surface density profile and opacity model. We therefore run models with α_DZ = 3 × 10⁻⁴ and 10⁻⁴ initialized with the pre-burst state computed via Sect. 3, and evolve them through a burst event and several kyr of viscous evolution thereafter.

Figure 8 shows the post-burst gas surface density and temperature profiles for models with different values of α_DZ. As expected, lower values of α_DZ result in more prominent spikes in Σ_g due to less efficient viscous diffusion during the burst evolution. The weaker viscous dissipation is also reflected in the lower temperatures for lower α_DZ values as well, albeit with sharper peaks due to the very high opacities reached at pressure maxima formed during the burst. The burst evolution is otherwise very similar among the different models, which is also expected as its dynamics are mostly governed by the value of α_MRI (see Sect. 3). As a result, the burst region extends to approximately the same radius of ∼1 au in all models, albeit with a weak inverse scaling of the size of the burst region with α_DZ, as the surface density in the pre-burst state is higher for lower α_DZ (see also Fig. 3).

However, lower values of α_DZ result in significantly longer viscous timescales (i.e., quiescent phase durations), weaker viscous diffusion across the substructures formed during a burst event, and lower turbulent velocities between dust grains. As a result, grains grow to larger sizes for lower α_DZ, can accumulate more efficiently on the less diffuse pressure maxima, and have significantly longer to pile up before the next burst event can occur. This is illustrated in Fig. 9, where we plot the same radial quantities as in Fig. 6 for models with different α_DZ values after ∼200/α_DZ yr of viscous evolution in the post-burst state.

It is worth noting that the dust-to-gas ratio (panel c) is notably enhanced at pressure maxima for lower α_DZ values, reaching values of ϵ ≈ 0.03 for α_DZ = 10⁻⁴ compared to ϵ ≈ 0.015 for α_DZ = 3 × 10⁻⁴. Furthermore, for the same α_DZ, the dust grows to sizes of ∼3 mm, transitioning from small- to intermediate-scale turbulence, evident by the increase in q_dust from ≈−3.75 to ≈−3.5 (panel f). The presence of larger grains and their efficient accumulation at long-lived pressure maxima could enable planetesimal formation via the streaming instability (Youdin & Goodman 2005; Johansen & Youdin 2007) in these regions during the quiescent phase, an aspect we will explore in more detail in Sect. 5.

Fig. 8 Post-burst states similar to Fig. 7 for models with different values of αDZ. Lower values of αDZ result in more prominent gas substructures due to the less efficient viscous diffusion.

5 Results: Planetesimal formation in the quiescent phase

In this section, we investigate if the dust concentration resulting from the outbursts is high enough to trigger the formation of planetesimals via the streaming instability (SI). Understanding the threshold for the SI to be active (clumping criterion), that is, quantifying the minimum dust concentration Z = Σ_d/Σ_g and dynamic pebble size (St) necessary to trigger the SI is still an active area of research (Youdin & Goodman 2005; Youdin & Johansen 2007; Johansen & Youdin 2007; Li & Youdin 2021). To investigate the formation of planetesimals in the quiescent phase, we consider two clumping criteria found in recent literature, namely the ones found in Lim et al. (2024) and Lim et al. (2026) (henceforth “L24” and “L25”, respectively). The first one is given by, $$\begin{array}{ll}\log Z_{\text{L24}}({\text{St}}_{\text{SI}},\alpha )\hfill & =0.15{(\log \alpha )}^2-0.24\log {\text{St}}_{\text{SI}}\log \alpha \hfill \\ \hfill & -1.48\log {\text{St}}_{\text{SI}}+1.18\log \alpha \hfill \end{array}$$(21)

for (St, α) ∈ ([10⁻², 10⁻¹], [10⁻⁴, 10⁻³]). The second is given by, $$\log Z_{\text{L25}}({\text{St}}_{\text{SI}})=0.10{(}^{\log }+0.07\log {\text{St}}_{\text{SI}}-2.36$$(22)

and was derived for St ∈ [10⁻³, 1] and in the absence of forced turbulence. Since these criteria are derived for a single dust size, in the following section, we assume St_SI = St(a_max). The planetesimal formation rate can then be calculated as (Miller et al. 2021), $$\frac{{\text{dΣ}}_{\text{plan}}^i}{\text{d}t}=-\zeta {\text{Σd}}_i^{}{\text{St}}^i{\text{ΩK}}_{},$$(23)

where ζ = 0.1 is the planetesimal formation efficiency per settling timescale and 𝒫 is the activation function (Miller et al. 2021) given by: $$𝒫=\frac{1}{2}\left[1+\tanh \left(\frac{\log (Z/Z_{\text{crit}})}{n}\right)\right],Z_{\text{crit}}\in \{Z_{\text{L24}},Z_{\text{L25}}\},$$(24)

where n = 0.03 is a smoothing factor. We compare the aforementioned clumping criteria with the dust concentration in our setups 100 yr after a burst, which can be seen in Fig. 10.

As we can clearly see, the only simulation/criterion pairing that predicts any significant amount of planetesimal formation is the α = 10⁻⁴ setup with the L25 threshold. To actually model the formation of planetesimals in this case, we evolved the disk from the hydrodynamical simulation starting with its after-burst state 700 yr after the end of the burst with an isothermal 1D simulation where we assumed a fixed T(R) (and therefore α(R)) for the duration of the simulation. This approach holds well for a large fraction of the quiescent phase. We used the TriPoDPy code (Kaufmann et al. 2025), which includes viscous gas evolution and uses the TriPoD method for the dust evolution and behaves equivalently to the hydro simulation in the quiescent phase (see Appendix D for a comparison), but is significantly faster as it omits both temperature and momentum evolution. The formation of planetesimals is modeled as a sink term given by Eq. (23) with the L25 clumping criterion. The resulting planetesimal surface density and the ratio of disk and clumping metallicity can be seen in Fig. 11.

After evolving the system for 50 kyr, we form 1.6 M_⊕ in planetesimals, which is split between the fast planetesimal formation at the pileups from the burst and the continuous formation at the inner rim. This is also neatly illustrated by the fact that the SI criterion is significantly surpassed in the bumps at the beginning, and the continuous formation of planetesimals keeps Z/Z_crit ≈ 1 at the inner edge.

We have to note that this formation of planetesimals triggered by outbursts is a very tentative result. Firstly, the criteria for SI were only met for one of the parameter setups (α = 10⁻⁴), considering the most favorable SI criterion. Additionally, there remains a large uncertainty with regard to the appropriate clumping criteria to use, as both the parameter space and setup where they were derived do not fully cover our simulations (e.g., considering the presence of multiple dust sizes). Also, connecting the SI being active to a planetesimal formation rate is still uncertain, with most models, including this work, assuming a constant formation rate per settling timescale (Miller et al. 2021; Schoonenberg et al. 2018; Lau et al. 2022).

Fig. 9 Post-burst states similar to Fig. 7 for models with different values of αDZ. While the gas surface densities are similar among the different models, the effect of lower αDZ in promoting dust growth and accumulation at pressure maxima is evident in panels c, e, and f.

Fig. 10 Comparison of the disk metallicity with clumping criteria from Eq. (21) and Eq. (22) after burst for the different viscosity values.

Fig. 11 Planetesimal surface density (panel a) and the ratio of disk metallicity to critical metallicity (panel b) throughout the TriPoDPy simulation.

6 Discussion

Here we discuss the implications of our results in the context of accretion outbursts and planet formation, as well as the limitations of our models. We also compare and contrast our models with more typical FU Orionis models in the literature.

6.1 Comparison to “traditional” FU Ori models

FU Orionis-type outbursts have been modeled extensively in the literature, with a major difference being the accretion rate and luminosity during the burst phase. In particular, typical FU Ori models feature peak accretion rates of ∼10⁻⁴–10⁻³ M_⊙/yr (e.g., Zhu et al. 2009a, 2010; Bae et al. 2013, 2014; Vorobyov & Basu 2015), resulting in accretion luminosities of ∼10²–10³ L_⊙ for solar-type stars. These values are significantly higher than those found in our models (see Fig. 5), and is mainly due to our choice of a much lower surface density profile. In fact, the aforementioned models leverage the high surface densities to drive the gravitational instability throughout the disk, which in turn supplies mass to the inner disk with an α_DZ that can be as high as 10⁻² (e.g., Zhu et al. 2009a; Vorobyov & Basu 2015), resulting in much higher pre-burst accretion rates and more massive burst events that can extend out to several au (see also dashed line in Fig. 3).

In fact, considering both the period and amplitude of the accretion bursts in our models (every ∼kyr, with L_acc ∼ 3L_acc0), our fiducial model results are unlike those of typical stellar outbursts seen in YSOs (see Fig. 3 in Fischer et al. 2023), but can be easily tuned to match with those of EX Lupi-type objects (Herbig 2007) or even FU-Ori-like objects by adjusting the disk mass and the values of α_DZ and α_MRI accordingly. An exact match with such objects is, however, not the focus of this work.

6.2 Constraints on turbulent dissipation and their implications

As we have seen both in our analytical estimates in Sect. 3 and in our numerical models in Sect. 4, the values of both α_DZ and α_MRI play key roles in determining the characteristics of accretion outbursts near the DZIE by setting, the pre-burst state, the extent of the burst region, and assigning a lifetime to both the quiescent phase and the radial substructures formed during the burst phase. In this work, we have adopted fiducial values of α_DZ and α_MRI following CF24, without actually arguing for the origin of a turbulent mechanism that would yield these values.

In particular, while the value of α_MRI = 0.1 we have adopted is well within the range of turbulent stress reported for the MRI from both ideal and non-ideal MHD models in the literature (e.g., Simon et al. 2011; Bai & Stone 2013a; Flock et al. 2017; Iwasaki et al. 2024), the origin of a non-zero α_DZ in the dead zone is far less clear. Mechanisms that are often invoked in the outer disk such as the vertical shear instability (VSI, Nelson et al. 2013; Flock et al. 2020), an ambipolar-diffusion-mediated MRI (e.g., Turner et al. 2014; Delage et al. 2022), or even the gravitational instability (GI, Toomre 1964; Gammie 2001) are not expected to operate efficiently within the inner 10 au of the disk due to prohibitively long cooling timescales (Lin & Youdin 2015), low ionization fractions (Bai & Stone 2013b), and high values of the Toomre parameter Q (Toomre 1964), respectively.

Possible sources of turbulence within a few au of the central star could include hydrodynamical instabilities such as the convective overstability (COS, Klahr & Hubbard 2014; Lyra 2014) or the zombie vortex instability (ZVI, Marcus et al. 2015, Marcus et al. 2016), both of which, however, are very sensitive to the disk’s thermal relaxation timescale, vertical stratification, and radial density and temperature profiles. As a result, the efficiency of these instabilities in driving adequate turbulence in the dead zone remains to be confirmed with high-resolution non-linear simulations including realistic thermodynamics. Finally, the turbulence driven by the streaming instability itself is expected to be very weak, with α ≲ 10⁻⁵ (Baronett et al. 2024).

Nevertheless, assuming that the bulk of the accretion is driven via a magnetothermal wind at the disk surface layers (Bai & Stone 2013b), all that is really required is some level of turbulent diffusion at the DZIE to enable the “ignition” of the MRI via viscous heating. To that end, Jankovic et al. (2021) and Iwasaki et al. (2024) have shown that the transition from the MRI-active inner disk to the dead zone may be mediated by a radial buffer zone with a nonzero α, which could potentially play the role of the “spark” that ignites the burst cycle. This could be modeled with a combination of imposing a laminar radial flow (e.g., Kimmig et al. 2020) and a different transition region between the active and dead zones via, for instance, a more detailed temperature dependence of α in Eq. (6) informed by the aforementioned works.

Although modeling the triggering of the MRI via a temperature threshold, T_MRI, is commonly used in simulations of outburst mechanisms (e.g., Zhu et al. 2009b; Bae et al. 2013; Kadam et al. 2020; Steiner et al. 2021; Cecil & Flock 2024), a more detailed determination of MRI-active and dead zones can be achieved by coupling the MRI activity to the ionization fraction and non-ideal MHD effects in the inner disk (Dzyurkevich et al. 2013; Mohanty et al. 2018; Jankovic et al. 2021; Delage et al. 2022). This can lead to a different position and shape of the DZIE, which then additionally depends on the structure and strength of the magnetic fields interacting with the disk material. Recent studies of steady-state solutions for large-scale magnetic fields in the inner disk (e.g., Steiner et al. 2025) in combination with tabulated diffusivities of non-ideal MHD effects (e.g., Desch & Turner 2015; Williams & Mohanty 2025) will be used in a forthcoming paper to investigate the influence of a more elaborate MRI activity description on the triggering and evolution of the outburst mechanism analyzed in this work.

It is worth noting that embedded planets will contribute to a significant fraction of the angular momentum transport within the inner disk via the excitation of spiral density waves (Goldreich & Tremaine 1979; Ogilvie & Lubow 2002), effectively driving a non-zero α (Goodman & Rafikov 2001). In a low-viscosity environment such as the dead zone, such planets could even open deep gaps whose edges will be prone to the RWI (Ziampras et al. 2025b), further enhancing the total Reynolds stress within the dead zone. Finally, shock heating by such planets can act as an additional, very efficient source of thermal energy especially in the inner few au of the disk (Rafikov 2016; Ziampras et al. 2020; Rowther et al. 2020; Ono et al. 2025; Okuzumi et al. 2026), which has the potential to eliminate the need for viscous heating altogether to trigger the MRI at the DZIE. The presence of giant, gap-opening planets in the inner disk could also lead to accumulation of gas near the DZIE, possibly leading to burst events of similar origin but different behavior to those studied here (Lodato & Clarke 2004). Investigating these possibilities could be the focus of future work.

6.3 Implications for planet migration

The inverted or otherwise heavily modified surface density profiles formed during a burst event and maintained during the majority of the quiescent phase (see top panels of Fig. 3), the presence of multiple pressure maxima for several kyr after a burst (see Fig. 8), and even the non-power-law-like temperature profiles for sufficiently high α_DZ (bottom panels of Fig. 3), could have significant implications for the migration of low-mass planets embedded in the inner disk. If indeed there is sufficient turbulent diffusion to prevent planets from opening gaps in the gas (e.g., Crida et al. 2006; Duffell 2015; Ziampras et al. 2025b), then such planets will be subject to type-I migration (Goldreich & Tremaine 1980; Ward 1997), a process known to be highly sensitive to the local radial gradients of both surface density and temperature (Paardekooper et al. 2010, 2011).

In this context, pressure bumps could act as planet traps (e.g., Coleman & Nelson 2014; Izidoro et al. 2017), halting inward migration and promoting the growth of planetary embryos via pebble accretion (Lambrechts & Johansen 2012; Lambrechts et al. 2014). Once planets have grown sufficiently massive to no longer migrate in the type-I regime, the DZIE itself would then act as a barrier, preventing further inward migration (Ataiee & Kley 2021; Chrenko et al. 2022). Our analytical prescriptions for the pre- and post-burst states from Sect. 3 can be used to model planet migration in the inner regions of protoplanetary disks qualitatively, without the need for full hydrodynamical simulations.

We stress that these expectations hinge on the assumption that the burst evolution is not strongly affected by an embedded planet, which is not straightforward to assume given that planetary spirals could either open gaps that can modify the local cooling properties of the disk or smear the burst-related pressure bumps, with either process affecting the evolution of the burst. The interplay between accretion outbursts and embedded planets will be the subject of future work.

6.4 Limitations of a 1D model

Our models are focused on the radial evolution of the inner regions of the protoplanetary disk during and after an accretion outburst, and as such are limited to a vertically integrated, axisymmetric (1D) framework. While this approach allows us to capture the key physical processes involved in the burst dynamics and dust evolution, it comes with some simplifying assumptions and limitations to what we can infer from our results.

Possibly the most significant limitation of our 1D approach is the inability to capture nonaxisymmetric instabilities that may arise during the burst phase, in particular the Rossby wave instability (RWI, Lovelace et al. 1999). The RWI is known to develop where the vortensity profile of the disk $\varpi \left(R\right)=\left(\nabla \times u\right)/\text{Σ}\cdot \widehat{z}$ features a local extremum (see also Chang & Youdin 2024), which is easily satisfied around the burst front even during the evolution of the outburst (see also Cecil & Flock 2024; Cecil et al. 2026). Effectively, the development of the RWI would erode the burst front into one or more vortices, weakening its sharpness and possibly the extent to which dust grains can accumulate at the resulting pressure maximum. The Reynolds stress generated by the vortices may also contribute to the overall angular momentum transport during the burst phase and while the vortices persist (Kuznetsova et al. 2022), which can last several thousand orbits at α = 10⁻⁴ (e.g., Rometsch et al. 2021), effectively increasing the value of α_DZ during that time and leading to faster turbulent diffusion of substructures formed during the burst. At the same time, however, the presence of vortices would promote dust accumulation by creating long-lived, nonaxisymmetric dust traps, possibly enhancing planetesimal formation via the streaming instability (e.g., Raettig et al. 2015). The interplay between the RWI, burst dynamics, and dust evolution will be explored in follow-up work (Ziampras et al., in prep.).

Regarding the vertical structure of the disk, our 1D models assume vertical hydrostatic equilibrium and a vertically isothermal temperature profile, that is, ignoring the hot corona of the disk atmosphere. We expect that this is a reasonable approximation for both the temperature structure and the gas–dust coevolution near the midplane, as that is where most of the mass is found and where viscous heating dominates the thermal budget, especially during the burst phase. However, we note that the 2D axisymmetric models of CF24 have shown vertical stirring of gas during the burst phase, which may affect the vertical distribution of dust grains for a brief period of time. The implications of this effect on dust evolution and planetesimal formation remain to be explored in future work.

Finally, it is entirely possible that the assumption of a smooth, power-law-like profile for the flaring angle θ(R) (see Eq. (8)) in the outer disk may not fully hold in reality, due to the presence of shadows cast by the puffed-up, directly illuminated inner rim (see also Dullemond & Monnier 2010; Flock et al. 2025). Such shadows could simply modify the irradiation heating term Q_irr in Eq. (1b), or even trigger additional thermal instabilities in the outer disk (Wu & Lithwick 2021; Melon Fuksman & Klahr 2022; Sudarshan et al. 2026). Nevertheless, given the relatively small radial extent of interest in our models (≲ 5 au), and the central role of viscosity in the burst evolution, we expect that the effects of radial modulations to the irradiation heating would be secondary to those of viscous heating and cooling.

7 Summary

In this work, we have presented vertically integrated, axisymmetric models of stellar outbursts via a thermal instability at the dead zone inner edge, including for the first time a fully integrated treatment of dust coagulation, dust–gas thermal and dynamical coupling, and radiative transfer via surface cooling, in-plane radiative diffusion, viscous heating, and stellar irradiation. While our models are limited to a 1D framework, several key messages emerge from our analysis.

We have developed a semi-analytical method that can be used to compute the pre- and post-burst states of a disk subject to accretion outbursts at the DZIE for a given set of stellar and disk parameters. We then demonstrated that our method yields results in excellent agreement with numerical simulations employing full radiation hydrodynamics and dust evolution. This method can be used to efficiently scan the parameter space of accretion outbursts without the need for expensive numerical simulations.

Dust evolution plays a key role in determining the radial extent of the burst region during an accretion outburst – moreso than radiative diffusion – due to dust fragmentation increasing the abundance of small, opacity-carrying grains. As a result, our models including dust evolution yield up to 1.5× larger burst regions and accreted mass onto the star compared to models without dust evolution, both of which have a noticeable effect on the burst amplitude and post-burst density distributions.

The level of turbulent diffusion in the dead zone influences dust growth and planetesimal formation in multiple ways. Regarding the gas, lower values of α_DZ result in more massive bursts that extend further out in radius, more pronounced and longer-lasting substructures due to the weaker viscous diffusion, and much longer quiescent periods, all of which act in favor of more efficient and longer-lived dust traps. As for the dust, lower α_DZ values lead to lower turbulent velocities and as a result larger grains, which can accumulate more efficiently at pressure maxima, all while further promoting trapping by reducing radial diffusion of dust grains.

Even though the outbursts lead to an accumulation of dust in several rings, it remains questionable whether the SI is triggered in these bumps. In our follow-up simulations, we found that only the most favorable SI criterion (L25) and α = 10⁻⁴ formed any planetesimals.

While our models combine several intertwined physical processes in a self-consistent manner, the limits of a 1D framework prevent us from capturing nonaxisymmetric instabilities such as the Rossby wave instability, which is expected to develop at the burst front during the outburst phase (Cecil et al. 2026). The implications of such instabilities on the burst dynamics and dust evolution will be explored in follow-up work.

Data availability

Data from our numerical models are available upon reasonable request to the corresponding author. Movie associated to Fig. 6 is available at https://www.aanda.org

Acknowledgements

AZ would like to thank Mario Flock and Zhaohuan Zhu for insightful discussions. AZ, TB, and NK acknowledge funding from the European Union under the European Union’s Horizon Europe Research and Innovation Programme 101124282 (EARLYBIRD). Views and opinions expressed are those of the authors only. TB acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. This research was supported in part by grant NSF PHY-2309135 to the Kavli Institute for Theoretical Physics (KITP). All plots in this paper were made with the Python library matplotlib (Hunter 2007).

References

Appendix A Comparison between PLUTO and Sect. 3

Here we compare the pre- and post-burst states obtained with the method described in Sect. 3 against our fiducial model from Sect. 4.1, which includes full radiation hydrodynamics and dust evolution. For this comparison, the dust opacity includes a snapshot of the dust distribution from the fiducial model at t = 3.37 kyr, just before the burst event begins. We therefore set a_max via Eq. (20), q_dust = −3.75, and a constant dust-to-gas ratio ϵ = 0.01. The latter two approximations hold very well during the viscous quiescent phase (see Fig. 9, blue curves).

Fig. A.1 Comparison of pre- and post-burst states in gas surface density (top) and temperature (bottom) obtained with PLUTO (solid lines) and with the method described in Sect. 3 (dashed lines). The two methods show excellent agreement overall.

Figure A.1 summarizes our findings, showing the pre- and post-burst states obtained with PLUTO (solid lines) and with the method outlined in Sect. 3 (dashed lines). We find excellent agreement in the overall behavior, with key differences being the presence of substructures in the PLUTO results due to the dynamic evolution of the burst, and slightly lower density/temperatures at the inner edge of the dead zone when using our method. The latter is most likely due to the change in a_max as the disk evolves viscously, a process not captured by our method, where a_max is fixed to the pre-burst profile. Nevertheless, this effect is minor, and does not affect the overall validity of our approach.

Appendix B In-plane radiative diffusion during the quiescent phase

To justify our choice of omitting in-plane radiative diffusion during the post-burst, quiescent phase in our models, we compared the cooling terms Q_cool and Q_rad (see Eqs. (3c) and (3d)) during different stages of the burst cycle in our fiducial model. The results are summarized in Fig. B.1, where we plot the two terms during the pre-burst, mid-burst, and post-burst phases.

We find that in-plane radiative diffusion is only relevant near the burst fronts formed during the burst phase, exceeding thermal cooling by a factor of 10 or even 100 near the two fronts at R ∼0.9 au and ∼0.5 au, respectively. However, Q_rad is otherwise negligible compared to Q_cool during the quiescent pre- and post-burst phases, with the exception of the outermost bump at ∼0.9 au post-burst, where Q_rad ∼ 0.3 Q_cool. This comparison highlights that in-plane radiative diffusion does have an effect on the disk temperature profile and the propagation of the burst fronts during the burst phase (see also Fig. 7), but also that our choice to omit it during the quiescent phase is well justified.

Fig. B.1 Comparison of the cooling terms Qrad (solid lines) and Qcool (dashed lines) during different stages of the burst cycle for our fiducial model. While in-plane diffusion is comparable to thermal cooling around the burst front during the burst phase (orange), it is negligible otherwise.

Appendix C growpacity: A computationally efficient dust opacity model suitable for coagulation models

growpacity is a Python module designed to provide temperature-dependent Rosseland and Planck mean opacities for a dynamically evolving dust size distribution, particularly suited for dust coagulation models such as TriPoD. To do so, it provides an interface to optool (Dominik et al. 2021), a tool that can compute accurate opacities for a given dust composition and size distribution, and then precomputes a grid of mean opacities as functions of temperature T, maximum grain size a_max, and size distribution exponent q_dust. These opacity tables are quite lightweight (typically ≲ 1 MB), and can be efficiently interpolated over during a hydrodynamical simulation to provide up-to-date opacities as the dust size distribution evolves.

For a given grain composition and assuming a grain size distribution characterized by a_min, a_max, and q_dust (see Eq. (4)), the optool package (Dominik et al. 2021) can compute the absorption and scattering opacities κ_abs(ν) and κ_sca(ν) (in cm²/g) as well as the asymmetry factor g(ν) (Henyey & Greenstein 1941) over a frequency grid ν. This calculation is done by default using the Distribution of Hollow Spheres method (DHS, Min et al. 2005). The Rosseland and Planck mean opacities κ_R and κ_P can then be computed as $$\begin{array}{l}{\kappa }_{\text{R}}(T)=\frac{{{\displaystyle \int }}_0^{\infty }{u}_{\nu }(T)\text{d}\nu }{{{\displaystyle \int }}_0^{\infty }{\left[{\kappa }_{\text{abs}}+(1-g){\kappa }_{\text{sca}}\right]}^{-1}{u}_{\nu }(T)\text{d}\nu },{u_{\nu }(T)=\frac{\text{d}{B}_{\nu }}{\text{d}T}|}_T,\hfill \\ {\kappa }_{\text{P}}(T)=\frac{{{\displaystyle \int }}_0^{\infty }{\kappa }_{\text{abs}}{B}_{\nu }(T)\text{d}\nu }{{{\displaystyle \int }}_0^{\infty }{B}_{\nu }(T)\text{d}\nu },\hfill \end{array}$$(C.1)

where B_ν(T) is the Planck function at temperature T. By fixing the grain composition and a_min, optool can be used to calculate the absorption and scattering opacities over a grid of a_max and q_dust. Equation C.1 can then be used to compute and tabulate κ_R and κ_P over a_max, q_dust, and T. We found that sampling a_max with 13 points per decade between 0.1 µm and 10 cm, q_dust with 9 points between −4.5 and −2.5, and T with 300 points between 1–3000 K results in accurate mean opacities during interpolation while keeping the size of the opacity tables to merely ∼ 300 kB. A verification test comparing the mean opacities computed with optool directly against those obtained via interpolation in the precomputed tables is shown in Fig. C.1, demonstrating very good agreement overall.

Fig. C.1 Comparison between exact calculations with optool and their interpolated counterparts for two different combinations of amax and qdust, representing a growing (fully grown) distribution with amax = 5 µm (amax = 5 mm) and qdust = −3.2 (qdust = −3.6). These values have been chosen to lie as far as possible from the sampled grid points, to best test the interpolation scheme. Overall, the interpolated opacities agree very well with the exact calculations, with deviations of < 7%.

To evaluate the opacity for a given q_dust, a_max, and T, we can interpolate within the 3D tables of κ_R(q_dust, a_max, T) and κ_P(q_dust, a_max, T). We choose to interpolate for log κ as a function of q_dust, log a_max, and log T, with regular sampling in this space (i.e., logarithmic spacing for a_max and T). This works best when the mean opacities follow a power-law with respect to temperature κ ∝ T^b ⇒ log κ ∝ b log T, which is a reasonable approximation and holds especially well for small grains (e.g., Bell & Lin 1994; Semenov et al. 2003). It also ensures that the interpolated opacities are always positive in case extrapolation would be needed, although by default we clamp the input parameters to the bounds of the precomputed tables as a safeguard.

We use a trilinear interpolation scheme, which is fast and simple to implement, and takes advantage of the fact that the arrays q_dust, log a_max, and log T are sorted and regularly spaced to efficiently locate the indices of the grid points that surround the point of interest. For each array x ∈ {q_dust, log a_max, log T} and for a target value x_t, we first find the index i such that x_i ≤ x_t < x_i+1 as i = ⌊(x_t − x₀) ∆x⁻¹⌋, where x₀ is the first (smallest) value in the sampling space and ∆x = x_i+1 − x_i is the (constant) sampling spacing. As this information is known a priori, this reduces the complexity of finding the required indices from 𝒪(log N) to 𝒪(1), where N is the number of grid points in x, and is especially efficient for larger grids. The interpolation scheme, implemented in both C and Python, is included in the growpacity package. An example of the interpolated mean opacities as functions of q_dust, a_max, and T is shown in Fig. C.2.

Fig. C.2 Heatmaps of interpolated Rosseland mean opacities as functions of pairs of parameters, with the third parameter otherwise fixed at a representative value.

We underscore that growpacity does not provide a new or more realistic dust opacity model, but rather one suitable for use in coagulation models, where dust densities and distributions can vary as a function of position and time—the applicability of the model depends on entirely user-defined choices. The method can be easily extended to include gas opacities (e.g., Semenov et al. 2003; Malygin et al. 2014) and prescriptions for the sublimation of dust species (e.g., Isella & Natta 2005), for a more complete opacity model in regimes where gas opacities are significant.

Appendix D Comparing TriPoDPy with PLUTO simulations during the quiescent phase

To ensure that the evolution in the quiescent phase with TriPoDPy matches the full hydro simulations, we ran a test case evolving the α = 10⁻⁴ simulation from t = 8 kyr to 20 kyr without planetesimal formation and compared the final state of the two simulations, which is depicted in Fig. D.1. As we can see, the gas surface density and the maximal grain size match perfectly, whereas there is a slight difference in the total dust surface density, which can be attributed to the fact that the fluids don’t exchange momentum in TriPoDPy and the flux through the MRI-active front is different in the two codes, which is also shown by the dust size distribution exponent in the cavity. For the purposes of investigating the planetesimal formation in the quiescent phase, the simulations behave the same, as the major differences are in the cavity, which does not contribute to the formation of planetesimals in any case.

Fig. D.1 Comparison between the state of the simulation at t = 20 kyr in TrPoDPy and PLUTO showing the gas surface density (top left), the total dust surface density (top right), the maximal grain size (bottom left) and the dust size distribution exponent (bottom right).

Appendix E Full duration radius–time heatmaps for the fiducial model

In Fig. E.1 we show the radius–time heatmaps of gas surface density, temperature, and dust surface density for the entire duration of the simulation of our fiducial model. The figure illustrates the clear periodicity of the burst cycle process, with the main features being the same as those shown in Fig. 4.

As shown by the features at the leftmost edge of Fig. E.1, the initial conditions are immediately unstable to the condition in Eq. (13) and a burst cycle begins at t = 0. This is simply a byproduct of the initial conditions in Eq. (15), which do not reflect the pre-burst state calculated in Sect. 3 and shown in Fig. 3. For this reason, as explained in Sect. 4.1, we discard this first outburst cycle and limit our analysis to t ≳ 3 kyr (as in Fig. (4)). We note that this is not the case with the models in Sect. 4.3, as there we initialize the disk appropriately following the method in Sect. 3 and verified using our fiducial model in Appendix A.

Fig. E.1 Radius–time heatmaps similar to those in Fig. 4, spanning the entire duration of the simulation.

All Figures

Fig. 1 “S-curve” of the gas surface density versus temperature at R = 0.3 au for our fiducial model, discussed in detail in Sect. 4.1. 1: runaway heating is triggered at the DZIE. 2: the burst front passes through. 3: dust sublimation caps temperature during burst. 4: dust recondensation; 5: dust cooling. Gray bands highlight the stable branches corresponding to viscous evolution in the quiescent (bottom) and burst (top) phases. Different colors denote reflares during the same burst cycle.

In the text

	Fig. 2 Sample dust size distribution, represented with a truncated power-law with qdust = −3.5 and reconstructed with the TriPoD method of Pfeil et al. (2024).
In the text

Fig. 3 Pre- (left) and post-burst (right) states of the gas surface density (top) and temperature (bottom) for different disk configurations computed using the method described in Sect. 3. Here, we assumed constant dust opacities of κd = 700 cm2/gdust (including fsubl from Eq. (7)). The mass accretion rate through the outer boundary is approximately ${10}^{-9}{\text{M}}_{\odot }\text{/yr}\times {\alpha }_{\text{DZ}}{\text{/10}}^{-3}\times {\text{Σ}}_{\text{g}0}/200{\text{\hspace{0.17em}g}}^{/}$. Dashed-dotted lines in the top panels have a slope ∝ R0.81—in good agreement with Eq. (11) – and help guide the eye, while in the bottom panels they represent the irradiation temperature.

In the text

Fig. 4 Radius–time heatmaps of the gas surface density (top), temperature (middle), and maximum dust grain size (bottom) as functions of radius for our fiducial model with αDZ = 10−3. The burst phase is highlighted with insets and marked with white boxes in the main panels, showcasing the emergence of reflares and substructures in the form of rings during the outburst. The rings diffuse away due to viscosity during the long quiescent phase following a burst. The behavior shown here repeats periodically, with the next burst cycle beginning at ∼6.8 kyr of evolution (see also Fig. 5).

In the text

	Fig. 5 Top: Accretion rate onto the star (through the inner radial boundary) and accretion luminosity as functions of time for our fiducial model. The inset zooms in on the narrow burst region highlighted in orange. The black line corresponds to exponential growth with an e-folding time of ∼2.2 kyr. Bottom: Total stellar luminosity as a function of time.
In the text

Fig. 6 Radial profiles of various gas- and dust-related quantities during the pre-burst (blue), burst (orange), and post-burst (green) phases for our fiducial model: gas surface density (panel a), temperature (b), dust-to-gas ratio (c), Rosseland mean opacity (d), maximum grain size (e), and dust size distribution exponent (f). The gray-shaded region in the burst phase highlights the area where dust is recondensing and recoagulating after the passage of the heating front. The associated movie is available online.

In the text

Fig. 7 Post-burst gas surface density (top) and temperature (bottom) profiles for models with different thermodynamical prescriptions: our fiducial model with dust evolution but no in-plane radiative diffusion (blue), a model without dust evolution (orange), and a model with in-plane radiative diffusion (green). Vertical lines mark the radial extent of the burst region in each model. A dashed line denotes the irradiation temperature profile.

In the text

	Fig. 8 Post-burst states similar to Fig. 7 for models with different values of αDZ. Lower values of αDZ result in more prominent gas substructures due to the less efficient viscous diffusion.
In the text

	Fig. 9 Post-burst states similar to Fig. 7 for models with different values of αDZ. While the gas surface densities are similar among the different models, the effect of lower αDZ in promoting dust growth and accumulation at pressure maxima is evident in panels c, e, and f.
In the text

	Fig. 10 Comparison of the disk metallicity with clumping criteria from Eq. (21) and Eq. (22) after burst for the different viscosity values.
In the text

	Fig. 11 Planetesimal surface density (panel a) and the ratio of disk metallicity to critical metallicity (panel b) throughout the TriPoDPy simulation.
In the text

	Fig. A.1 Comparison of pre- and post-burst states in gas surface density (top) and temperature (bottom) obtained with PLUTO (solid lines) and with the method described in Sect. 3 (dashed lines). The two methods show excellent agreement overall.
In the text

	Fig. B.1 Comparison of the cooling terms Qrad (solid lines) and Qcool (dashed lines) during different stages of the burst cycle for our fiducial model. While in-plane diffusion is comparable to thermal cooling around the burst front during the burst phase (orange), it is negligible otherwise.
In the text

Fig. C.1 Comparison between exact calculations with optool and their interpolated counterparts for two different combinations of amax and qdust, representing a growing (fully grown) distribution with amax = 5 µm (amax = 5 mm) and qdust = −3.2 (qdust = −3.6). These values have been chosen to lie as far as possible from the sampled grid points, to best test the interpolation scheme. Overall, the interpolated opacities agree very well with the exact calculations, with deviations of < 7%.

In the text

	Fig. C.2 Heatmaps of interpolated Rosseland mean opacities as functions of pairs of parameters, with the third parameter otherwise fixed at a representative value.
In the text

	Fig. D.1 Comparison between the state of the simulation at t = 20 kyr in TrPoDPy and PLUTO showing the gas surface density (top left), the total dust surface density (top right), the maximal grain size (bottom left) and the dust size distribution exponent (bottom right).
In the text

	Fig. E.1 Radius–time heatmaps similar to those in Fig. 4, spanning the entire duration of the simulation.
In the text