Viscous heating as the dominant heat source inside the water snowline of V883 Ori

FU Orionis-type objects(FUors) are embedded protostars that undergo episodes of high accretion, potentially indicating a widespread but poorly understood phase in the formation of low-mass stars. Gaining a better understanding of the influence exerted by these outbursts on the evolution of the surrounding protoplanetary disc may hold significant implications for the process of planet formation and the evolution of disc chemistry. The heating due to outbursts of high accretion in FUors pushes the snowlines of key volatiles farther out in the disc, so they become easier to observe and study. Among the known FUors, V883 Ori is of particular interest. V883 Ori was the first FUor to show indirect evidence of a resolvable snowline beyond 40 au. By introducing a radial-dependent model of this source including viscous heating, we show that active heating is needed to reproduce the steep thermal profile of dust in the inner disc of V883 Ori. Our disc modeling combines the effect of stellar irradiation and the influence on the disc shape caused by the outburst of accretion. The accuracy of our model is tested by comparing synthetic ALMA images with continuum observations of V883 Ori, showing that the model successfully reproduces the 1.3 mm emission of V883 Ori at high spatial resolution. Our final predictions underline the importance of viscous heating as a predominant heat source for this type of object, changing the physical conditions (shape and temperature) of the disc, and influencing its evolution.


INTRODUCTION
FU Ori objects(FUors) are young stellar objects with episodes of outbursts that occur when the accretion rates of the embedded protostars are enhanced by orders of magnitude (Hartmann et al. 1996).During these events, they reach accretion rates up to 10 −5 , 10 −4 M⊙ yr −1 (Hartmann et al. 2016b;Audard et al. 2014;Fischer et al. 2023).In particular, accretion outbursts also explain the low luminosity problem in the stellar formation process of low-mass stars.Assuming steady accretion, the inferred accretion rates from the luminosity in low-mass protostars are too low for them to reach their final masses within the circumstellar disc lifetime (Dunham & Vorobyov 2012;Kóspál et al. 2017).The short episodes of high accretion rates solve this problem and explain the low luminosities of the circumstellar discs during the low accretion rate states.Nevertheless, understanding the stellar assembly is an ongoing effort (Fischer et al. 2023).
Due to the outbursts of high accretion rate in FUors, we expect circumstellar discs to suffer significant changes in their physical conditions during this stage (Armitage 2011;Hartmann et al. 2016a).MacFarlane & Stamatellos (2017) model circumstellar discs with radiative and episodic feedback, finding that episodic accretion changes the thermal and density profiles of the disc.Moreover, Zhu et al. (2009) analyze how young stellar objects accrete, and the changes induced on thermal profiles by different accretion rates and stellar masses.Due to the high accretion rates reached by these objects, it is important to examine the relevant heat sources and how these sources shape the circumstellar discs and change their emission and evolution.
Observations with the Atacama Larga Millimeter Array(ALMA) of V883 Ori, a FUor, have shown that V883 Ori has a very bright and optically thick core in the submm continuum surrounded by an optically thin halo with an abrupt change of the spectral index at ∼40 au in millimeter wavelengths (Cieza et al. 2016).Cieza et al. (2016) associated the spectral index shift as an indirect tracer of the water snowline due to the sublimation of icy grains.Usually, the location of the water snowline for a passive disc of a T Tauri star is at a couple of au (Lecar et al. 2006).However, V883 Ori's snowline beyond ∼40 au is an order of magnitude farther than the expected value, which cannot be explained by passive heating.i.e., only from stellar irradiation.Moreover, using molecular line emission, the snowline location in V883 Ori has been proposed to be even farther out at 80 au (Tobin et al. 2023) and ∼100 au (van 't Hoff et al. 2018;Leemker et al. 2021).
Understanding the thermal structure of V883 Ori is essential for the interpretation of an increasing number of molecular line data on the object.The ongoing outburst in the V883 Ori disc has sublimated enormous amounts of volatiles typically trapped on ice mantles, rendering the source as a unique laboratory for the study of disc chemistry (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).
There are multiple origins that can potentially explain the outward shift of the snowline in FUors.Harsono et al. (2015) used radiative transfer simulations to study the snowline location shifts when a flattened envelope is included in the disc.However, they found that an envelope would only shift the snowline by a couple of au.A shift of a few au is not large enough to explain the case of V883 Ori.Another explanation of the spectral index change can be linked to a more significant free-free emission in the inner disc due to higher temperatures instead of sublimation of icy grains (Liu et al. 2017).Nonetheless, even in the presence of freefree emission, its main contribution in a circumstellar disc occurs at longer wavelengths than the ALMA band 6 observations (λ = 1.3 mm), making this scenario improbable for V883 Ori.Another possible reason for the snowline shift could be found in the accretion rate, which heats up the disc and raises the dust thermal emission.The inclusion of viscous dissipation can produce shifts of tens of au (Min et al. 2011;Lecar et al. 2006), significantly boosting the brightness temperature and emission in the inner regions of actively accreting discs (Labdon et al. 2021;Ueda et al. 2023).Therefore, we can explain the steep thermal profile in V883 Ori by including viscous heating in the disc.
The accretion outburst does not only change the temperature in FUors, it also changes the disc's physical structure.Bell et al. (1997) show that viscous heating causes the flaring index to be variable in FUors, i.e., their aspect ratio has considerable radial variation, changing the incident angle of the stellar irradiation.For high accretion rate scenarios, the inner part of the disc follows a flared disc structure produced by the sharp radial gradient in viscous heating.However, once the viscous heating stops being significant in the outer disc, the shape starts to change, producing a self-shadowing effect (Bitsch et al. 2014(Bitsch et al. , 2013)).Hence, we need to include the study of a non-constant flaring index structure in the modeling of this type of object.A variable flaring or shape in the disc can be key in explaining the brightness temperature and the steep dust continuum radial profile of V883 Ori.It will also have important effects on the chemical evolution of the disc by moving the water snowline (Molyarova et al. 2018;Notsu et al. 2022).
Here, we present the first model that reproduces the steep thermal profile for V883 Ori in millimeter dust continuum emission through the inclusion of active heating sources and heat diffusion in the disc which can play an important role in heat transport in discs (Ziampras et al. 2023).Our model incorporates viscous heating as the main heat source in the inner disc, where the millimetric optical depth (λ = 1.3 mm) in V883 Ori has been estimated to be high enough to be emitting as a blackbody at radii up to ∼40 au.Since the thermal structure and the shape are dependent on each other, we developed a radial-dependent algorithm that fits the millimeter dust continuum emission profile at each location.The algorithm fits the continuum intensity by modeling the physical temperature, its scale height, the dust distribution of the disc, and its accretion rate.
This paper is structured as follows.In Section 2, we describe a two-layer analytical model for the disc using viscous heating and we present the model that we used as input for the radiative transfer code RADMC-3D (Dullemond et al. 2012).In Section 3, we describe the results of the radiative transfer model for the disc structure and the dust temperature together with the predictions for the millimeter emission.We discuss limitations in our model and further effects of viscous heating in Section 4. Finally, we summarize our findings in Section 5.

