Surviving the Heat: multi-wavelength analysis of V883 Ori reveals that dust aggregates survive the sublimation of their ice mantles

Investigating the response of icy dust aggregates to water ice sublimation is essential for understanding the formation and properties of planetesimals in protoplanetary discs. However, their fate remains unclear, as previous studies suggest aggregates could either survive or completely fall apart to (sub){\mu}m-sized grains. Protoplanetary discs around stars undergoing accretion outbursts represent a unique laboratory to study the ice sublimation process, as the water snowline is pushed outward to regions accessible to current observatories. In this work, we aim to understand the aggregates' response to ice sublimation by focusing on V883 Ori, a system currently undergoing a powerful accretion outburst. We present new analysis of archival high resolution ALMA observations of the protoplanetary disc of V883 Ori at 0.88, 1.3, 2.0, and 3.1 mm, and derive new radial spectral index profiles, which we compare with predictions from one-dimensional dust evolution simulations. In the region of V883 Ori where water ice has sublimated, we find lower spectral indices than previously obtained, indicating the presence of cm-sized particles. Coupled with our dust evolution models, we find that the only way to explain their presence is to assume they formed before the outburst, and survived the sublimation process. The resilience of dust aggregates to such intense events leads us to speculate that it may extend to other environments with more gentle heating, such as pebbles drifting through the water snowline in quiescent protoplanetary discs. In that case, it may alter the formation pathway of dry planetesimals interior to the snowline.


