© The Authors 2025

1. Introduction

In protoplanetary disks, snowlines mark the locations where volatile ice sublimates. This process dictates volatile distributions and alters grain properties in the disk (e.g., Öberg et al. 2023). In this context, outburst systems are ideal laboratories for studying the effect of snowlines. Due to an abrupt increase in accretion, these young stellar objects brighten by many orders of magnitude, heating up the disk for tens to hundreds of years and enabling large-scale sublimation of ices (e.g., Fischer et al. 2023). In particular, the snowline is pushed to ∼10–100 au, bringing it within reach of facilities such as the Atacama Large Millimeter/submillimeter Array (ALMA).

V883 Ori is one of the best-studied FU Orionis outburst systems. Cieza et al. (2016) identified a “dark annulus” in ALMA band 6 continuum, suggesting that its water snowline lies at ≈40 au. This shoulder-like feature in the intensity profile was initially attributed to the accumulation of small grains inside the snowline, produced by the disintegration of dust grains after ice sublimation (Saito & Sirono 2011; Aumatell & Wurm 2011). However, multiwavelength analysis has revealed continuity in the spectral index across the intensity shoulder, implying no abrupt change in the grain physical properties by ice sublimation (Houge & Krijt 2023; Houge et al. 2024). The exact snowline location is further complicated by the detection of H${}_2^{18}$O (Tobin et al. 2023) and indirect tracers such as methanol or HCO⁺ (van ’t Hoff et al. 2018; Leemker et al. 2021), which suggest a snowline location ≳80 au.

In this Letter, we demonstrate that these observations can be consistently explained by a model in which the snowline retreats along with the decline of bolometric luminosity. At the inward-moving snowline front, the intensity shoulder arises from the increase in surface density of solids due to ice recondensation, while in the outer disk the sluggish recondensation of vapor preserves a relic of a previously hot state.

2. Model

We set up the disk at the outburst active phase, when the bolometric luminosity reaches its peak value (see Sect. 2.2), assuming the dust-vapor mixture is in hydrostatic balance and thermal equilibrium. The main focus of this work is to model the subsequent dimming phase, when L_bol declines and vapor recondenses onto grain surfaces. A key assumption is that dimming proceeds rapidly (∼200 yr; see Sect. 2.2) so that dynamic processes (gas advection, pebble drift, and diffusion) can be neglected.

2.1. Disk structure

The gas surface density profile is taken as a power law,

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{g}}(r)={\mathrm{\Sigma }}_c{\left(\frac{r}{r_c}\right)}^{-\gamma }.\end{array}$$(1)

Following the fit of dust thermal emission by Cieza et al. (2018), we take r_c = 31 au and Σ_c = 160 g cm⁻², while γ is a treated as a parameter given its poor constraint (van ’t Hoff et al. 2018). The surface density of pebbles is taken as

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{p}}=Z_{\mathrm{peb}}{\mathrm{\Sigma }}_{\mathrm{g}}(r_c){\left(\frac{r}{r_c}\right)}^{-{\gamma }_{\mathrm{peb}}},\end{array}$$(2)

where Z_peb is the metallicity at the reference position. Beyond the snowline, pebbles are assumed to be 50% refractory and 50% water ice by mass, with the material densities of constituent species following DSHARP (Birnstiel et al. 2018).

Following Wang et al. (2025) we constructed a 2D model in the disk’s R − z plane, where a standard viscous disk with α = 10⁻³ (Shakura & Sunyaev 1973) is assumed. For solids, we adopted an Mathis-Rumpl-Nordsieck (MRN) size distribution (Mathis et al. 1977), which is appropriate when dust growth is limited by fragmentation (Birnstiel et al. 2011). The size distribution is characterized by an upper bound, s_max, corresponding to a midplane Stokes number, St₀, a lower bound of s_min = 0.1 μm, and is assumed to be unchanged by the outburst (Houge et al. 2024). Pebbles of different sizes are vertically distributed according to their scale height, H_p, i(H_g, St_i, 0, α), where H_g is the gas pressure scale height and St_i, 0 the Stokes number of pebble i at the midplane (see Appendix B).

The vapor condensation rate (ℛ_cond) is determined by its saturation condition and the total surface area of particles (e.g., Ros & Johansen 2013; Schoonenberg & Ormel 2017). Specifically, the vapor density at any point in the disk reduces at a rate of