Two-layer Model Approach
To reproduce the dust emission in V883 Ori, our first modeling approach is done with a passive disc, i.e., only with irradiation heating.The model is resolved with an MCMC optimization of the parametric model from Cieza et al. (2018).The parametric model uses the self-similar solution of Lynden- Bell & Pringle (1974) for the dust surface density: where Σc is the dust surface density at the characteristic radius Rc.This first approach is done considering a constant flaring index ψ, i.e., the scale height of the disc is determined by the following relationship: We started with a standard flaring index ψ = 1.15 different from the fit from Cieza et al. (2018) because a flaring index lower than one causes the disc to constantly self-shadow and cool down.This low flaring index is mostly produced by the overcompensation of the fit due to the lack of active heating sources in the disc.A more detailed explanation for the parametric disc model and the results of the MCMC optimization can be found in Cieza et al. (2018).Before adding the viscous heating to the disc, we took the passive disc parameters from the MCMC optimization in Table 1 as the starting point of our model, finding that they were unable to recreate the centrally peaked and steep emission in the inner disc even with a variable aspect ratio.Thus, we add the viscous heating and let the accretion rate be variable with radius, so our model also predicts the vertical structure of the disc and the radial dependence of the accretion rate.We estimate the radial parameters and associated uncertainties using a bootstrapping method sampling values from a distribution from a self-iterative method (more details in Appendix A ).
V883 Ori is a young Class I source with a massive disc and a thin envelope, so the stellar size is quite uncertain.Nevertheless, we set its stellar size to R * = 2.5 R⊙ following the model from Gramajo et al. (2014).V883 Ori is part of the Northern region of the L1641 cloud cluster (Reipurth 2008), which is close to the Orion Nebula Cluster (ONC).Thus, we adopted a distance of 388 pc to V883 Ori (Kounkel et al. 2017;Connelley & Reipurth 2018).
The dynamical stellar mass from PV diagrams is consistent with M * = 1.3 M⊙ (Cieza et al. 2016;Ruíz-Rodríguez et al. 2022), while the starting disc mass M d = 0.6 M⊙ was obtained from the MCMC optimization made by Cieza et al. (2018).V883 Ori bolometric luminosity, ∼ 400 L⊙ (Strom & Strom 1993), is converted to a stellar accretion rate of Ṁ * ∼ 7 × 10 −5 M⊙ yr −1 .Liu et al. (2022) found a slightly larger accretion rate ∼ 1.1 × 10 −4 M⊙ yr −1 by modeling the observed energy spectrum.These observed accretion values are similar to the one inferred for FU Ori itself, which is ∼ Ṁ * ∼ 4×10 −5 M⊙ from a similar fitting procedure (Pérez et al. 2020).Despite the fact that it has been proposed that V883 Ori started its outburst of accretion over a century ago (Strom & Strom 1993), recent observations estimate its luminosity to be 212 L⊙ (Connelley & Reipurth 2018), which is probably showing the typical decay of FUors over time.Regardless of its luminosity decay, V883 Ori's luminosity is still very high for a 1.3 M⊙ protostar, luminous enough to consider viscous heating as an extra heat source.

Assumptions in our two-layer model
We follow a two-layer model to reproduce V883 Ori's emission.One layer is the surface or photospheric layer, exposed to stellar irradiation; and the second layer is the midplane.We list the mathematical symbols used in the formulation of our model in Table 2.In our modeling, we assume the following conditions: (i) The disc is geometrically thin enough to be considered locally plane-parallel.
(ii) The disc is axisymmetric.
(iii) The disc is in local thermodynamic equilibrium and the viscous heating is distributed following the density vertical structure.
(iv) We assume the radiation field to be isotropic, so we can use the Eddington approximation, i.e., the net radiation flux in the midplane is zero.
(v) For simplicity of the calculations we do not consider scattering in the analytical formulation.