INTRODUCTION
Planets form in discs of dust and gas surrounding young stars.Inherited from the collapsed molecular cloud, the initial reservoir of solids is made of m-sized dust grains, which coagulate into larger and larger aggregates through sticky collisions at low velocities (Weidenschilling & Cuzzi 1993).When particles reach ≈ cm-size, socalled pebbles, the growth is halted for example due to fragmenting collisions or efficient inward radial drift (Blum & Münch 1993;Weidenschilling 1997;Blum & Wurm 2008).The streaming instability (SI) is a promising mechanism to overcome the growth barrier and achieve the formation of km-sized planetesimals, by forming highly concentrated clumps of pebbles that become gravitationally unstable (Youdin & Goodman 2005;Johansen et al. 2007Johansen et al. , 2009;;Bai & Stone 2010).But to be triggered, the SI requires regions of protoplanetary discs with a pre-existing enhanced dust-to-gas ratio, as compared to the typical value of 1 per cent.(Johansen et al. 2009;Bai & Stone 2010;Carrera et al. 2015;Li & Youdin 2021).
In that context, the water snowline is thought to play an important role, as planetesimals have been argued to form preferentially (a) interior to the snowline, because of the reduced stickiness of dry silicates as compared to ice-rich grains which leads to a traffic jam ★ E-mail: ah1006@exeter.ac.uk effect (Dominik & Tielens 1997;Supulver et al. 1997;Wada et al. 2013;Gundlach & Blum 2014), and (b) just exterior to the snowline because of the cold finger effect increasing the ice surface density (Stevenson & Lunine 1988;Cuzzi & Zahnle 2004;Ros & Johansen 2013;Drążkowska & Alibert 2017;Schoonenberg & Ormel 2017;Lichtenberg et al. 2021).Formation route (a) has been recently called into question, with the traditional view of a lower  frag for silicates as compared to ice-rich particles being challenged by new laboratory experiments (e.g.Gundlach et al. 2018;Musiolik & Wurm 2019).If, however, icy pebbles happen to disintegrate into m-sized grains upon drifting through the snowline (many-seeds scenario, see e.g.Schoonenberg & Ormel 2017), their drift velocity would still be reduced, offering another pathway for the traffic jam effect and the formation of dry planetesimals.Interestingly, the release of m-sized grains would also boost formation route (b), as small particles can diffuse back across snowline and enhance further the solids-to-gas ratio (Saito & Sirono 2011;Ida et al. 2016;Schoonenberg & Ormel 2017;Hyodo et al. 2019Hyodo et al. , 2021)).In summary, the response of icy dust aggregates to sublimation plays an important role in dust evolution and planetesimal formation around the snowline.However, the water snowline is generally located relatively close to its host star (e.g.Mulders et al. 2015), in regions inaccessible to current observatories at millimetre wavelengths, which prevents us to probe the state of the dust size distribution.
A unique opportunity to study how the water snowline and ice sublimation affect the dust population is offered by protoplanetary discs undergoing FUor-type accretion outbursts.During such events, the luminosity and mass accretion rate of the central protostar increase by ∼2−3 orders of magnitude, and can remain elevated for hundreds of years (Kenyon et al. 1990;Dunham & Vorobyov 2012;Audard et al. 2014;Fischer et al. 2023).As a result, the dust and gas content in the surrounding protoplanetary disc significantly warms up, leading for example to the crystallization of amorphous silicates in the inner disc (Ábrahám et al. 2009), and sublimation of ice species on large scales as snowlines effectively move outwards (e.g.Jørgensen et al. 2013;Molyarova et al. 2018;Wiebe et al. 2019).
A few protoplanetary discs are known to be currently undergoing FUor-type accretion outbursts (Fischer et al. 2023).Notably, V883 Ori is a young 1.3 M ⊙ star surrounded by a 0.2 − 0.6 M ⊙ protoplanetary disc still embedded in its envelope (Cieza et al. 2016(Cieza et al. , 2018;;Leemker et al. 2021), undergoing a powerful FUor-type accretion outburst that likely began before 1888 (Pickering 1890).The system received a lot of attention with ALMA in the last few years, including a suite of spectral line observations (e.g.van't Hoff et al. 2018;Lee et al. 2019;Leemker et al. 2021;Ruíz-Rodríguez et al. 2022;Tobin et al. 2023).A first analysis of the water snowline of V883 Ori was made by Cieza et al. (2016), using dust continuum emissions exclusively in ALMA Band 6.The continuum data revealed a dark annulus at 42 au as well as an abrupt change in optical depth and millimetre spectral index.These features matched prediction of steady state (i.e.non-outbursting) dust and snowline models (Banzatti et al. 2015), prompting Cieza et al. to conclude the snowline in V883 Ori had been moved significantly outwards by the outburst.More recently, Tobin et al. (2023) estimated the water snowline to be located at 80 au based on direct measurements of the column density of HDO and H 2 18 O.
Studies of the dust size distribution in discs undergoing accretion outbursts suggests different scenarios for the response of icy aggregates to the sublimation of their ice content.Based on ALMA Band 6 observations of V883 Ori presented in Cieza et al. (2016), Schoonenberg et al. (2017) developed a monodisperse dust coagulation and sublimation model and obtained close agreement to the measured spectral index with a model in which: (1) icy aggregates are completely disrupted following sublimation at the outburst onset, and (2) their m silicate seeds rapidly re-coagulate to the local fragmentation threshold of 300 m.Conversely, Liu et al. (2021) employed multi-wavelength observations to infer the presence of large pebbles (≳ 1.6 mm) in the inner ∼ 10 au, 400 K hot disc around FU Ori, also currently undergoing an accretion outburst, indicating that here pebbles had survived the sublimation process.An alternative way of studying the question of the dust response to ice sublimation comes from laboratory experiments, in which mixtures of dust/ice can be heated and studied in controlled settings.Previous laboratory studies also found a variety of outcomes, including complete disintegration of the dust aggregates back to m silicate seeds (Aumatell & Wurm 2011), or survival of the refractory matrix (Spadaccia et al. 2022).
In this paper, we aim to shed some light on how icy dust aggregates respond to water ice sublimation by performing a new multiwavelength analysis of V883 Ori, using archival ALMA data (see Sect. 2) and the latest insights concerning the location of the water snowline (Tobin et al. 2023).We select archival ALMA data at 0.88 mm (Band 7), 1.3 mm (Band 6), 2.0 mm (Band 4) and 3.1 mm (Band 3).Given that the snowline is generally located at ≈ 3 au in quiescent conditions, we focus on characterising the dust size distribution of V883 Ori in between 3 <  < 80 au, where the dust population underwent sublimation at the onset of the outburst.From the multi-wavelength dataset, we robustly retrieve the radial profile of the dust emission spectral index (see Sect. 3), a key tracer for dust particle size1 .We additionally perform 1D dust coagulation and evolution simulations using the state-of-the-art DustPy package (Stammler & Birnstiel 2022), adapted to include the effects of an accretion outburst.Given the variety of outcomes uncovered by laboratory studies (Aumatell & Wurm 2011;Spadaccia et al. 2022), we incorporate two end-member scenarios for the response of icy particles to water ice sublimation: the many-seeds scenario, in which all particles fall apart to m-sized grains after sublimation; and the resilient scenario, in which dust particles survive (see Sect. 4).In Section 5, we present the spectral index radial profiles from the dust evolution models in the two scenarios, and compare them qualitatively with the observational results.We discuss our main findings and their implications in Section 6, followed by our conclusions in Section 7.

OBSERVATIONS AND DATA REDUCTION
To characterise the dust size distribution in V883 Ori, it is crucial to perform a multi-wavelength analysis with high resolution data.We gathered from the ALMA archive the highest resolution observations available in Band 7 (project code: 2016.1.00728.S, PI: Lucas Cieza), Band 6 (project code: 2015.1.00350.S, PI: Lucas Cieza), Band 4 (project code: 2017.1.00388.S, PI: Hauyu Baobab Liu), and Band 3 (project code: 2016.1.00728.S, PI: Lucas Cieza).In particular, the Band 6 data correspond to those previously used by Cieza et al. (2016) to obtain the first estimate of the spectral index in V883 Ori, and discover the dark annulus at 42 au.A detailed analysis of the Band 7 and 3 data will be presented in Cieza et al. (in prep.).As the high resolution data in Band 7 and Band 6 did not well cover the shorter baselines, we additionally retrieved lower resolution archival data in these two bands (resp.project code: 2015.1.00041.S, PI: John Tobin, and project code: 2013.1.00710.S, PI: Lucas Cieza).Additional details about the observations can be found in Table 1.We note that while the observing dates differ by up to 2 years, we do not expect the thermal dust emission of FU Orionis objects like V883 Ori to vary on such short timescales (e.g.Wendeborn et al. 2020).
The data were initially calibrated by ALMA staff, some manually and some using the ALMA pipeline.We used CASA (Common Astronomy Software Applications v6.4.1;McMullin et al. 2007) for further inspection, self-calibration, and imaging.As a first step, we aligned the data by fitting a Gaussian to the image of each individual execution and shifting the phase center of each measurement set.Mismatches in the amplitude scales of different observations may also appear.We therefore rescaled the fluxes of each band following the methods in Andrews et al. (2018), considering one of the executions at each band as reference and scaling the others by a factor using the task gencal.The Band 7 observations have slightly different spectral setups between the long and short baseline data (see Table 1), so in this band we only scaled the fluxes between the three short baseline executions.
We then performed self-calibration on each individual dataset; except on the short baselines of Band 7, where we started from the concatenation of the three executions.For the self-calibration, we employed the task tclean with the hogbom deconvolver, a briggs weighting with  = 0.5, and  = 1, i.e. assuming a linear dependency for the emission frequency.We started with phase self-calibration following an iterative approach, decreasing step by step the solution intervals until it no longer improved the signalto-noise ratio (S/N).This was followed by a round of amplitude self-calibration ( =  ), which did not provide any improvements of the S/N and was therefore not considered.Overall, self-calibration on individual datasets provided improvements of the S/N ranging in between 10 and 220%.We then combined the short and long baseline executions of Band 7 and 6.
After combining the short and long baseline executions, we performed an additional cycle of phase self-calibration to ensure a good alignment of the data.This was followed by a round of amplitude self-calibration, especially important for the Band 7 dataset whose fluxes were not rescaled due to the frequency offset in the spectral windows.We used the task tclean now using  = 2, i.e. assuming the frequency dependency follows a Taylor expansion to quadratic terms.This second round of self-calibration provided an improvement of the S/N in Band 7 and 6 respectively by 13 and 32%.
We then built the final images with the task tclean, using the mtmfs deconvolver (Rau & Cornwell 2011),  = 2, multiple scales from the point-source to 5 times the beam size, and a Briggs weighting with different  values to explore higher resolutions.In addition to the continuum images, we derived the dust emission spectral index  =  ln   / ln , which is a widely used metric to trace the size distribution of dust particles in protoplanetary discs (e.g.see reviews by Testi et al. 2014;Andrews 2020).To do so, we built spectral index maps by combining pairs of bands in the task tclean, which derives the spectral index at each pixel based on the observed visibilities, again using the mtmfs deconvolver and  = 2.This allowed us to obtain three spectral index maps from the four bands, while also benefiting from a better  coverage as compared to deriving the spectral index from the intensity measured in individual bands (see also Huang et al. 2018).More details about the final images can be found in Table 2.

Final ALMA images
We present in Fig. 1 our final images of the dust thermal emission of the protoplanetary disc around V883 Ori at 0.88 mm (Band 7), 1.3 mm (Band 6), 2.0 mm (Band 4) and 3.1 mm (Band 3).In Band 6, our image is obtained from the same dataset presented in Cieza et al. (2016), although here combined with another execution to better cover the shorter baselines.In all bands, the dust emission appears azimuthally symmetric.A ring-shaped structure appearing as a dark annulus is well visible at 42 au in Band 7 and Band 6.This feature was also seen in Band 6 by Cieza et al. (2016).The dark annulus is also slightly visible in Band 4, but absent in Band 3, which we explain by the lower resolution of this observation (Table 2).

Spectral index radial profiles
We present in Fig. 2 the radial profiles of the dust emission spectral index at 0.88 − 1.3 mm (red curve), 1.3 − 2.0 mm (orange curve), and 2.0 − 3.1 mm (green curve).They are obtained from the spectral index maps mentioned in Sect.2, where we average the spectral index in elliptical rings of width 5 mas with an eccentricity given by the disc inclination  = 38.3• and position angle P.A. = 32.4• (Cieza et al. 2016).At each radius, the spectral index is given by the mean ᾱ in the ring, and the uncertainty is calculated as the error of the mean following (Macías et al. 2021) where   is the standard deviation in the ring,   is the area of the ring, and  beam is the area of the beam.We additionally show the azimuthal average of the spectral index error maps produced by tclean in fainter shades.In the lower right of each panel, we illustrate as solid grey lines the absolute flux calibration systematic errors of ALMA, assuming the nominal errors of 10% in Band 7, and 5% in Band 6, 4, and 3.Such systematic uncertainty could increase or decrease the spectral index at all radii up to a maximum shift indicated by the grey lines, but it would not affect the relative shape of the individual profiles.Our multi-wavelength dataset also suggests that these systematic errors should not be significant, as an important error in one band would lead to shifts in the other profiles, which is not visible in Fig. 2.
The three spectral index profiles display common features.The inner disc ( < 20 au) is characterised by a region with  ≲ 2, distinctive of optically thick emissions with high dust albedo, in which case the emission is fainter than the black-body radiation (e.g.Birnstiel et al. 2018).It is followed by a transition (20 <  < 50 au) into the optically thin outer disc ( > 50 au) characterised by  ≈ 2.5.In this region, we notice that the dark annulus visible in the final images around 42 au (Fig. 1) has no strong imprints in the spectral index profiles.Then, even though higher than in the optically thick inner part, the values at  > 50 au are small, especially in comparison to what was obtained by Cieza et al. (2016) (see solid grey line in the left panel of Fig. 2).In fact, they found values of about  ≈ 3.5 at radii  > 50 au.The difference can be explained by the fact that they only used data from Band 6 (see Sect. 2), while we are computing our spectral indices with multiples bands, allowing us to probe a larger range of wavelengths and leading to a more robust estimate of the dust spectral index.This difference has important implications for the inferred particle sizes.Schoonenberg et al. (2017) best reproduced this larger spectral index (grey line in the left panel of Fig. 2) with a model featuring the many-seeds scenario and the rapid re-coagulation of the grains into 300 m aggregates.Instead, the lower value of the spectral index obtained here ( ≈ 2.5), as we will see further, implies the presence of considerably larger particles.
Individual profiles also display distinct features.The radial extent of the inner region characterised by  ≲ 2 varies with the wavelength, extending up to 30 au in the spectral index profile at 0.88 − 1.3 mm, while reaching 10−20 au at 1.3−2.0mm and 2.0−3.1 mm.Due to the steep slope that follows the inner  ≲ 2 region, its exact radial extent may be affected by the absolute flux calibration systematic errors of ALMA.The value of the spectral index in the outer disc also differs, being slightly higher ( ≈ 2.75) at 1.3 − 2.0 and 2.0 − 3.1 mm.These wavelength-specific features can be used to characterise the properties of dust particles, such as the dust temperature or surface density, which will be the focus of future work.For  > 125 au, the spectral index at 1.3 − 2.0 and 2.0 − 3.1 mm seems to decrease with the distance to the star.This could indicate the presence of larger aggregates, however, at such large distances, we argue that it might just be due to the larger errors.
Overall, our spectral index radial profiles suggest the presence of larger particles than previously assumed in V883 Ori ( ≈ 300 m, Schoonenberg et al. 2017).In fact, assuming an MRN size distribution, it would require maximum particle sizes in excess of 1 cm to reach spectral indices as low as  = 2.5 in the optically thin regime (e.g.Fig. 4 in Testi et al. 2014).However, fragmentation, radial trans- port, and the accretion outburst itself may significantly alter the size distribution from a standard MRN distribution (Brauer et al. 2008;Birnstiel et al. 2011;Houge & Krĳt 2023).To further understand what our observational results imply concerning the physics of dust aggregates in outbursting discs, we will dedicate the next part of this manuscript to study two specific scenarios (many-seeds vs. resilient) with dust evolution simulations that solve self-consistently for the evolution of the entire dust size distribution.Note that obtaining a detailed fit of the observed spectral indices is beyond the scope of this paper, given the intricacies of all the different dust evolution models, and unknowns about disc parameters and initial conditions.Rather, we seek to model how the two end-member cases (many seeds vs. resilient) impact the dust distribution and evolution of the spectral index, to distinguish which scenario is the most appropriate to describe the dust content of V883 Ori.Specifically, our main goal is to identify which model can reproduce the newly-inferred low spectral index ( ≈ 2.5) in the optically thin region ( > 50 au).

MODELING THE DUST EVOLUTION IN V883 ORI
To explore the evolution of the dust population in V883 Ori, we use the code DustPy 2 (Stammler & Birnstiel 2022) which models the coagulation, fragmentation, and transport of dust particles in a 1D radial disc profile.In this Section, we describe our disc and dust models, and how we implemented an accretion outburst and the two 2 We used DustPy v1.0.3 for the simulations presented in this manuscript.
sublimation scenarios into DustPy.We present the results of our dust evolution simulations in Sect.4.5.

V883 Ori disc model
We modify the initial disc model of DustPy based on observational constraints of V883 Ori.V883 Ori is a T Tauri star of mass 1.3 M ⊙ (Cieza et al. 2016) surrounded by a disc of estimated mass  disc ≈ 0.3 M ⊙ (Cieza et al. 2018;Leemker et al. 2021).We assume the initial gas surface density profile to follow a tapered power-law (Lynden-Bell & Pringle 1974;Hartmann et al. 1998) where  c = 75 au (van't Hoff et al. 2018;Leemker et al. 2021),  = 1.25, and We simulate the disc over a discretized grid of 100 radial cells logspaced between 1 au and 200 au.We initialize the simulation during the past quiescent phase of V883 Ori, where we assume the disc to be passively irradiated, characterised by a midplane temperature profile where  = 0.05 is the disc opening angle and  * = 6 L ⊙ is the stellar luminosity (Cieza et al. 2016).Assuming the sublimation The white curve indicates the dark annulus at 42 au.The color bar in each panel goes from about twice the rms up to the maximum intensity.More information about the final   values used, the rms, peak S/N, and beam shape can be found in Table 2.

Dust model
We initialize the simulation with a population of grains following a size distribution characteristic of the ISM (Mathis et al. 1977), with a power-law () ∝  −3.5 from 0.5 m monomers up to 100 m aggregates (assuming dust growth in the early stages of the disc evolution, e.g.Bate 2022).We adopt a dust composition following Birnstiel et al. (2018), where grains are made of a mixture of silicates, troilite, and refractory organics, with a relative mass abundance respectively set to 0.411, 0.093, and 0.496.With a respective material density of 3.3, 4.83, and 1.50 g cm −3 , the overall bulk density of dust grains is  s = 2.11 g cm −3 .Outside the water snowline, H 2 O molecules are frozen out onto dust particles such that the water fraction is 20% of the total mass, comparable to measurements of comet 67P/Churyumov-Gerasimenko (Pätzold et al. 2016).Using a water ice density of 0.92 g cm −3 , it leads to icy dust grains with a bulk density  s = 1.67 g cm −3 .We set the fragmentation velocity to  frag = 1 m s −1 inside the water snowline, and 10 m s −1 outside, in agreement with laboratory experiments finding enhanced stickiness for water-rich dust particles (e.g.Supulver et al. 1997;Gundlach & Blum 2014).Regarding the dust structure, we assume dust grains to coagulate following a compact aggregation model, remaining spherical and homogeneously mixed in dust and ice at all times.

Outburst event
We let the dust population evolve in quiescent conditions until  otb = 10 5 yr, marking the onset of the accretion outburst.This time, chosen to be long enough for aggregates to reach the fragmentation barrier The shaded areas are the error of the mean at each radii, and in fainter colours the azimuthal average of the spectral index error maps produced by tclean.The bottom right inset shows the geometric mean of the beam size for each spectral index map (Table 2) and the vertical grey lines represent the systematic error introduced by the flux calibration uncertainty.Horizontal dotted lines indicate  = 2 (optically thick emission in the Rayleigh-Jeans limit) and  = 3.7 (typical value for micrometer dust grains in the ISM, see Finkbeiner et al. 1999).The vertical dashed lines represent respectively the dark annulus visible in Fig. 1 and the water snowline position (Tobin et al. 2023).
but short enough for radial drift to not have removed a considerable amount of dust (see however Sect.5.1), is comparable to the typical time span between outburst events (Contreras Peña et al. 2019), and, within a few factors, in agreement with the estimated age of V883 Ori (Cieza et al. 2016).After that time, we mimic an accretion outburst by raising the dust and gas temperature.We assume the outburst temperature profile to be constant throughout the accretion outburst, and rather than using Eq. 4 with the current bolometric luminosity of V883 Ori, we choose to fix the new outburst temperature profile to where the power-law index is identical to the one used in Schoonenberg et al. (2017) and consistent with the one of a passively irradiated disc (Kenyon & Hartmann 1987); and the normalisation has been updated so that the water snowline position corresponds roughly to 80 au, as suggested by 3 Tobin et al. (2023).We discuss the assumption of a constant temperature profile in the context of outburst decay and variability in Sect.6.5.We assume the temperature to increase instantaneously to the outburst value, as the heating timescale of the dust is generally smaller than the typical outburst duration (Johnstone et al. 2013).We further assume the disc vertical structure to instantaneously adapt to a new hydrostatic equilibrium when the temperature is modified.This is a valid assumption when the gas thermal timescale (e.g.Equation 7in Ueda et al. 2021) is short compared to the dynamical timescale and the outburst duration, which is the case in our disc model at all relevant radii.The onset of the accretion outburst in V883 Ori was not directly observed, but based on observations of its reflection nebula in 1888 by Pickering (1890), Strom & Strom (1993) argued that it has been active for at least 135 yr, and likely began even before that observation.Due to the uncertainty regarding the exact onset time, we keep the outburst active for 10 4 yr.This duration is considerably longer than what is expected for FUor-type accretion outbursts (Audard et al. 2014), but is chosen to allow us to study the evolution of dust particles on longer timescales as they readjust to the new coagulation/fragmentation equilibrium (see Sect. 4.5).We end the simulation after that time ( otb + 10 4 yr), and do not model the post-outburst phase and the recondensation of water vapor (for that matter, see Vorobyov et al. 2022;Houge & Krĳt 2023).

Dust response to sublimation
In the midplane between the initial quiescent snowline (3 au) and outburst water snowline (80 au), water is initially in ice phase deposited onto the surface of refractory particles, before sublimating when the outburst starts at  =  otb .This zone will be further referred to as the sublimation region.Similarly to previous works (Cieza et al. 2016;Schoonenberg et al. 2017), we assume water ice to sublimate instantaneously at the outburst onset.It is a valid assumption in warm outbursting conditions, as the desorption of volatile molecules is very sensitive to the temperature (e.g.Supulver & Lin 2000).Moreover, in V883 Ori, water vapor and signs of an active water chemistry with HCO+ is observed up to 100 au (Leemker et al. 2021;Tobin et al. 2023), providing additional evidence that dust particles up to such radii underwent heating and sublimation.
In the sublimation region, the composition of dust particles is modified at the outburst onset, leading to a change in bulk density  s and in the fragmentation velocity  frag , decreasing from 10 to 1 m s −1 (Sect.4.2).In addition, even in the absence of collisions, the sublimation of water ice may lead directly to the disintegration of icy dust aggregates (see Sect. 1).Laboratory experiments by Aumatell & Wurm (2011) investigated the sublimation of mm-sized pure ice pebbles, and found that they break up in a large number (≈ 104 ) of smaller pieces, supporting the idea that sublimation disrupts icy dust aggregates into many m-sized silicate grains.More recently, however, Spadaccia et al. (2022) studied mm-sized pebbles made of a mix of dust and ice and found that they are able to survive sublimation, especially if the ice content is around 15% of the pebble mass, which is around what we expect for icy particles in protoplanetary discs (see Sect. 4.2).In agreement with their findings, we implement two scenarios for the response of icy dust aggregates to sublimation at the outburst onset: (a) the many-seeds scenario, where we disrupt all dust particles in the sublimation region back to the monomer size (0.5 m); and the resilient scenario, where particles survive and are left unchanged.Note that in the resilient scenario, individual particles masses would also decrease by a fraction relative to their ice content (20% in our dust model, see Sect.4.2), but as shown in Houge & Krĳt (2023), it is the compositional change and subsequent lowered fragmentation velocity that dominates the appearance and behaviour of resilient aggregates.Furthermore, we do not extend the many-seeds model to the general case of pebbles drifting through the water snowline at any time of the simulation.At the snowline location, we therefore keep the standard implementation of DustPy, where pebbles survive as they drift through the snowline, but quickly fragment as they become less resistant to fragmentation.

Many-seeds scenario
We begin with the results from the many-seeds scenario.We present in Fig. 3 the dust density distribution as a function of particle size at different key times of the simulation: right before the outburst (topleft); after an outburst duration of 10 yr (top-right); after an outburst duration of 100 yr, which corresponds approximately to the minimum duration of the outburst in V883 Ori following the measurements of Pickering (1890) (bottom-left); and after an outburst duration of 10 4 yr (bottom-right).As mentioned in Sect.4.3, the latter duration is much longer than what is typically expected for an FUor-type outburst (≈ 100 yr), but is included to allow us to discuss how dust readjust to the new coagulation/fragmentation equilibrium.The blue and green lines in each panel in Fig. 3 represent respectively an estimate of the current fragmentation and radial drift barriers, after which dust growth through coagulation is halted (Birnstiel et al. 2012).
Before the outburst (top-left panel in Fig. 3), the dust population has evolved for 10 5 yr and reached the fragmentation or radial drift barrier up to ≈ 100 au, reaching sizes ≥ 1 cm.When the outburst starts (top-right panel in Fig. 3), water sublimates in the sublimation region (3 <  < 80 au), upon which all dust particles fall apart to m-sized silicate grains, in agreement with the many-seeds scenario (Sect.4.4).Due to the change in composition in the sublimation region, the fragmentation barrier (solid blue line) decreases by about two orders of magnitude, as  max ∝  2 frag (Birnstiel et al. 2011).The fragmentation barrier also slightly decreases throughout the disc due to the rise in temperature, as it results in larger relative velocities between dust particles ( max ∝  −1 , see Birnstiel et al. 2011).This temperature-induced fragmentation can be seen in the inner disc ( < 3 au) on the top-right panel of Fig. 3, where pebbles that are now too large (i.e.above the blue line) fragment to reach the new fragmentation barrier.Thanks to high density in the inner disc and the small amplitude of the size change required, it only takes a few years for the inner disc population to adapt to the new coagulation/fragmentation equilibrium.Note that these particles did not disrupt to m-sized grains because they are located inside the quiescent water snowline, where water is in a vapor state at all times.We can also see a thin patch of large particles (>1 cm) directly at the quiescent snowline, corresponding to pebbles that just drifted through the snowline before the outburst started, and have not yet fragment.The radial drift barrier (solid green line) also decreases by factor of a few, which only impacts particles outside of the outburst snowline ( > 80 au).
After 100 yr past the outburst onset (bottom-left panel in Fig. 3), we notice that dust aggregates are re-coagulating in the sublimation region.However, the collision rates decrease with increasing radial distance, which leads to a growth front propagating outward with time since the outburst onset (see also Vorobyov et al. 2022).After 100 yr, it is located at ≈ 5 au: dust particles inside of that radius already reached the new fragmentation barrier, while particles outside of that radius are still smaller than 10 m.Even after 10 4 yr (bottomright panel in Fig. 3), the dust population has not reached the new fragmentation barrier outside of 20 au.After such extended outburst duration, we also start to notice the impact of radial drift, replenishing the region between 30 − 80 au in mm-sized particles although most of the dust mass in this region is still stored in grains with sizes below 10 m.
Overall, in the many-seeds scenario, m-sized grains require considerable timescales to re-coagulate, and no dust aggregates over a millimetre size can be found in the sublimation region (3 <  < 80 au) shortly after the outburst.The only exception comes from pebbles drifting through the outburst water snowline ( > 80 au), but these only appear after ≈10 4 yr.

Resilient scenario
We now present in Fig. 4 the dust density distribution as a function of particle size in the resilient scenario, at identical times as Fig. 3.The pre-outburst distribution is identical to Fig. 3, with the dust growth halted by the fragmentation or radial drift barrier up to ≈ 100 au.When the outburst starts (top-right panel in Fig. 4), the resilient scenario (Sect.4.4) assumes dust aggregates survive sublimation and the dust distribution is therefore similar to the preoutburst state.However, as the composition changed to dry silicates, the fragmentation barrier has decreased by two orders of magnitude ( max ∝  2 frag , see Birnstiel et al. 2011).Pebbles that are now too large, i.e. above the blue line, start fragmenting, trying to settle in the new (local) coagulation/fragmentation equilibrium.Nevertheless, it requires time to fragment all the pebbles, and after 100 yr (bottomleft panel in Fig. 4), fragmentation still dominates.Similarly to the many-seeds scenario, the dust population responds to the outburst inside-out, adapting faster closer to the host star thanks to larger dust densities and collision rates.It leads to a fragmentation front propagating outward with time since the outburst onset 4 .Even after 10 4 yr (bottom-right panel in Fig. 4), pebbles are still fragmenting and a considerable amount of dust with sizes above the local fragmentation barrier (i.e.above the blue line) persists at radii between ≈ 10 and 80 au.
In summary, at several tens of au, where the collision timescales are fairly long as compared to the outburst typical duration (≈ 100 yr), the evolution of the dust size distribution during an outburst is set  primarily by how individual icy dust aggregates respond to ice sublimation at the outburst onset.

COMPARISONS BETWEEN OBSERVATIONS AND DUST SIMULATIONS
In this section, we convert the results from our dust evolution models into spectral indices, and compare them with the observed values (Fig. 2).Given the dust composition detailed in Sect.4.2, we take the optical constants from Draine (2003) for astronomical silicates, Henning & Stognienko (1996) for troilite and refractory organics, and Warren & Brandt (2008) for water ice.We then compute the mixed dielectric function using the Bruggeman effective medium theory (Bruggeman 1935).Next, we compute the absorption opacity and scattering opacity5 using the DSHARP-OPAC package from Birnstiel et al. (2018).Finally, using the dust density distribution at each radius (self-consistently derived from our dust evolution models), we calculate the total absorption opacity  abs  and total scattering opacity  sca  .From the total absorption and scattering opacity, we can calculate the dust emission intensity   at each radius assuming a 1D slab model following (Sierra et al. 2019;Carrasco-González et al. 2019)  (Sierra et al. 2019;Carrasco-González et al. 2019)  where   = √ 1 −   .This 1D slab model is a good approximation for protoplanetary discs, given that, at sub-mm wavelengths, we mainly trace large pebbles highly settled into the disc midplane.From Eq. 6, we derive the dust emission spectral index as where  1 and  2 are chosen in agreement with the central wavelength of the images at each band (Table 2).Note that the spectral index is computed from our dust evolution models, sampled on the log-spaced radial grid described in Sect.4.1.As a result, the spectral index from our models is not convolved with the resolution of the ALMA images, which would only affect our results at relatively small spatial scales given the high angular resolution of the observations (see Table 2).
We begin with the spectral index radial profiles of the dust evolution model in the many-seeds scenario, presented in the upper panels of Fig. 5 at different times before and after the outburst onset.Following the goal of this work, we focus on the sublimation region of V883 Ori (3 <  < 80 au).The pre-outburst spectral index is relatively close to the observed data, featuring an optically thick inner disc ( < 15 au) with values around  ≈ 2, followed at larger radii by optically thin emissions with a spectral index  ≈ 2.5 at 0.88 − 1.3 mm.The relatively low spectral index is due to the presence of large cm-sized pebbles in these regions during the preoutburst phase (Fig. 3).Comparing different wavelengths, the shape of the pre-outburst profile is similar, although it shifts to slightly higher spectral indices at larger wavelengths due to varying optical depth and albedo.The model does not manage to reproduce the values  < 2 in the inner disc, which could be caused by unknowns about the disc parameters, or by the fact that we do not consider a vertical temperature profile, which was shown to decrease the spectral index in the optically thick regime (Sierra & Lizano 2020).When the outburst starts, dust particles disrupt to m-sized grains, and the spectral indices increase significantly, reaching values close to  = 3.7, similar to small grains in the ISM (Finkbeiner et al. 1999).With time, the spectral index begins to evolve in the innermost and outermost part of the sublimation region.This is due respectively to the re-coagulation of dust particles (see growth front in Sect.4.5.1) and the inward drift of pebbles that survived outside of the outburst water snowline ( > 80 au).However, these effects are not sufficient to replenish in large particles the sublimation region of V883 Ori within the typical range of outburst duration (see also Schoonenberg et al. 2017, for a discussion on drift timescale and outburst duration).
We then present the spectral index radial profiles of the resilient scenario in the lower panels of Fig. 5.The pre-outburst spectral index is identical to the many-seeds case.However, when the outburst starts, the spectral index remains relatively close to the pre-outburst profile, much closer to the observed data (Fig. 2), with values  ≈ 2.5.Over time, more pebbles are fragmenting causing the spectral index to vary and move away from the observed data.This is particularly true in the inner part of the sublimation region (see fragmentation front in Sect.4.5.2).In our dust evolution models, the resilient scenario is therefore the most appropriate to describe the dust population of V883 Ori.

Varying the turbulence parameter
As discussed in Sect.3.2, the spectral index and its temporal evolution are highly sensitive on uncertain disc and dust parameters.Notably, the turbulence parameter  T is rather difficult to estimate in protoplanetary discs (see Rosotti 2023, and references therein), but has a great impact on dust evolution: stronger turbulence leads to higher relative velocities between particles, which shifts the fragmentation barrier to smaller sizes  max , reducing the speed of radial drift (Birnstiel et al. 2011).Moreover, high relative velocities result in higher collision rates, so that coagulation proceeds more rapidly.To investigate the impact of the choice of  T on our results, we perform additional dust evolution simulations using  T = 10 −2 and  T = 10 −4 .We present their spectral index radial profile at 0.88 − 1.3 mm in Fig. 6, and refer the interested reader to Appendix A for more details about how the turbulence influences the dust density distribution.
When comparing the pre-outburst spectral index for different turbulence values (dashed red and dashed blue curves in Fig. 6), we notice the significant impact of the turbulence parameter on the spectral index radial profile.For example, for lower turbulence  T = 10 −4 (right panels in Fig. 6), pebbles are larger and more efficiently transported, preventing them to fragment and redistribute mass in smaller particles.When the outburst starts at  otb = 10 5 yr, the reservoir of mm-sized particles has significantly decreased, leading to a larger spectral index inside 50 au (see Appendix A for more details).However, when comparing the resilient and many-seeds scenarios for identical turbulence values (top vs. bottom panels in Fig. 6), we reach similar conclusions as in the previous section: the spectral index in the many-seeds case is too high, and does not well represent the observed dust population in V883 Ori (Fig. 2).
To conclude, based on comparisons between observations (Fig. 2) and dust evolution models (Fig. 5 and Fig. 6), we find that large dust particles must be present in V883 Ori to explain the observed spectral indices.Given that, in the many-seeds scenario, mm-sized aggregates require timescale significantly larger than the typical outburst duration to re-appear, the only way to explain the dust population in V883 Ori is that the resilient scenario applies, hence that dust particles survived the sublimation of their water ice content at the onset of the accretion outburst.

How does dust respond to an accretion outburst?
Previous works have studied the dust size distribution in detail in two protoplanetary discs undergoing accretion outbursts: FU Ori and V883 Ori.For FU Ori, Liu et al. (2021) employed multi-wavelength observations coupled with MCMC analysis to infer the presence of large pebbles (>1 mm) surviving the water-sublimating temperature, in line with the resilient scenario.For V883 Ori, however, based exclusively on ALMA Band 6 data obtained by Cieza et al. (2016), Schoonenberg et al. (2017) best reproduced the observed spectral index by assuming the fast re-coagulation of refractory dust grains after their total disruption to m size at the outburst onset, i.e. the many-seeds scenario.
In this paper, thanks to new multi-wavelength analysis of V883 Ori in Band 7, 6, 4, and 3 (see Sect. 3), we argue instead that the dust particles in V883 Ori also survived the sublimation event at the outburst onset.With the same scenario fitting FU Ori and V883 Ori, we conclude that icy dust aggregates follow the resilient scenario, and survive to water sublimation during accretion outbursts.As discussed in Houge & Krĳt (2023), these results imply that the impact of accretion outbursts on the dust size distribution is limited, and will mostly concern the composition of pebbles.In fact, the largest pebbles will require considerable time after the outburst to retrieve their pre-outburst water content (between 500 and 2000 yr at 5 au, see Fig. 13 in Houge & Krĳt 2023).These timescales are generally longer than the time required for a cloud of pebbles to gravitationally collapse into planetesimals (∼100 yr, see Visser et al. 2021).Therefore, if the formation of planetesimals is triggered by the outburst, it may lead to the formation of water-poor planetesimals.Note, however, that the collapse of a pebble cloud first requires the formation of a highly concentrated clump of material (Johansen et al. 2007(Johansen et al. , 2009)), e.g.induced by the streaming instability (Youdin & Goodman 2005), which itself takes hundreds of orbital periods to form (Carrera et al. 2015).It should thus be already present for the outburst to have an impact on planetesimal formation.
We note that in this work, we adopt two end-member scenarios: (1) the many-seeds case, where all aggregates fall apart, and (2) the resilient case, where all aggregates survive.However, depending on the complexity of the pre-outburst size distribution, a hybrid scenario can be envisioned.For example, Spadaccia et al. (2022) found that the ability of icy dust aggregates to survive sublimation depends on their water mass fraction, and on the mineralogy and size distribution of the grains they are made of.One could imagine these parameters to vary within the dust population, leading to an hybrid scenario which could be studied if based on additional inputs from laboratory experiments.However, to explain the relatively low values of the spectral indices we measure in V883 Ori (Fig. 2), cm-sized pebbles must still survive sublimation.Whether or not smaller grains follow the same trend would not impact our findings and conclusions.

What about pebbles drifting through the water snowline?
We can compare the behaviour of icy aggregates in the sublimation region of outbursting discs to that of pebbles drifting across the (static) water snowline in a quiescent disc.Accretion outbursts represent an important stress for dust particles, as the temperature changes abruptly to higher values (≳ 100 K in a few years).In comparison, the temperature of pebbles drifting through the water snowline in a quiescent disc increases slowly with the inward radial motion.We can express the heating rate of drifting pebbles as where  r is the radial drift velocity.Considering the upper limit case of pebbles characterised by St = 1, for which the radial drift velocity reaches the maximum value, and using Eq. 4 for the temperature profile in a quiescent disc, we find d/d < 1 K yr −1 at the snowline location (3 au), representing a considerably less abrupt heating than an accretion outburst.For the drift scenario, Piso et al. (2015) found that, as a result of the shallow temperature profile, the desorption timescale of H 2 O ice mantles near the water snowline can even be longer than the drift timescale, suggesting aggregates can traverse a significant distance before shedding their mantles.Based on the  Figure 5. Spectral index radial profiles in the many-seeds scenario (upper panels) and resilient scenario (lower panels) after different outburst duration at 0.88 − 1.3 mm (red curves), 1.3 − 2.0 mm (orange curves), and 2.0 − 3.1 mm (green curves).We include at each wavelength the pre-outburst spectral index as a dashed line, and the observed spectral index with a solid grey line (Fig. 2).
gentler heating and sublimation process, we speculate that if pebbles survive sublimation in an outbursting disc, they should also survive when drifting through the water snowline.The survival of pebbles drifting through the water snowline may have important implications for the formation of dry planetesimals interior to the snowline.In fact, the appearance of the traffic jam effect at that location either requires (Sect.1): (1) that ice-covered particles are more resistant to fragmentation than bare silicates, which has been contested by recent laboratory experiments (e.g.Gundlach et al. 2018;Musiolik & Wurm 2019); or (2) that dust particles fall apart to m-sized monomers, which seems not to be the case as suggested by our findings.As a consequence, if future investigations confirm that ice-rich particles do not benefit from enhanced resistance to fragmentation (e.g. using the approach discussed in Sect.6.4), then the formation pathway of dry planetesimals interior to the water snowline may need to be revised.In addition, without the release of many m-sized grains able to diffuse through the snowline and boost the solids-to-gas ratio, the formation of icy planetesimals exterior to the snowline may also be lowered (Saito & Sirono 2011;Ida et al. 2016;Schoonenberg & Ormel 2017;Hyodo et al. 2019Hyodo et al. , 2021)).

Implications for the internal structure of dust aggregates
The fact that dust aggregates survive the sublimation of their water content can give us constrains about the actual structure of dust particles in protoplanetary discs, and how water ice is distributed on their surface.For icy porous aggregates, it is generally assumed that water ice covers the entire surface of monomers, bonding them to each other (e.g.Okuzumi et al. 2012;Krĳt et al. 2016).However, this approach contradicts our findings.In fact, if water ice glues refractory monomers to each other, then the sublimation of water ice should lead to a disruption of the aggregates, i.e. our manyseeds scenario.Finding large surviving pebbles in V883 Ori and FU Ori teaches us that: either (1) water ice does not cover the entire surface of monomers, which has important consequences for the collisional properties of porous aggregates, by changing the values of key parameters like the fragmentation limit  frag or the rolling energy  roll , or (2) water ice does cover the surface of monomers, but, after sublimation, monomers are able to immediately collide and reconnect the aggregate.The latter could be tested with N-body simulations or laboratory experiments.
Alternatively, our findings concerning the resilient scenario could indicate that icy dust particles are rather made of a large silicate skeleton embedded in a water ice matrix, which is released at the outburst onset, as sublimation should not change the structure of the dust within the matrix.In that case, we speculate that the large dust skeleton would need to form before being covered in ice, meaning that dust particles primarily grow either: 1) in the disc atmosphere, 2) during the warm Class 0 phase (Steinacker et al. 2010;Ormel et al. 2011), or 3) inside the water snowline, before being redistributed to the outer disc.For the latter, we note that there are evidences that efficient outward transport should have occurred in the solar nebula, .Spectral index radial profiles in the many-seeds scenario (upper panels) and resilient scenario (lower panels) after increasing outburst duration at 0.88 − 1.3 mm for  T = 10 −2 (left panels, red curves) and  T = 10 −4 (right panels, blue curves).We include the pre-outburst spectral indices as dashed lines, and the observed spectral index with a solid grey line (Fig. 2). in order to explain the distribution of CAIs in Solar System objects (e.g.Ciesla 2007;Marschall & Morbidelli 2023).