$$\begin{array}{c}\hfill -\frac{\mathrm{d}{\rho }_{\mathrm{vap}}}{\mathrm{d}t}={\mathcal{R}}_{\mathrm{cond}}=⟨\pi s_{\mathrm{p}}^2{n}_{\mathrm{p}}⟩v_{\mathrm{th},\mathrm{vap}}{\rho }_{\mathrm{vap}}(1-\frac{P_{\mathrm{eq}}}{P_{\mathrm{vap}}}),\end{array}$$(3)

where $⟨\pi s_{\mathrm{p}}^2{n}_{\mathrm{p}}⟩$ is the surface area of pebbles per unit volume (Appendix B, equation (B.3)).

2.2. Luminosity during outburst dimming phase

We modeled the dimming of V883 Ori following the published luminosities, accounting for observational uncertainties arising from dust extinction and instrumental saturation (see Appendix A). Specifically, the time-dependent bolometric luminosity reads

$$\begin{array}{c}\hfill L_{\mathrm{bol}}(t)=L_{\star }+(L_{\mathrm{pk}}-L_{\star })exp[-0.5{\left(\frac{t-t_{\mathrm{beg}}}{50\mathrm{yr}}\right)}^2],\end{array}$$(4)

where L_⋆ is the stellar luminosity (6 L_⊙; Cieza et al. 2016). We adopted a peak luminosity of L_pk = 4000 L_⊙ (comparable to that used in Leemker et al. 2021), which pushes the water snowline to 80–100 au. The parameter t_beg indicates the year when the outburst begins to fade. As is shown in Fig. 1, choosing t_beg = 1880 A.D. makes the dimming curve consistent with the lower limit of the estimated L_bol of V883 Ori (400 L_⊙, Strom & Strom 1993; and 218 L_⊙, Furlan et al. 2016). By shifting t_beg forward by 40 yr (A.D. 1920), the curve reaches the upper limit (647 L_⊙) given by Liu et al. (2022). Recently, Carvalho et al. (2025) fit the spectrum at 0.4–4.2 μm and obtained an accretion luminosity of 458 L_⊙, lying well within the assumed range. Through equation (4), the model output is only a function of t′=t − t_beg, the time since L_bol starts to decline.

Fig. 1. Evolution of the bolometric luminosity of V883 Ori during the outburst dimming phase. The symbols denote measurements from different studies, with horizontal bars showing the time range of adopted data. The blue band represents the adopted luminosity curve, which carries an uncertainty of ∼40 yr due to the uncertainties in estimating Lbol.

2.3. Temperature model

We adopted a two-stream radiation transfer (2sRT) method to calculate the temperature (e.g., Hubeny 1990). The heating sources at each grid cell are due to irradiation (q_irr), viscous dissipation (q_vis), and latent heat exchange (q_latent). We computed q_irr from the bolometric luminosity (equation (4)), while the viscous heating rate follows from the α-disk viscosity ($q_{\mathrm{vis}}=\frac{9}{4}\rho \alpha c_{\mathrm{s}}^2\mathrm{\Omega }$) and the latent heating rate is proportional to the condensation rate (equation (3), see Wang et al. 2025). In the outburst active phase, the disk adjusts to a temperature and density structure determined by the bolometric luminosity, L_pk, assuming hydrostatic balance. In the dimming phase, following Flock et al. (2017), we let the disk adapt to a new equilibrium temperature obtained by the 2sRT (T_eq) over a thermal relaxation timescale, dT/dt = (T_eq − T)/t_relax (see Appendix C for the calculation of t_relax). The Rosseland mean opacity per unit gas mass, κ_R – a key quantity that determines t_relax – is a free model parameter. In our simulation (see Sect. 3), due to the low metallicity and disk viscosity, the disk temperature is determined by irradiation, with q_vis and q_latent contributing little.

2.4. Fitting ALMA observations