Heat Sources present in Active discs
We consider three heat sources in our model.The first one is the irradiation heating that comes from the central star.
The second heat source is the viscous heating produced by the accretion rate, and the third one is the background component coming from the interstellar radiation field.The last term is negligible compared to the other two heating sources.Nonetheless, we add it to set a constraint in the lower limit of the disc temperature.
The first heat source, irradiation heating, takes into account stellar irradiation and its dependence on the disc shape.Using the thin disc assumption and considering that the disc is being irradiated on both sides, the irradiation received by each side is approximately (Fukue 2013): where r is the local radius from the star, R * the stellar radius, h the local scale height and hp the photospheric scale height.
We introduce the variable Tirr in Eq. ( 3 ) and its local thickness ( h r ), so it changes with the shape of the disc.In Equation (3), we add the factor hp, or the height of the visible photosphere.The photospheric scale height was used by Chiang & Goldreich (1997) to describe the location of the photosphere in circumstellar discs around protostars.In our model, the photospheric scale height is proportional to the scale height of the disc.At radii between 1 to 100 au, it is considered to be between 4 − 5 times the scale height of the disc (hp ≈ 4 − 5h).In particular, our model uses hp = 4h.
Besides the irradiation heating, we add the expression associated with viscous dissipation.The viscous heating for a standard α-disc (Shakura & Sunyaev 1973) for each side of the disc in a ring of width δr is (Frank et al. 2002): where Σ is the surface density, ν the viscosity, r the radius and ΩK = GM * /r 3 the Keplerian angular velocity.
Taking the flux energy density of the viscous heating in Eq. ( 4), we derive the expression: Different radial and vertical prescriptions can be found depending on the disc structure assumed for the parameters Σ and ν.These parameters are difficult to constrain observationally, so we use the accretion rate instead, which also diminishes the number of free parameters in our model.To express the viscous heating as a function of Ṁ , we use the equation relating the accretion rate with the viscosity and gas surface density for an accretion disc: The viscous heating then becomes a function of the accretion rate and the Keplerian angular velocity.Combining Equations ( 5) and ( 6) the final equation for the viscous heating is the following: In Equation ( 7), Tacc represents the viscous heating as an effective accretional temperature.In this case, most of the viscous dissipation will be embedded inside the disc, below the photosphere.By using the diffusion approximation, the viscous heating will eventually reach the upper layers, heating them up.Given that it comes from the potential energy released by the gas being accreted, we assume that the vertical distribution of the viscous heating follows the density distribution.We have to multiply our viscous heating expression in Equation ( 7) by two to account for both disc halves.Therefore, the viscous energy dissipation per unit volume per unit mass in our model is: where ρ(r, z) is the density and Σ(r) the surface density of the disc.The surface density in our model is a free parameter, but we expand it vertically assuming hydrostatic equilibrium.Hence, we assume a Gaussian distribution for the vertical structure of the disc using the scale height as the dispersion, i.e., ρ(r, z) As expected, the accretion rate is closely linked to viscous heating.Equations ( 7) and ( 8) show a linear dependence between the viscous heating and the accretion rate.If we consider a constant accretion rate, we would have a radial dependence of the form q + ∝ Ω 2 K ∝ r −3 .Because of the cubic dependence of the viscous heating, its energy input will decay faster with radius than the stellar heating that follows the square-power law. Figure 1 shows the comparison between stellar irradiation heating and viscous heating for a star with T * =7000 K and R * = 2.5R⊙.Under highly viscous conditions, accretion rates of FUors can reach values > 10 −4 M⊙ yr −1 (Zhu et al. 2009;Audard et al. 2014).When the accretion rate is of the order of 10 −4 M⊙• yr −1 , the viscous heating overrules any other heat source to distances of the order of tens of au or even a hundred au.
We also define a temperature set by the background radiation.This background radiation field sets a lower limit for the disc temperature, and is added as an external heat source as follows: where we set T bg = 10 K.
Given the large thermal gradients produced by viscous heating and the high density of the disc, we included radial and vertical diffusion in the disc midplane as key energy transport mechanisms in addition to the previously mentioned heating sources.We describe the implementation of the radial and vertical diffusion in further detail in Sections 2.1.4and 2.1.5respectively.

Effective Temperature
We assume that the dust emission is thermal, i.e.Sν = Bν .Thus, we consider the effective temperature as the brightness temperature at millimetric wavelengths.For the optically thick regions (τ ≳ 1) we use the brightness temperature as the dust temperature, while for the optically thin region, we correct the continuum attenuation by optical depth effects.
The input energy at a given radius is the viscous heating ) and the viscous heating for three different constant accretion rates in a 7000 K star with 1.3 M ⊙ .We observe that viscous heating dominates in the inner regions of discs with strong accretion (r <10 au).
from the disc's interior plus the irradiation heating from the star and what we have called the background radiation from the medium.The energy output at a certain location is just given by T eff : The left-hand side of Eq. ( 10) is the input of energy in the disc and the right-hand side is the continuum emission.From the energy conservation equation (Eq.( 10)) we obtain the relation for the radiative equilibrium of the disc.
For active discs, such as V883 Ori, viscous heating comes from the release of potential energy produced by the accretion rate.The fact that the dust emission is so steep and passive heating is not able to explain the brightness temperature means that T 4 acc ≫ T 4 bg + T 4 irr .We assume that the brightness temperature(TB) of the emission in the optically thick region should match the accretional temperature of the dust within a small margin, i.e., we enforce the condition that T eff ≈ Tacc as our initial estimation of the local accretional temperature at each radius.So, we model the dust emitting as a blackbody at a temperature T = Tacc, i.e., Then, we use the accretional temperature to solve the accretion rate using the viscous dissipation in Eq. 7. Therefore, the accretion rate that we infer from the brightness temperature in the inner disc is: With that parameter, we estimate a local surface accretion rate necessary to achieve the dust emission once the effects of irradiation and background heating are considered simultaneously for all radii, including the disc shape dependence of the irradiation term.
After we obtain the viscous dissipation from the brightness temperature, we calculate the midplane temperature using radiative transfer equations with the assumptions previously detailed to calculate the vertical and radial diffusion.
We then use the midplane temperature to study the disc shape and scale height by assuming hydrostatic equilibrium.