Dust emission substructures in outbursting discs
Our findings have implications for the characterisation of outbursting discs using dust emission alone.First, the resilient scenario complicates the viability of using dust emission to detect the outburst water snowline.In fact, in the many-seeds scenario, we could have observed a sharp transition at the location of the snowline, as, inside, all particles fell apart to m grains, while outside, water ice remained stable and large pebbles survived (Fig. 3).In the favoured resilient case, however, we do not expect such a transition (Fig. 4).In addition, models predicting the signature of the water snowline based on the pile-up of pebbles at that location, such as Banzatti et al. (2015), cannot readily be applied to the snowline of outbursting systems, as the timescale for dust pile-up is orders of magnitude higher than the outburst duration (Schoonenberg et al. 2017).
Second, in the resilient scenario, we discussed in Sect.4.5.2 that the fragmentation of pebbles occurs inside-out from the quiescent snowline, propagating outward with time.As fragmentation destroys particles with sizes  > 1 mm and redistributes their mass in smaller grains, substructures may appear at millimetre wavelengths outside of the quiescent water snowline (see Fig. 5 and Fig. 6).The observation of the fragmentation front could provide key constrains, such as the position of the quiescent snowline or the time passed since the outburst onset (see also Houge & Krĳt 2023).More importantly, it could provide a further understanding concerning the fragmentation velocity  frag of ice-rich and ice-poor dust particles.In fact, there are uncertainties concerning the fragmentation velocity of water-rich particles, with laboratory experiments finding  frag = 10 m s −1 (e.g.Supulver et al. 1997;Gundlach & Blum 2014), while recent works rather find  frag = 1 m s −1 similarly to the value associated with dry silicates (e.g.Gundlach et al. 2018;Musiolik & Wurm 2019).As the fragmentation front is only present if there is a modification of the fragmentation velocity when water ice sublimates, its observation could provide evidence that water-rich particles are indeed more resistant then dry silicates towards fragmentation.Note, however, that it would require high resolution observations, as in our model of V883 Ori the fragmentation front propagates slowly.After 10 000 yr, it is still within ≈ 10 au for  T = 10 −3 (Fig. 4) and ≈ 30 au for  T = 10 −2 (Fig. A1).
Finally, we note that the dark annulus visible in Fig. 1 around 42 au does not have a counterpart in the spectral index profiles (see Fig. 2), suggesting that it is not related to a local change in dust particle size or properties.Instead, the dark feature could be caused for example by a dip in the dust temperature or surface density.In the context of the models presented here, such features might be attributed to (a) shadowing of a disc region by the dense inner disc which contains a large population of small grains due to the fragmentation front (see Fig. 4), or (b) radial drift clearing an annular region of mm-sized particles.Detailed radiative transfer calculations, including evolving dust size distributions and covering the optically thick inner disc, are necessary to test these ideas, and discard its direct association with the water snowline and/or ongoing planet formation.