The intensity shoulder in V883 Ori is identified in both ALMA band 6 (0.038″, project code: 2015.1.00350.S, PI: Lucas Cieza) and band 7 (0.028″, project code: 2016.1.00728.S, PI: Lucas Cieza) high-resolution continuum images. We converted the continuum images into radial profiles by azimuthally averaging the data within concentric annuli of radial width 1/5 of the respective beam size. To fit the continuum profiles, we generated synthetic intensity profiles in both bands by conducting 2sRT with profiles for the dust temperature and density taken from the simulation results. The 2sRT was performed in a plane-parallel manner in the disk’s vertical direction, with the optical depth corrected for disk inclination i = 38.3° (Cieza et al. 2016). Pebble opacities for different sizes and ice fractions were computed with optool (Dominik et al. 2021), adopting the DSHARP optical constants (Birnstiel et al. 2018). After being convolved with the beam size, the synthetic intensity profiles were compared with the continuum data at 45–85 au, where the intensity shoulder lies.

Tobin et al. (2023) measured the column density of water isotope H${}_2^{18}$O in band 5 (0.126″). Given the limited data sensitivity, we fit only the mean column density between 80–120 au in our simulations to the observed value, assuming an isotopic ratio ¹⁸O/¹⁶O = 1/560 (Wilson & Rood 1994). In reality, the observed water vapor abundance may be influenced by processes beyond simple recondensation (Sect. 3.2). Therefore, we introduced an additional free parameter, f_cond, such that −dρ_vap/dt = f_condℛ_cond.

To incorporate observational constraints from both continuum and water emission, we applied a log-normal likelihood for each dataset. Each likelihood was then weighted and combined by the number of beams covered by the area of interest (see equations (E.1) and (E.2)).

3. Results

We conducted Markov chain Monte Carlo (MCMC) analysis with emcee (Foreman-Mackey et al. 2013) to obtain the posterior distribution (see Appendix E). The parameter’s posterior values are given in Table 1.

3.1. Emergence of intensity shoulder

The best-fit results are shown in Fig. 3, where the synthetic intensity profiles agree well with the dual-band data, especially capturing the intensity shoulder feature at 50–70 au (gray-shaded region). To understand this, we plot the time evolution of the corresponding surface density profiles in Fig. 2. As the disk cools during the outburst dimming phase, vapor gradually recondenses on pebbles, creating a bump in surface density at the condensation front, which appears in continuum as the intensity shoulder. Matching the observed position of the shoulder requires the snowline to retreat to ≈60 au at the observational epoch (t′ = 121 yr). This places a stringent constraint on the thermal relaxation time, or IR opacity, κ_R, of the disk. If the disk had cooled more rapidly (lower κ_R), the snowline would already have retreated inside 50 au, shifting the shoulder much closer to the star, within the expected observational time.

Fig. 2. Time evolution of the surface density profiles of the best-fit model. Silicate surface density remains constant with time. The thick line denotes the best fit, t′ = 121 yr.

Fig. 3. Comparison of the synthetic intensity profiles with ALMA continuum. Crosses show ALMA data adopted from Houge et al. (2024), spaced according to the beam sizes; shaded regions indicate the rms error. Solid lines represent the best-fit retreating snowline model, while dotted lines correspond to the best-fit static model. Enhanced emission toward the inner regions arises likely from intense viscous heating (Alarcón et al. 2024), which is not accounted for in the model.

Simultaneously fitting the dual-band continuum yields St₀ ≈ 0.01, corresponding to a maximum grain size of around a centimeter. The need to have centimeter-sized particles in the disk follows from the low spectral indices at 50–70 au (α = 2.5 − 3.0), a result consistently found across ALMA band 3-7 in V883 Ori by Houge et al. (2024).

3.2. Extended water emission

The mean column density of H${}_2^{18}$O at 80–120 au, as estimated from our models assuming optically thin emission, is presented in Fig. 4. In all simulated cases, the water abundance drops sharply within the first 110 years, and transitions to a slowly evolving stage afterward. This is due to different condensation timescales in the disk’s midplane and the upper layers (see Appendix D). After the snowline retreats, vapor in the midplane recondenses rapidly. Conversely, in the disk’s upper layer, where pebble densities are much lower, the recondensation time can reach ∼10² − 10³ yr (Visser et al. 2015; Rab et al. 2017). Therefore, due to the lag caused by the finite recondensation time, the snowline inferred from molecular emission would consistently appear further away than its actual position during the outburst dimming phase. Our modeling shows that an enhanced condensation rate (f_cond ≈ 4) is required to match the inferred water vapor abundance, which may indicate a shallower size distribution and/or porous grains (Birnstiel 2024), both of which increase the total dust surface area at the disk atmosphere. In addition, the destruction of water by photodissociation, heterogeneous nucleation, and different outburst durations may also influence the observed gas-phase water abundance (Fischer et al. 2023; Ros & Johansen 2024; Romero-Mirza et al. 2024).