Vertical Heat Diffusion
As opposed to passive discs, the midplane temperature in accretion discs can reach temperatures >1000 K (D' Alessio et al. 2005).To calculate the effect of viscous heating in different layers of the disc, we also consider the optical depth as a driver of heat diffusion.By using the dust surface density together with the Rosseland mean opacity, κrs, we calculate the optical depth at each radius of the disc.We consider the temperature-dependent opacity prescription from Semenov et al. (2003) for the heat diffusion in the disc.Once the opacity and dust surface density are accounted for, the optical depth to the midplane is calculated using τ mid = κΣ 2 .For the optically thick inner regions, the disc is much hotter in the midplane (where the viscous heating has been input) than in the disc atmosphere, showing an expected thermal inversion (Calvet et al. 1992;D'Alessio et al. 1998).With such a vertical thermal gradient, the disc is heated up from the midplane to the surface instead of typical passively heated discs.
To solve the vertical structure of the disc temperature we employ the radiative transfer equations for a gray atmosphere in the vertical axis, as in Hubeny (1990).Considering µ = cos(θ), with θ the angle of the incident radiation, we use the standard intensity moments: Jν as the mean intensity, Hν as the Eddington flux, and Kν as the vertical projection of the radiation stress tensor or K-integral, i.e., and Using the plane-parallel approximation for a thin disc with no scattering term, the radiative transfer equations become: where ρ is taken as the gas density, κν the opacity, and jν the dust emissivity.We apply the Eddington approximation (radiation field is isotropic on each hemisphere) to substitute Kν by the mean intensity Jν with an Eddington factor fK = K/J = 1/3, i.e., Then, we employ the energy conservation equation between matter and radiation for a given (r, z) to solve the radiative transfer equation and get the Eddington flux: where Γ(r, z) is the vertical distribution of the viscous heating given by Eq. ( 8) and κJ = κν Jν dν/J.We resolve the vertical structure in the disc by integrating Equation ( 17) in frequency space combined with Equation ( 20).The integrated-frequency version of Equation ( 17) is the mean vertical energy flux, which is Γ(r, z)/4π.Integrating Equation ( 17), from the disc surface down to the midplane, leads us to the Eddington flux as a function of vertical depth in the disc: We use the definition for optical depth in our vertical reference system from the surface to the midplane, to solve Eq. ( 23) as a function of τ for simplicity.Before finding the vertical thermal structure of the disc we derive the expression for the Eddington flux at each optical depth: where H(0) is the Eddington flux at the surface, and τz, is the optical depth measured from the disc surface to the midplane.H(0) is given by the energy balance between the effective flux from the energy being diffused from the disc's interior from the viscous heating diffusion, Fin; and the external energy input at the disc surface from the stellar irradiation and the background field, Fout.Using Equation 10, the expression for H(0) is then: Then, we use the optical depth definition in Equation ( 24) and the Eddington approximation between K and J in Eq. 19 to integrate Equation ( 18) through all frequencies to solve for J: Then, expanding H(τ ) from Eq. ( 25) in Equation ( 30) we get the following relationship: Our mean intensity field J at the midplane, which is directly related to the vertical diffusion Fz, results in: To solve for the temperature at the midplane we apply the following boundary condition at the surface of the disc including external heating sources (Fukue 2013).The boundary condition assumes a Lambertian reflective surface, i.e. the reflection is isotropic, and it takes into consideration the irradiation effects on the disc.1 : The final expression for the vertical diffusion flux given our assumptions is: To calculate the optical depth in the millimeter regime of the inner disc, we use a mean opacity value without scattering κmm = 0.02 cm 2 • gr −1 (Beckwith et al. 1990).This opacity value is used to calculate the physical temperature form the observed brightness temperature.

Radial Diffusion
In the thick midplane, heat diffusion has a radial and a vertical component, any azimuthal component would be zero given our assumption of axisymmetry.We still have radial heat diffusion as an additional component of our model.Following the derivation in Casassus et al. (2019), the radial heat diffusion at a given radius r is: where ρ0 is the density in the midplane, κrs the Rosseland opacity, dT r the radial thermal gradient at a given radius r, T the midplane temperature as a function of r, and h the local scale height.Radial diffusion attempts to smooth strong thermal gradients and discontinuities in the disc that can be caused by other effects, such as abrupt changes in dust opacity or physical star-disc interactions.