Accretion luminosity variability
In this paper, following previous modelling efforts (e.g.Schoonenberg et al. 2017;Wiebe et al. 2019;Vorobyov et al. 2022;Houge & Krĳt 2023), we assume a step-like function for the temporal evolution of the disc temperature: the temperature instantaneously rises from the quiescent (Eq.4) to the outburst value (Eq.5), and remains constant throughout the outburst.However, the accretion luminosity and disc temperature during FUor-type accretion outbursts is expected to decay with time passing from the onset (Audard et al. 2014).In V883 Ori, the current bolometric luminosity is ∼ 218L ⊙ (Furlan et al. 2016), while it was measured to be ∼ 400L ⊙ about 20 years ago (Sandell & Weintraub 2001), and was likely higher at its onset, before 1888 (Pickering 1890;Allen et al. 1975;Strom & Strom 1993).In addition to a long-term decay, FU Orionis objects can exhibit low-amplitude variability on short timescales (∼weeks, see Siwak et al. 2018;Takagi et al. 2018), although Wendeborn et al. (2020) found that low-amplitude variability in the millimetre emission of FU Orionis objects is rare and, if present, associated with free-free emission.Extending models like the ones presented here to include a more complex luminosity and gas and dust temperature evolution could lead to interesting behavior, and will be the subject of future works.For example, the heating and cooling rates of dust particles are dependant on the particle size (e.g.Wolf 2003;Gavino et al. 2023), so that the onset of ice sublimation may occur at different times depending on the particle size.Moreover, the freeze-out timescale of water molecules can be long (∼100−1000 yr), which could lead to a decoupling in between the position at which the dust temperature equals the sublimation temperature of water, and the location where water molecules are actually found dominantly in the ice phase (Jørgensen et al. 2013;Visser et al. 2015;Hsieh et al. 2019).