Fig. 4. Average column density of H${}_2^{18}$O at 80–120 au under different condensation rates (fcond). The vertical gray band indicates the time range when observing the intensity shoulder. The dashed black line denotes the average column density derived from ALMA band 5 data by Tobin et al. (2023) with 1σ uncertainty shaded in gray. Dotted gray line denotes the best-fit value in static snowline model.

3.3. The need for a retreating snowline model

A retreating snowline naturally explains the distant intensity shoulder and the extended water emission, relics of the earlier outburst. In contrast, a static snowline model fails to explain the observations (Fig. 3 and Fig. E.2). To illustrate this, we performed another MCMC simulation, where the disk structure is set by a constant L_bol (treated as a free parameter) rather than evolving with the luminosity curve. In other words, the temperature and vapor equilibrium instantly adapt to the stellar luminosity. Because of the instantaneous adjustment, the static model fails to maintain a sufficient amount of vapor in the outer disk: the best-fit water number density is ≈7σ away from the observed value (Fig. 4). Furthermore, in the static model, the snowline is located at ≈30 au, failing to reproduce the shoulder feature at 50–70 au in the observed intensity profile.

4. Conclusions and discussions

We investigated V883 Ori’s continuum and water line emission with a retreating snowline model. The key findings are:

1. The dimming stage of a stellar outburst results in a receding water snowline. Vapor condensation at the snowline front creates a bump in surface density, marking the transition from icy to ice-free pebbles.

2. Vapor in the disk upper layer takes ∼10² − 10³ yr to recondense, leading to substantial amounts of water vapor in the outer disk, even after the snowline has retreated inward. The observed H₂¹⁸O outside the current snowline region (50–70 au) is a leftover of a much hotter past phase.

3. The retreating snowline model simultaneously reproduces the intensity shoulder feature in ALMA band 6/7 continuum and the extended water emission. In contrast, a static model fails to maintain the snowline at 50–70 au under the present-day luminosity.

In our modeling, both the disk cooling time (t_relax), which controls the snowline retreat speed, and the millimeter opacity, which sets the continuum emission, depends on dust opacities. Thus, the thermal history of outburst systems constrains dust properties. For example, the retreating snowline model requires κ_R ≈ 4.6 cm² g⁻¹, whereas a self-consistent calculation from the DSHARP opacity yields κ_R ≈ 0.4 cm² g⁻¹. This mismatch implies excess IR opacity relative to millimeter opacities, possibly arising from a non-MRN size distribution, carbon-rich small grains (Woitke et al. 2016), or porous dust (Zhang et al. 2023). Under the assumed dimming curve of luminosity, our model predicts that the intensity shoulder will move by 10 au over 25 years (Fig. 2), a trend testable with future observations. With ALMA’s growing sample of outburst systems at high resolution, we expect the intensity shoulder to emerge as a ubiquitous feature of large-scale water phase transitions in disks.

Acknowledgments

The authors thank for useful discussion with Vardan Elbakyan, Shangjia Zhang, Hanpu Liu. This work is supported by the National Natural Science Foundation of China under grant Nos. 12233004 and 12473065.

References

Appendix A: Bolometric luminosity measurements

The onset of the outburst in V883 Ori was not directly observed. But it likely starts before A.D. 1888, when the reflection nebula IC 430 was observed to be illuminated by the outburst (Pickering 1890), as argued by Strom & Strom (1993). By fitting the spectral energy distribution, the bolometric luminosity of V883 Ori was estimated to be 400 L_⊙ (Strom & Strom 1993) in the 1980s (Allen et al. 1975; Nakajima et al. 1986) and had significantly decreased to 218 L_⊙ (Furlan et al. 2016; Connelley & Reipurth 2018) in about 20 years. Given the typical month-to-year rise time of FUors (e.g., Fischer et al. 2023), it is very likely that V883 Ori has been in the dimming phase since A.D. 1888.