Midplane temperature in the optically thick regions
Once we have all the heating and diffusion terms, we calculate the midplane temperature considering both, the vertical term from Equation 36 and the radial diffusion term from Equation 38 balancing the energy input and output in the disc midplane.Thus, the midplane temperature, T mid , is: Once T mid is calculated, the local scale height is updated accordingly, which will change the calculated values for the aspect ratio, h r , and the flaring index.We iterate the procedure in order to get self-consistent results until |T model − T brightness | <0.1.

Outer disc, the Optically Thin Regime
In the outer disc, the measured brightness temperature is lower than the temperature of an irradiated disc, confirming that the disc has to be optically thin in this region (τmm ≲ 1).Therefore, we model dust emission differently.Given that the disc is more diffuse at this region, we apply a simplifying approach assuming that the disc is vertically isothermal, i.e., the midplane and surface temperature are the same.
We measure the optical depth in the outer disc by solving the radiative transfer equation Iν = Bν (1 − e −τν ), from the previous fitting of V883 Ori's physical structure in the passively heated region.Therefore, the expression for the millimeter optical depth, τmm, is:

Radiative Transfer Simulations
We run radiative transfer simulations using the derived twolayer model including active heating to get the thermal structure of V883 Ori.We simulate the radiative transfer using RADMC-3D (Dullemond et al. 2012) adding an additional heat source in the disc, as it was implemented by Hord et al. (2017).This approach includes the amount of energy added by the external heat source at every cell of the grid.We include the extra heat source using the expression in Equation ( 7) and the inferred accretion rate from the two-layer model fitting.We then expand it vertically as a 3D flux density using Equation ( 8).Despite the fact that the vertical distribution is not well known, changing its vertical distribution between a uniform and a gaussian did not significantly change the synthetic millimeter continuum emission.The numerical simulations are run in a numerical grid of 500 × 256 × 1 cells in (r, θ, ϕ) respectively.The radial cells follow a logarithmic sample ranging from 0.8 au to 180 au.The colatitudes values are linearly sampled from 0.1 to π 2 assuming mirror symmetry to the other side of the midplane.The input density field of the radiative transfer simulations is the one obtained from the two-layer model.The gas and dust are vertically distributed with a gaussian centered at z = 0 and the scale height as the dispersion.We use a variable flaring including the shadowing effect obtained from the self-consistent fit.
We calculate the thermal equilibrium using 10 8 photon walkers for each run.Dust opacity is modeled using two different dust compositions.The first dust composition corresponds to dry grains with mass fractions of 50% organics, 41% astrosilicates, and ∼ 8.8% troilite.The second composition, the icy grains, has a water mass fraction of 20%, 40% in organics, 33% in astrosilicates, and ∼ 7% in troilite.The opacity values were taken from Henning & Stognienko (1996) for the troilite and organics composition, from Draine (2003) for the astrosilicates, and Warren & Brandt (2008) for water ice.We set the radiative transfer with three iterations to update the dust populations for self-convergence.Wherever the dust temperature is lower than 150 K, the dust population is composed of icy grains, and if the temperature is higher than 150 K, it is composed of dry grains.
The dust size distribution follows a standard power law proportional to the dust grain size: with a the size of the dust particle (Mathis et al. 1977).The dust grain sizes range from amin = 0.05 µ m to amax = 2.5 mm for the icy large grains, while amax = 0.5 mm for the dry grains inside the snowline, simulating an enhanced dust growth due to the water ice stickiness in the outer disc.
Once the thermal equilibrium is calculated, we generate a synthetic continuum image for λ = 1.3 mm.After that, we compare it with ALMA band 6 observations of V883 Ori by convolving the resulting image with the same elliptical beam of the observations.The elliptical beam of the observations has a major axis of ≈ 37 mas, a minor axis of ≈ 27 mas, and a position angle of 52 degrees.

Radial Dependence of Flaring and Aspect Ratio
The solution in Figure 2 shows the self-shadowed regions of the disc and its variable flaring.There is a decline in the aspect ratio beyond the snowline where the accretion rate also starts to decrease.The change in the disc morphology is mainly caused by the decrease in the effect of viscous heating and the change in the opacity of the dust grains.As the disc gets hotter, it also gets geometrically thicker so the flaring index actually decreases when the viscous heating starts to decay.When the flaring index is less than one, the disc does not receive direct stellar irradiation.Thus, behind the hot and thick region heated by viscous heating, the disc gets colder due to self-shadowing, having additional strong implications for the photochemical reactions in these regions.in The flaring index in Figure 2 shows abrupt changes.Those jumps can be explained by the significant opacity changes which are temperature-dependent (Semenov et al. 2003).One of those changes is strongly related to the water snowline.As water sublimates from icy grains, the opacity suffers drastic changes which translates into an abrupt jump in the disc shape as well, which is illustrated by the flaring index.
In the self-shadowed region, Figure 2 shows a thick disc between 10 and 30 au.This thick region causes the outward zones to be shadowed by stellar radiation.Additionally, our After 50 au, the viscous heating is negligible so the disc behaves as a passive disc and the accretion becomes uncertain at that point.We also include a comparison between the measured in-fall velocity and the gas free-fall velocity to verify that gas and dust are coupled so dust temperature is representative of the viscous dissipation.
model predicts a wiggle in the accretion rate.The accretion rate reaches its peak in the plateau and starts an abrupt decay afterward, right at 40 au.The change in the disc's accretion rate could mean that there is an accumulation of material at the location of the spectral index shift.Such an effect could be key to understanding the variability of FUors and the effect of snowlines on dust evolution and planet formation.Moreover, shadowing regions produce variability through thermal-wave instability and radial evolution in the location of the snowlines (Okuzumi et al. 2022).