SUMMARY AND OUTLOOK
We present new analysis of archival high resolution ALMA observations at 0.88 mm (Band 7), 1.3 mm (Band 6), 2.0 mm (Band 4) and 3.1 mm (Band 3) of the continuum emission in the protoplanetary disc around V883 Ori, currently undergoing a powerful accretion outburst.From the multi-wavelength data, we obtain new images of the disc (Sect.3.1), and constrain the radial profile of the dust spectral index at 0.88−1.3mm, 1.3−2.0mm, and 2.0−3.1 mm (Sect.3.2).We performed 1D dust evolution simulations, using DustPy, to constrain the state of the dust population in V883 Ori, considering multiple scenarios for the response of dust particles to sublimation of their water ice mantle (Sect. 4.5.1 and Sect. 4.5.2).We then qualitatively compared results of the dust evolution models with the observed spectral indices (Sect.5).Our findings are summarized as follows: 1. We obtain new images of the dust thermal emission in V883 Ori in four bands, and observe a dark annulus around 42 au in Band 4 and Band 7 (Fig 1), as was also seen in Band 6 by Cieza et al. (2016).2. The inner 20 au of V883 Ori are characterised by spectral indices  < 2, indicative of high optical depths dominated by self-scattering.Outside of 50 au, the emission is optically thin with spectral indices around  ≈ 2.5 (Fig. 2).The dark annulus visible in Fig. 1 around 42 au has no strong imprints in the spectral index profiles.3. Our relatively low spectral indices suggest the presence of aggregates larger than 1 cm, unlike what was previously assumed in V883 Ori based on Band 6 measurements (≈ 300 m, see Cieza et al. 2016;Schoonenberg et al. 2017).4. Detailed dust evolution models show that, in the sublimation region (3 <  < 80 au), the evolution of the dust size distribution during an outburst is set primarily by how dust aggregates respond to ice sublimation at the outburst onset.The dust distribution is not in coagulation/fragmentation equilibrium, and particles will either be smaller than 10m (many-seeds scenario, Fig. 3) or larger than 1 cm (resilient scenario, Fig. 4). 5.The resilient scenario is best suited to reproduce the spectral indices observed in V883 Ori ( Fig. 5), as the many-seeds scenario inevitably leads to spectral indices that are too high (i.e. ≈ 3.5), even with additional simulations performed with different turbulence parameters (Fig. 6).As large particles were also found in the outbursting disc of FU Ori (Liu et al. 2021), we conclude that dust particles survive the sublimation of their water ice mantle at the onset of accretion outbursts (Sect.6.1).
Our findings illustrate how constraints on particle size distributions in outbursting systems can be used to study the mechanical and collisional properties of dust aggregates.As FUor-type accretion outbursts are powerful heating events, our findings may have implications for the ice sublimation process in other environments where heating is more gentle.For example, this could apply to icy pebbles accreting within the envelopes of protoplanets (Johansen et al. 2023), or for icy pebbles drifting through the water snowline in quiescent protoplanetary discs, which may have consequences for the formation pathway of dry planetesimals (Sect.6.2).Future studies could lead to a more complete picture, for example by focusing on N-body simulations and laboratory experiments to explore how sublimation may vary for different particle sizes, porosities, and water contents, and assessing the plausibility of the reconnection hypothesis (Sect.6.3).Alternatively, higher resolution or longer wavelengths observations contrasting quiescent and outbursting discs may bring key constraints on the fragmentation velocity of ice-rich vs. ice-poor particles (Sect.6.4).