However, estimating the bolometric luminosity of V883 Ori is fraught with substantial uncertainty. As a young stellar object, it suffers severe extinction and reddening at near-infrared wavelengths from dusty proto-stellar envelopes. In the mid-infrared, where extinction is less severe, observations are hindered by heavy saturation in modern surveys like WISE (Contreras Peña et al. 2020). Liu et al. (2022) derived a much higher bolometric luminosity of ≈647 L_⊙ by accounting for the reddening. Unlike most works adopting discrete photometric data, recently, Carvalho et al. (2025) fits the medium resolution spectrophotometry within a range of 0.4–2.5 μm and found excellent match to even broader spectrum across 0.4–4.2 μm. They yielded an accretion luminosity of 458 L_⊙, which lies well inside the assumed luminosity curve (Fig. 1).

Appendix B: Dust model

At each radius, we modeled the pebble size distribution by constructing several size bins spanning radii from s_max and s_min. For small pebbles (e.g., s_p < 0.1 cm), logarithmic spacing was used, whereas for larger ones, the size bins are more evenly spaced. The pebble surface density for each size bin (Σ_p, i) was then obtained by dividing Σ_p into these size bins according to the MRN distribution $\mathrm{d}{N}_{\mathrm{p}}/\mathrm{d}{s}_{\mathrm{p}}\propto s_{\mathrm{p}}^{-3.5}$. For each size bin i, its Stokes number (St_i) is obtained by scaling the pebble size at the bin center (s_p, i) with s_max, assuming that the pebbles follow the Epstein drag law: i.e., St_i, 0 = St₀s_p, i/s_max. Then the 2D pebble density distribution is obtained by adopting normal distribution in the vertical direction:

$$\begin{array}{c}\hfill {\rho }_{\mathrm{p},i}(z)={\rho }_{\mathrm{p},0}exp[-\frac{1}{2}{\left(\frac{z}{H_{\mathrm{p},i}}\right)}^2].\end{array}$$(B.1)

Here H_p, i is the scale height of pebbles in the size bin i, following Dubrulle et al. (1995):

$$\begin{array}{c}\hfill H_{\mathrm{p},i}(z)=H_{\mathrm{g}}{[1+\frac{{\mathrm{St}}_{i,0}}{\alpha }\left(\frac{2H_{\mathrm{g}}^2}{z^2}\right)(exp\left(\frac{z^2}{2H_{\mathrm{g}}^2}\right)-1)]}^{-1/2}.\end{array}$$(B.2)

In equation (B.2), the vertical variation in the Stokes number due to gas stratification was accounted for, ${\mathrm{St}}_i(z)={\mathrm{St}}_{i,0}{\rho }_{\mathrm{g}}(z)/{\rho }_{\mathrm{g}}(0)\approx {\mathrm{St}}_{i,0}exp(z^2/2H_{\mathrm{g}}^2)$. The gas scale height, H_g, was determined in the outburst active phase, after iterating between the gas density distribution, assuming hydrostatic balance, and the temperature structure calculated from 2sRT. Given Σ_p, i, the midplane pebble density, ρ_p, 0, can be obtained by integrating the vertical density distribution following equation (B.1).

Recently, Houge et al. (2024) preformed dust evolution modeling and demonstrated that, to match the observed spectral index in V883 Ori, pebble sizes should not evolve appreciably during the swift outburst event. Therefore, we maintain a constant pebble size distribution throughout the outburst dimming phase.

Given the size distribution, the total surface area of pebbles per unit volume can be obtained by summing over all size bins,

$$\begin{array}{c}\hfill ⟨\pi s_{\mathrm{p}}^2{n}_{\mathrm{p}}⟩=\pi {\displaystyle \sum _i}{s}_{\mathrm{p},i}^2{n}_{\mathrm{p},i},\end{array}$$(B.3)

which is used to compute the vapor recondensation rate ℛ_cond. The number density of particles in each size bin is given as n_p, i = ρ_p, i/m_p, i, where $m_{\mathrm{p},i}=4\pi {\rho }_{\bullet }{s}_{\mathrm{p},i}^3/3$ is the particle mass. We adopt ${\rho }_{\bullet }=1.5{\mathrm{g}\mathrm{cm}}^{-3}$ as the internal density of pebbles, assuming a compact ice-silicate mixture with half ice in mass (e.g., Schoonenberg & Ormel 2017).