Effects of Viscous Heating on Thermal Structure
The first consequence of the addition of viscous heating is the increase of the dust temperature to values ∼ 1000 K close to the star (see Figure 2 ).The high gas density and the uniform vertical distribution of viscous heating in the disc produce an almost constant thermal structure at different layers of the disc.
The results of the two-layer fit in Figure 2 show that the viscous heating is substantially more important at the inner 40 au.When the viscous heating is included in V883 Ori, the  midplane becomes hotter than the surface layers.We also observe that the viscous heating increases the temperature at the disc midplane as far as 40 au from the star.For the outer part of the disc (r > 70 au), the viscous heating stops being significant and the disc is almost isothermal in the vertical direction, so our initial assumption for the outer disc seems reasonable.
Figure 3 shows the results of the fitting employing the bootstrapping method (more details in Appendix A2) with the respective uncertainties.Overall, the bootstrapping method matches the fit from the self-consistent approach, with a very good constraint of the dust temperature in the midplane.We also notice that even though the surface density matches the assumed input, the uncertainties are larger, and therefore the dust distribution is less constrained in the optically thick region of the disc.
We illustrate the input dust density meridional structure in Figure 4 showing the radially variable structure of the dust distribution in the V883 Ori.From comparing the effects of viscous heating in the disc temperature in Figure 5, we observe that the 115 K layer is shifted by 20-30 au, almost matching the spectral index change in the disc.The role of viscous heating in pushing the condensation fronts outwards is proved by the radiative transfer models.Moreover, it is observed that the viscous heating affects the disc mostly in the inner regions(<40 au), while the outer regions remain without significant thermal changes.At those radii, the irradiation heating dominates and the typical vertical structure of passive discs is recovered.
The simulations without viscous heating produce a dust temperature lower than the brightness temperature at almost every layer.Since the brightness temperature is a lower bound for the dust temperature, the conclusion of an extra heat source for the best fit of the data is reinforced.Irradiation heating by itself is not enough to explain the spiky radial continuum profile or the high brightness temperatures in the inner regions of V883 Ori (r≲ 40 au).
Due to the high midplane temperature caused by viscous heating and the high surface density of V883 Ori, we compare the different heating terms in order to gauge the relevance of radial diffusion in our model.Figure 6 shows the comparison between the radial heat diffusion (see Eq. ( 38)) and the other heating terms, i.e., the viscous heating, q + , and the irradiation heating, σSBT 4 irr .The Figure proves the fact that radial diffusion is comparable to vertical diffusion when the radial thermal profile of the disc is steep, attempting to counteract radial temperature differences and smoothing the radial thermal profile in both directions.The radial heat diffusion smooths strong thermal gradients that can be caused by obscured regions by self-shadowing effects.Nevertheless, as the disc thermal gradient flattens due to viscous heating decay and the disc starts to become optically thin, the radial diffusion flux becomes negligible.

Accretion Rate as a Radial Function
We show the local accretion rate in the disc as a function of radius in Figure 2. We observe that the accretion rate reaches a peak value ∼ 10 −3 M⊙• yr −1 at the location of the plateau (30 < r < 40 au).Beyond the plateau, the accretion rate and the viscous heating are less strong, and the disc is passively heated.The fact that the accretion rate is not continuous throughout the disc means that there must be a pile-up of material.However, given that V883 Ori is a FUor, we expect that it has gone through a sequence of accretion outbursts, explaining the differences in the local accretion rate.
As a safe check, we also measure the radial velocity from the accretion rate, which is shown in Figure 2. Our measured radial velocity is below the free fall velocity, even at the location where the peak of accretion occurs.Therefore, the gas motion is not related to virialization and the assumption of dust-gas coupling seems reasonable.So, the assumption of using the continuum emission as a tracer of viscous heating is reliable.We calculate the accretion luminosity produced by viscous heating in the inner disc as another one of our safety checks.This accretion luminosity should not be higher than the measured bolometric luminosity for V883 Ori ,L * ∼ 400 L⊙ (Sandell & Weintraub 2001).To calculate the accretion luminosity in the disc, we integrate the viscous heating through concentric rings at each radius, i.e., The integrated accretional luminosity from Equation 44is Lacc ≈ 80 L⊙.Even if we consider that roughly half of the accretion luminosity should be emitted very close to the star inside the dust cavity, our integrated luminosity is less than 20% of V883 Ori's bolometric luminosity.Thus, it is possible then that we are underestimating the contribution of other emitting components in V883 Ori that are not probed by millimeter wavelengths.Another probability is that the accretion luminosity in the inner disc is highly variable, which might explain why different measurements taken at different epochs give different values (e.g.L bol = 200 vs 400 L⊙).