Figure 1 .
Figure 1.Final images of the dust thermal emission of the disc around V883 Ori at 0.88 (top-left), 1.3 (top-right), 2.0 (bottom-left) and 3.1 mm (bottom-right).The white curve indicates the dark annulus at 42 au.The color bar in each panel goes from about twice the rms up to the maximum intensity.More information about the final   values used, the rms, peak S/N, and beam shape can be found in Table2.

Figure 2 .
Figure2.Spectral index radial profiles of the dust emission of V883 Ori at 0.88 − 1.3 mm (red curve), 1.3 − 2.0 mm (orange curve), and 2.0 − 3.1 mm (green curve).The profiles are retrieved by azimuthally averaging the spectral index maps obtained from the multi-bands images (Sect.2).We also show the spectral index obtained from the Band 6 data derived byCieza et al. (2016) (grey curve).More information, including the error bars, can be found in the original paper.The shaded areas are the error of the mean at each radii, and in fainter colours the azimuthal average of the spectral index error maps produced by tclean.The bottom right inset shows the geometric mean of the beam size for each spectral index map (Table2) and the vertical grey lines represent the systematic error introduced by the flux calibration uncertainty.Horizontal dotted lines indicate  = 2 (optically thick emission in the Rayleigh-Jeans limit) and  = 3.7 (typical value for micrometer dust grains in the ISM, seeFinkbeiner et al. 1999).The vertical dashed lines represent respectively the dark annulus visible in Fig.1and the water snowline position(Tobin et al. 2023).