Appendix C: Thermal relaxation timescale

Following Flock et al. (2017), we combine the optically thin and optically thick limits to determine the thermal relaxation time:

$$\begin{array}{c}\hfill t_{\mathrm{relax}}=t_{\mathrm{thin}}+t_{\mathrm{thick}}=\frac{l_{\mathrm{mfp}}^2}{3D_{\mathrm{rad}}}+\frac{H_{\mathrm{phot}}^2}{D_{\mathrm{rad}}},\end{array}$$(C.1)

where l_mfp = 1/(κ_Rρ_g) is the mean free path of photons, $D_{\mathrm{rad}}=16{\sigma }_{\mathrm{b}}{T}^3/(3{\kappa }_{\mathrm{R}}{\rho }_{\mathrm{g}}^2{C}_{\mathrm{v}})$ is the radiation diffusion coefficient and H_phot is the height of disk photosphere (τ_R = 2/3), as measured from the simulation. Here C_v = 8.2 × 10⁷ erg g⁻¹ K⁻¹ is the specific heat capacity of the He-H₂ mixture (assuming 71% H₂ and ideal gas) and κ_R is the Rosseland mean opacity per gas mass, which is a free parameter.

Appendix D: 2D structure of retreating snowline

To illustrate the evolution of disk’s thermal structure during the outburst dimming phase, we present 2D snapshots of our best-fit simulation results at different times since L_bol declines in Fig. D.1. We focus on three epochs: t′ = 110 yr (the observational epoch) and two additional snapshots taken 50 yr before (when recondensation just begins) and after (when most vapor has recondensed in the outer disk) this point (see Fig. 4). At each time, the disk temperature, vapor condensation timescale (t_cond) and number density of water vapor (n_H2O) are shown. The condensation timescale is defined as t_cond = ρ_vap/ℛ_cond.

Fig. D.1. Disk temperature, vapor condensation timescale (tcond), and number density of water (nH2O) during the outburst dimming phase of the best-fit model. Here t′=t − tbeg is the time lapse since the bolometric luminosity started to decline. The value of Lbol corresponding to this time is labeled. In the top row, the range of 140-150 K is shaded in white to highlight the snowline region. At t′ = 121 yr – the present epoch – the emission surfaces (τmm = 2/3) of ALMA band 5, 6 and 7 and disk photosphere (see Sect. C) are plotted along with the temperature contour. In the middle row, the region where no condensation occurs is plotted in white. In the bottom row, two lines denoting ice(vapor)-to-gas ratio of 10−3 are indicated to highlight the snowline region.

As L_bol decreases, the temperature drops in the entire disk, with the midplane cooling slower due to the large infrared optical depth. This difference in cooling time renders the midplane hotter than the upper layer, which is also found in Laznevoi et al. (2025). Concurrently, vapor recondenses and the snowline retreats, as shown in the bottom row of Fig. D.1. The middle row of Fig. D.1 shows the condensation timescale. In the upper layer of the disk, t_cond is significantly longer than near the midplane, leaving uncondensed vapor plumes in the disk atmosphere (bottom row). This explains the slowly declining tail in the time evolution of the water column density (Fig. 4). This “inertia” in vapor’s response to changing luminosity was recently observed in DQ Tau, where the cold water abundance shows little reaction to variations in accretion luminosity on timescales of days to weeks (Kóspál et al. 2025). Radially, t_cond reaches a minimum value at ∼80 − 100 au, increasing outward due to lower pebble densities and inward due to higher temperature. Notably, this increase in condensation timescale was also reported by Smith et al. (2025) in their modeling of EX Lup, implied by the higher water abundance toward the outer disk (their Fig. 10). The ice freezing-out timescale in EX Lup, constrained by multi-epoch observations of Spitzer and JWST (Smith et al. 2025), is one magnitude shorter (∼10 yr) than calculated from our simulations. This difference is expected, given the much larger pebble number density in EX Lup’s inner disk.