Synthetic Predictions for V883 Ori
The synthetic predictions of the continuum images ( λ = 1.3 mm ) in Figures 7 and 8 match the observations.The Figures Midplane Temperature (K) Figure 6.Ratio between the radial flux inside and the heating terms considered in this work.For the self-shadowed region, the radial heat diffusion will be comparable to the heating terms so it is expected that it will raise the disc temperature.Radial heat diffusion seems to be important right outside the snowline as well.However, our opacity values change drastically in and out of the snowline due to the ice sublimation, causing the radial heat diffusion to be more significant.We observe that our model with viscous heating is the best match for the observations while the passive heating is not enough to reach the observed emission.
illustrate that viscous heating causes a considerable change in dust emission when it is compared to the passive heating model.The viscous heating produces a centered increment in the intensity that decays at outer radii.The effects of viscous heating translate into a steeper emission radial profile, such as the one observed in V883 Ori.
The MCMC fit done by Cieza et al. (2018) considers a model with a constant flaring index and without the addition of viscous heating.In general, the omission of the viscous heating for sources with typical accretion rates Ṁ <10 −6 M⊙ yr −1 will not produce considerable changes in the thermal structure beyond ≈5 au.However, for rapidly accreting objects like FUors, in particular for V883 Ori, viscous heating will strongly affect the disc shape, their thermal structure, and therefore the dust continuum emission.
Overall, we reproduce the inner disc emission.However, we do not explain the wiggle at the border of the snowline.The wiggle could be produced by changes in the dust density or different dust populations due to dust evolution.Wiggles and phenomena like these have been predicted before in the literature (Banzatti et al. 2015;Okuzumi et al. 2016;Zhang et al. 2015Zhang et al. , 2016)).Those studies report the presence of density wiggles and dips before and after the snowlines, finding a greater concentration of dust just outside of them, and a small depletion on the inside.Such radial density fluc-tuations could then explain the observed wiggle in the dust emission.However, understanding the specific changes in the dust composition around the snowline is beyond the scope of this paper.
Our model also has a sharp increment of the optical depth at r∼ 40 au.Besides a possible change in the opacity regime caused by the melting of the ice grains, another possible treatment for this issue could be a differentiation in the dust species.Notsu et al. (2022) show that there is a chemical enrichment of the ices behind the water snowline due to the recondensation of many volatiles, changing their composition.The change in dust population after the snowline produces a local change in the dust opacity.A thinner disc inside the snowline locally decreases its emission and heats up the dust right beyond the snowline by increasing its exposure to stellar irradiation.Thus, the intensity wiggles at 1.3 mm can also be explained by a dust population composed of larger ice grains beyond the snowline.

DISCUSSION
From the analysis of molecular line emission, van 't Hoff et al. (2018); Tobin et al. (2023) have placed the location of the snowline between 80 and 100 au, while Cieza et al. (2016) estimated it ∼ 40 au from the dust continuum emission.The discrepancies in the location of the water snowline might be related in part to distinct scale heights traced by different observations together with the complex thermal structure of an outbursting source such as V883 Ori as shown by the results of our modeling.
We have shown that V883 Ori, as a FUor, has significant changes in its intrinsic structure compared to passively heated discs.The disc shape, accretion rates, and the shift of the snowline make FUors particularly interesting for the study of their chemical and physical processes (Molyarova et al. 2018).Changes in the shape and the additional heating sources have significant effects on the disc emission even beyond the main region of influence.Figure 8 shows that viscous heating also increases the millimeter emission beyond 100 au.When viscous heating is not included, matching the brightness temperature causes the disc to be severely flared in the inner disc to match the emission.This compensation produces an important shadowing effect toward the outer disc making it cooler.In addition, a more extended hotter inner disc diffuses heat toward the outer disc, heating it further.
The V883 Ori disc structure and physical properties make it a suitable source for the study of the composition of the material that takes part in the planet formation process and the effects caused by disc evolution.For example, van 't Hoff et al. (2018) and Lee et al. (2019) were able to characterize the detection of organic molecules far out in the disc, and Ruíz-Rodríguez et al. ( 2022) did a comprehensive study of the emission of multiple chemical tracers in V883 Ori, while Tobin et al. (2023) resolved the deuterium hydrogen oxide(HDO) emission coming from the disc.Multiwavelength studies have the potential to shed light on the effect that the snowline has on dust populations and their possible role in the formation of pebbles and planet formation.
Our results show a variable accretion rate with radius.
The recovered accretion ranges between different values, from by using SED fitting (Liu et al. 2022), which are slightly lower than Ṁ * ∼ 7 × 10 −5 M⊙ yr −1 using the bolometric luminosity to ∼ 1.1 × 10 −4 M⊙ yr −1 the peak of our retrieved accretion rate profile.A variable accretion rate with a radius least to pile-ups of material in the local minima of the radial profile.Variable accretion rates, hence, pile-ups of material are still consistent with the characteristic accretion outburst behavior of FUors.Pile-ups can trigger the MRI and the Thermal Instability.Moreover, theoretical models predict that the outbursts are highly variable and episodic in nature (Zhu et al. 2009;Zhu et al. 2010;Bae et al. 2013;Hartmann & Bae 2017).
Even though viscous heating explains the steep radial trend in emission at the inner regions of the disc, the continuum emission has a wiggle at the boundary of the snowline.Changes in the composition and size distribution of the dust populations could be a possible solution to such variation in the radial emission profile.The sublimation of the icy mantles on dust grains makes dust growth less efficient, favoring the formation of smaller grains inside the snowline.Therefore, we expect changes in the dust opacity due to different dust populations or modulations of the surface density in and out of the snowline, producing the wiggle in the dust continuum emission.Water ice can increase dust stickiness, enhancing dust growth and changing the dust size distribution in the disc.The thermal inversion can produce distinct footprints in molecular absorption in infrared wavelengths in the inner disc (Calvet et al. 1992) or in the brightness radial profiles of different ALMA continuum bands for example.
One of the caveats in our modeling approach is that we assumed the vertical density distribution of a passive disc.However, the vertical distribution may not follow a Gaussian distribution.Uzdensky (2013) described the vertical distribution of an active disc powered by the Magneto-Rotational Instability(MRI).In their assumptions, assuming hydrostatic equilibrium, the distribution of the active disc is proportional to ∝ (1 − z 2 z 2 c ) 3 .Even though the assumed prescription for the gas distribution will produce changes in the vertical temperature as well, the midplane temperature will only have changes within the same order of magnitude, while preserving a diffusive approach.Changing the heating input uniformly or following the gas density distribution did not make significant changes to the millimeter brightness temperature of the disc or the overall reported trends in this work, although it will impact other disc parameters locally.
Another source of uncertainty in our work is the distance to V883 Ori.The Gaia catalogs flag the astrometric solution of V883 as unreliable.Thus, we adopted the distance of 388 pc from Connelley & Reipurth (2018); Kounkel et al. (2017), also consistent with the value derived by Lee et al. (2019) from neighboring sources.Despite the points mentioned above, using the Gaia distance to V883 Ori (d∼ 270 pc) does not qualitatively change our results.Such change would have mainly affected the values of the disc and stellar parameters, such as their masses, sizes, and the location of the snowline.If we have used Gaia's distance as the distance to V883 Ori, the snowline would have been at ∼ 26 au, still farther than expected values.Thus, the need for an extra heat source besides passive heating still remains.