Figure 3 .
Figure 3. Dust density distribution as a function of particle size in the many-seeds scenario.The colormap refers to the surface density of dust in a given size range at a given radius.The four panels indicate different key times of the simulation: before the outburst (top-left), after an outburst duration of 10 yr (top-right), after an outburst duration of 100 yr (bottom-left), and after an outburst duration of 10 000 yr (bottom-right).The vertical white dotted line (H 2 O) specifies the current snowline position.The blue and green solid lines respectively represent the fragmentation and radial drift growth barriers.The solid white line indicates the size at which particles have a Stokes number St = 1.

Figure 4 .
Figure 4. Dust density distribution as a function of particle size in the resilient scenario.The colormap refers to the surface density of dust in a given size range at a given radius.The four panels indicate different key times of the simulation: before the outburst (top-left), after an outburst duration of 10 yr (top-right), after an outburst duration of 100 yr (bottom-left), and after an outburst duration of 10 000 yr (bottom-right).The vertical white dotted line (H 2 O) specifies the current snowline position.The blue and green solid lines respectively represents the fragmentation and radial drift growth barriers.The solid white line indicates the size at which particles have a Stokes number St = 1.
Figure6.Spectral index radial profiles in the many-seeds scenario (upper panels) and resilient scenario (lower panels) after increasing outburst duration at 0.88 − 1.3 mm for  T = 10 −2 (left panels, red curves) and  T = 10 −4 (right panels, blue curves).We include the pre-outburst spectral indices as dashed lines, and the observed spectral index with a solid grey line (Fig.2).

Figure A1 .
Figure A1.Dust density distribution as a function of particle size in the many-seeds (top) and resilient scenarios (bottom) for  T = 10 −2 (left) and  T = 10 −4 (right).The colormap refers to the surface density of dust in a given size range at a given radius.In each case, the four panels indicate different key times of the simulation: before the outburst (top-left), after an outburst duration of 10 yr (top-right), after an outburst duration of 100 yr (bottom-left), and after an outburst duration of 10 000 yr (bottom-right).The vertical white dotted line (H 2 O) specifies the current snowline position.The blue and green solid lines respectively represent the fragmentation and radial drift growth barriers.The solid white line indicates the size at which particles have a Stokes number St = 1.

Table 1 .
Summary of the ALMA observations.

Table 2 .
Final images and spectral index maps.