Finally, we plot the emission surfaces in our 2sRT calculation (upper row, middle panel), where the millimeter optical depth (τ_mm) equals 2/3. There is a clear elevation of the emission surface across snowline at ≈60 au, rising from the recondensation of vapor. As the emission surfaces trace different disk layers thus different temperature, we expect that the continuum emission will not only depend on the surface density distribution but also be modulated by the disk temperature structure. For molecular emission of water at band 5, because τ_mm ≈ 0.7 at the midplane of 80 au, we only expect minor dust attenuation on measuring the mean water column density at 80-120 au.

Appendix E: MCMC analysis

To fit the ALMA observations, we adopt a log-normal likelihood for each dataset:

$$\begin{array}{c}\hfill ln{\mathcal{P}}_n^i(x)=-\frac{1}{2}{\left(\frac{ln{x}_n^i-ln{x}_{n,\mathrm{syn}}^i}{{\sigma }_{ln,n}^i}\right)}^2-\frac{1}{2}ln(2\pi {{\sigma }_{ln,n}^i}^2),\end{array}$$(E.1)

where $x_n^i$ is either the azimuthally averaged continuum intensity of the n-th annulus in band i (I) or the mean water column density between 80-120 au ($N_{{\mathrm{H}}_2^{18}\mathrm{O}}$). Correspondingly, $x_{n,\mathrm{syn}}^i$ is the synthetic data from the simulations. Following Viscardi et al. (2025), ${\sigma }_{ln,n}^i=\sqrt{B_i/A_n}\times {\sigma }_{x,n}^i/x_n^i$ refers to the fractional root mean square (rms) error of the mean within annulus n, where B_i is the beam area of band i and A_n is the area corresponding to annulus n. For the ${\mathrm{H}}_2^{18}\mathrm{O}$ emission, the area spanning 80 to 120 au is treated as a single annulus. Then the joint log-likelihood reads:

$$\begin{array}{c}\hfill ln\mathcal{P}={\displaystyle \sum _{i=\mathrm{B}6,\mathrm{B}7,n}}ln{\mathcal{P}}_n^i(I)+ln{\mathcal{P}}_1^{\mathrm{B}5}(N_{{\mathrm{H}}_2^{18}\mathrm{O}}).\end{array}$$(E.2)

With this joint likelihood, MCMC analysis was performed for both the retreating snowline model and the static snowline model. In the static model, the parameters t′ and f_cond are replaced by a single L_bol, which controls the temperature structure together with the opacity κ_R in the IR band (Sect. 2.3). Gaussian priors are used for t′ (σ_t′ = 10 yr) and L_bol (σ_Lbol = 110 L_⊙), assuming that the uncertainties in the luminosity curve span 2σ range (i.e., the luminosity curve spans 40 yr). Other parameters adopt uniform priors. For Z_peb, St₀, f_cond and κ_R, sampling is carried out in logarithmic space.

We ran MCMC until all parameters converged over > 50 correlation timescales. The resulting posterior distributions are shown in Figs. E.1 and E.2, with values summarized in Table 1. In the retreating snowline model, several correlations can be seen among disk parameters (Z_peb, γ, γ_peb and St₀) and the gas density slope γ remains poorly constrained. These degeneracies are expected since the adopted disk model parameters do not directly correspond to the observational diagnostics. Specifically, the following features are noted:

Fig. E.1. Corner plot of retreating snowline model. The dashed lines indicate 16th, 50th, and 84th percentiles of the posterior distribution.

Fig. E.2. Corner plot of static snowline model.

1. The parameter St₀, which determines the maximum grain size, is positively correlated with Z_peb and negatively correlated with γ and γ_peb. First, the Stokes number of pebbles decreases in denser gas. A shallower (larger) γ, which implies a higher gas density at 60 au, thus requires a smaller St₀ to maintain the same s_max of ∼ cm, as suggested by the spectral indices (see Sect. 3.1). Second, pebble size affects the continuum emission differently across the disk. In the inner disk where s_max ≫ mm, increasing the pebble size significantly reduces the millimeter opacity. In the outer disk where s_max ≈ mm, however, opacity becomes less sensitive to s_max (e.g., Birnstiel et al. 2018). Consequently, increasing St₀ both reduces the overall continuum intensity and flattens its radial slope, necessitating larger Z_peb and a steeper surface density profile (lower γ_peb) to fit the data.