SUMMARY
In this work, we report that when viscous heating is included in V883 Ori, the midplane temperature in the optically thick region is higher than the surface temperature.In the first few au, viscous heating increases the temperature to the order of 1000 K and higher values.In fact, viscous heating is the principal heating source at the inner disc, more important than the stellar or flaring terms.However, for stars with lower accretion rates ( Ṁ < 10 −6 M⊙ yr −1 ), the stellar heating will still be the main heating source in most parts of the disc (see Fig. 1).
The results of our semi-analytical approach show a strong initial flaring, and then a plateau of the aspect ratio of the disc in the inner 30 au.After the initial 30 au, the disc starts to self-shadow.In that region, the disc has a flaring index lower than one, i.e., the aspect ratio diminishes with radius, as it is shown in Fig. 2. We conclude that viscous heating plays a significant role in shaping the disc structure.
We compare the radiative transfer predictions of the model with the observations in Figure 8. From such comparisons we can conclude that viscous heating provides the extra energy budget to reach the emission levels, showing that the addition of viscous heat increments the continuum emission as a central source close to the star by changing the radial slope of the thermal profile, making it steeper.
In summary, classical passive heat sources are not capable of reproducing the emission profile of V883 Ori, an additional heating source is needed.This extra source has the particularity of being significant in the inner part, centrally peak and it decays quickly after a few au.High accretion rates are typical of FUors and their associated viscous heating matches the necessary conditions to explain the dust continuum emission profile of V883 Ori.Viscous heating not only heats up the disc but also changes its physical structure.V883 Ori in particular also presents a wiggle right at the snowline boundary where dust evolution could play a major role in creating dust substructures.
); Tirr is the equivalent temperature for the incoming flux from the irradiation heating at one side of the disc.The first term in the right-hand side of the equation, is associated with the stellar size.The second term, ), is related to the flaring of the disc ( dh dr

Figure 1 .
Figure 1.Comparison between the stellar irradiation flux (σ SB T 4irr ) and the viscous heating for three different constant accretion rates in a 7000 K star with 1.3 M ⊙ .We observe that viscous heating dominates in the inner regions of discs with strong accretion (r <10 au).

Figure 2 .
Figure2.Physical conditions of our model.The midplane temperature looks hotter than the surface temperature, as expected given the optical depth effects.The flaring index and the aspect ratio show the disc self-shadowing and thickening generated by viscous heating.The flaring index shows a wiggle behavior caused by abrupt changes in the temperature-dependent opacity.After 50 au, the viscous heating is negligible so the disc behaves as a passive disc and the accretion becomes uncertain at that point.We also include a comparison between the measured in-fall velocity and the gas free-fall velocity to verify that gas and dust are coupled so dust temperature is representative of the viscous dissipation.

Figure 3 .
Figure3.Midplane temperature and dust surface density as a function of radius in V883 Ori.We show a comparison between the distributions obtained from the bootstrapping algorithm and the self-iterative method, which are in agreement with each other.

Figure 4 .
Figure 4. Dust density distribution in our simulations including viscous heating.It illustrates the variable disc shape and the selfshielded regions caused by including viscous heating.

Figure 5 .
Figure 5. Dust temperature in our models turning on and off the viscous heating.The blue line is an isothermal contour at T = 115 K in the model with viscous heating added and the cyan the same isothermal contour in the model without viscous heating.We use 115 K as the temperature at which Cieza et al. (2016) locate the spectral index shift in V883 Ori.Left: Image with viscous heating.Right: Image without viscous heating.We observe in the images that when we add viscous heating the snowlines get shifted by almost 20 au.

Figure 7 .Figure 8 .
Figure 7.Comparison of dust continuum images at λ = 1.3 mm between the observations and our models with and without viscous heating.Left: Observations.Center: Image with viscous heating.Right: Image without viscous heating.The images show that the addition of viscous heating matches the central peak emission of the observations, whereas the image without viscous heating is not able to match the central emission in the inner au of the disc.

Table 1 .
Parameters of Parametric disc in V883 Ori.

Table 2 .
List of symbols used in this work