2. The gas density slope γ is negatively correlated with Rosseland mean opacity κ_R. This arises since to reproduce the intensity shoulder, the snowline should retreat to ≈60 au at the expected observation time (t′), following the assumed luminosity curve. Because the gas density and κ_R together govern the disk’s cooling time (t_relax), that is, the retreat speed, γ and κ_R are naturally related. Furthermore, the posterior of κ_R essentially follows the assumed Gaussian prior of t′.

3. Large γ (toward the upper bound) is preferred. This can be understood as follows: during the outburst dimming phase, the outer disk cools faster than the inner regions, steepening the temperature gradient over time (Laznevoi et al. 2025). A higher γ (shallower density profile) therefore tends to equalize the cooling rate across the disk, "tilting" the intensity profile toward a plateau, as being observed. A value γ > −1 may fit the data even better, but such flat profiles are implausible.

4. The factor f_cond, which controls the condensation rate, is negatively correlated with κ_R. An increase in κ_R slows the retreat of the snowline (larger t′), allowing a reduced f_cond to achieve similar condensation levels.

From Fig. 3, we find that the model underestimates the emission within 40 au and correspondingly, struggles to capture the declining spectral index toward inner disk. These discrepancies likely arise from intense viscous heating (Alarcón et al. 2024) and optically thick environment (Tobin et al. 2023) in V883 Ori’s inner disk., which is included in the model. Simple power-law profiles, as adopted in this work, cannot simultaneously reproduce both the inner and outer disk, implying a sharp transition in the disk properties of V883 Ori.

In the static snowline model, many parameters are driven toward the edges of their prior ranges, indicating that the model cannot adequately reproduce the data. In particular, the disk struggles to reach sufficiently high temperatures. Consequently, κ_R and γ are pushed to their minimum value, and L_bol is constrained to the upper wing of its Gaussian prior, trying to obtain maximum penetration of stellar photons into the disk midplane regions.

All Tables

All Figures

	Fig. 1. Evolution of the bolometric luminosity of V883 Ori during the outburst dimming phase. The symbols denote measurements from different studies, with horizontal bars showing the time range of adopted data. The blue band represents the adopted luminosity curve, which carries an uncertainty of ∼40 yr due to the uncertainties in estimating Lbol.
In the text

	Fig. 2. Time evolution of the surface density profiles of the best-fit model. Silicate surface density remains constant with time. The thick line denotes the best fit, t′ = 121 yr.
In the text

Fig. 3. Comparison of the synthetic intensity profiles with ALMA continuum. Crosses show ALMA data adopted from Houge et al. (2024), spaced according to the beam sizes; shaded regions indicate the rms error. Solid lines represent the best-fit retreating snowline model, while dotted lines correspond to the best-fit static model. Enhanced emission toward the inner regions arises likely from intense viscous heating (Alarcón et al. 2024), which is not accounted for in the model.

In the text

Fig. 4. Average column density of H${}_2^{18}$O at 80–120 au under different condensation rates (fcond). The vertical gray band indicates the time range when observing the intensity shoulder. The dashed black line denotes the average column density derived from ALMA band 5 data by Tobin et al. (2023) with 1σ uncertainty shaded in gray. Dotted gray line denotes the best-fit value in static snowline model.

In the text

Fig. D.1. Disk temperature, vapor condensation timescale (tcond), and number density of water (nH2O) during the outburst dimming phase of the best-fit model. Here t′=t − tbeg is the time lapse since the bolometric luminosity started to decline. The value of Lbol corresponding to this time is labeled. In the top row, the range of 140-150 K is shaded in white to highlight the snowline region. At t′ = 121 yr – the present epoch – the emission surfaces (τmm = 2/3) of ALMA band 5, 6 and 7 and disk photosphere (see Sect. C) are plotted along with the temperature contour. In the middle row, the region where no condensation occurs is plotted in white. In the bottom row, two lines denoting ice(vapor)-to-gas ratio of 10−3 are indicated to highlight the snowline region.

In the text

	Fig. E.1. Corner plot of retreating snowline model. The dashed lines indicate 16th, 50th, and 84th percentiles of the posterior distribution.
In the text

	Fig. E.2. Corner plot of static snowline model.
In the text