Impact of Size-dependent Grain Temperature on Gas-Grain Chemistry in Protoplanetary Disks: the case of low mass star disks

Grain surface chemistry is key to the composition of protoplanetary disks around young stars. The temperature of grains depends on their size. We evaluate the impact of this temperature dependence on the disk chemistry. We model a moderately massive disk with 16 different grain sizes. We use POLARIS to calculate the dust grain temperatures and the local UV flux. We model the chemistry using the 3-phase astrochemical code NAUTILUS. Photoprocesses are handled using frequency-dependent cross-sections, and a new method to account for self and mutual shielding. The multi-grain model outputs are compared to those of single-grain size models (0.1 $\mu$m), with two different assumptions for their equivalent temperature. We find that the Langmuir-Hinshelwood (LH) mechanism at equilibrium temperature is not efficient to form H$_2$ at 3-4 scale heights ($H$), and adopt a parametric fit to a stochastic method to model H$_2$ formation instead. We find the molecular layer composition (1-3 $H$) to depend on the amount of remaining H atoms. Differences in molecular surface densities between single and multi-grain models are mostly due to what occurs above 1.5 $H$. At 100 au, models with colder grains produce H$_2$O and CH$_4$ ices in the midplane, and warmer ones produce more CO$_2$ ices, both allowing efficient depletion of C and O as soon as CO sticks on grain surfaces. Complex organic molecules (COMs) production is enhanced by the presence of warmer grains in the multi-grain models. Using a single grain model mimicking grain growth and dust settling fails to reproduce the complexity of gas-grain chemistry. Chemical models with a single grain size are sensitive to the adopted grain temperature, and cannot account for all expected effects. A spatial spread of the snowlines is expected to result from the ranges in grain temperature. The amplitude of the effects will depend on the dust disk mass.


Introduction
As planets form in protoplanetary disks, studying the physical and chemical evolution of their gas and dust content along the whole disk lifetime is important in order to investigate how planetary formation proceeds. In the early stages of disks (the socalled protostellar phase) small dust grains are essentially dynamically well coupled to the gas. When grain growth occurs (up to at least cm-sized particles) larger particles, which dynamically decouple from the gas-phase, settle toward the disk midplane. This vertically changes the gas-to-dust ratio in the disk. Moreover, grain growth also allows UV to penetrate deeper in the disk, changing its chemistry.
As a consequence of the complex interplay of physical and chemical processes, it is now demonstrated that the disk structure can be separated vertically in three distinct regions: from top to bottom, the irradiated tenuous disk atmosphere (above 3 − 4 H or scale heights) where the gas can be hot and is ionized or in atomic form, the warm denser molecular layer (1 − 3 H) where most molecules reside and below 1 H, the midplane. Beyond the Send offprint requests to: Sacha Gavino e-mail: sacha.gavino@nbi.ku.dk CO snowline, which is typically located at a radius about 20 au in a TTauri disk, the midplane is very dense and cold (≤ 20 K). This area is essentially shielded from the stellar radiation, with both low ionization and turbulence levels. In this cold region, most molecules are stuck onto dust grains which exhibit icy mantles that can still be processed by cosmic rays.
Thanks to ALMA, several high angular resolution and sensitive observations of disks orbiting either low or intermediate mass young stars (TTauri and Herbig Ae stars) have quantitatively confirmed this general scheme. CO observations of HD 163296 have clearly shown the CO gas-phase depletion in the dark, cold disk midplane (de Gregorio-Monsalvo et al. 2013), revealing also that the CO gas is located in the intermediate warm molecular layer. More recently, the ALMA observations of the edge-on TTauri disk called the Flying Saucer  directly show the molecular layer at intermediate scales (1 − 2H) with CS, the denser gas tracer being located at about 1H, below the CO emission (1−3H) which better traces the whole envelope of the molecular layer.
In the last 10 years, several astrochemical models have emerged that attempt to incorporate the observed complexity of disk physics in order to provide more accurate molecular abun-Article number, page 1 of 39 arXiv:2106.05888v1 [astro-ph.GA] 10 Jun 2021 dance and surface density predictions. While the very first models included only the gas-phase chemistry, gas-grain coupling was added to take into account firstly adsorption and desorption for a few molecules (e.g. Woitke et al. 2009) and then reactions on the grain surface (e.g. Semenov & Wiebe 2011;Walsh et al. 2014). Another step was finally achieved with models which include 3 material phases: the gas-phase and 2 phases for the gasgrain chemistry with grains modeled by an icy surface layer and an active mantle (Ruaud & Gorti 2019;Wakelam et al. 2019). Thermochemical models which calculate the density and thermal structures in a self-consistent manner have been also developed (e.g. Hollenbach & Gorti 2009). Recently, Ruaud & Gorti (2019) have also coupled such a thermochemical model (Gorti et al. 2011) to a 3-phase chemical model.
An important improvement in chemical models was to introduce more realistic dust disk structures by taking into account grain growth and dust settling. Grain growth reduces the dust cross-section and therefore increases the UV penetration while dust settling locally changes the gas-to-dust ratio. Grain growth was first introduced in disk chemistry by Aikawa & Nomura (2006). Dust settling is often mimicked by adding larger grains in the midplane (e.g. Wakelam et al. 2019).
Several authors have already studied the effect of multiple grain sizes onto chemistry. As an example, Acharyya et al. (2011) used five different grain sizes in their models but the grain temperature was the same for all grains. Following this work, Pauly & Garrod (2016) tested the impact of a simple temperature distribution, with T (a) ∝ a −(1/6) , a being the grain size, while Ge et al. (2016) explored a small range of temperatures (14.9-17.9) for the Cold Neutral Medium. More recently, Iqbal & Wakelam (2018) also developed a model handling different grain sizes.
On another hand, Akimkin et al. (2013) developed a disk model where they couple the dust grain time evolution (grain growth and settling) to chemistry. Another improvement was achieved with the development of grids of thermochemical models allowing multiwavelength fitting of dust and molecular lines such as DIANA (Woitke et al. 2019). However, these elaborate approaches use a single grain temperature which is either fixed or derived from a dust radiative transfer simulation while the temperature of a grain depends on its size (e.g. Wolf 2003), larger grains being colder. Chapillon et al. (2008) found in the disk of CQ Tau, whose temperature is well above the CO snowline temperature, an important CO depletion. They suggested that CO may remain trapped onto larger grains that are cold enough to prevent thermal CO desorption. Heese et al. (2017) have investigated the dust temperature variations with grain size and position (radius and altitude above the midplane) in a typical TTauri disk using the 3D dust radiative transfer code Mol3D (Ober et al. 2015). The dust temperature variations in the molecular layer appear to be significant enough to affect the disk chemistry. This dependency questions the applicability of elaborate, selfconsistent, thermochemical models that rely on a single grain size and temperature.
We explore here the impact of the variations of the grain temperature with their size onto the chemistry of a representative TTauri disk. Our goal is not to build a complete chemical model that incorporates the best of all previous studies, neither to make grid models that could provide relevant comparisons with the observations. Instead, we only intend to check the limits of an assumption that is used in most previous models, and explore the main consequences of dust grain temperature dependence with size.
For this purpose, we have coupled the 3D Monte-Carlo continuum radiative transfer code POLARIS (Reissl et al. 2016) to the NAUTILUS Multi Grain Code (Hersant et al. 2009;Ruaud et al. 2016) using the 3-phase version which also takes into account a grain size distribution (Iqbal & Wakelam 2018). We present the two codes and how we parametrize the gas and dust disk in Section 2. Section 3 deals with the disk model description and global results. We discuss in more detail the results in Section 4 and we state our conclusions in Section 5.

Model description
To estimate the impact of the grain temperature onto the chemistry, we couple the radiative 3D Monte-Carlo code POLARIS with the astrochemistry code NAUTILUS. This requires to independently assume a gas and dust disk structure, a gas temperature, a dust distribution and settling. Moreover, the coupling between the two codes requests for NAUTILUS a new method to calculate the dust extinction, the self and mutual shielding using the UV flux calculated by POLARIS in each disk cell. Figure 1 is a scheme summarizing how we proceed.

Disk Model
Our disk structure is derived from a simple disk model where truncated power laws describe most of the quantities: surface density, midplane temperature, velocity field, with sharp inner and outer radii. A detailed model description is given in Appendix A We adopt a power-law grain size distribution dn(a) = Ca −d da with dn(a), the number of grain of size a and C a normalization constant that can be derived from the dust mass (see Appendix A.2). The resulting distribution is discretized into 16 logarithmically distributed intervals, with grain sizes ranging from a min = 5 nm to a max = 1 mm. The fraction of mass and area per bin are given in Table A.1.
Our goal being to evaluate the impact of dust grain temperatures, we assume a common gas temperature structure for all our simulations (including single and multi-grain models), in order to avoid chemical effects that could be due to gas only. We selected values derived from the observations of the Flying Saucer TTauri disk (see Dutrey et al. 2017), whose edge-on orientation allow a direct measure of the 2D gas temperature distribution from ALMA observations of the optically thick CO J=2-1 line. Values are given in Table 1.
The vertical distribution of gas is self-consistently derived from the imposed temperature law. However, to account for dust settling, we assume grains of a given size follow a simple Gaussian vertical profile, whose scale height is simply related to the gas hydrostatic scale at the disk midplane temperature according to a simple prescription (see Appendix A.2). Our model is therefore not fully consistent, since very small grains should follow the vertical profile of the gas, but should be adequate for the larger grains observed in disks. To compute the settling, we assume an α viscosity of 0.01 and a Schmidt number of S c = 1. Although often assumed in disk modeling, these values are somewhat arbitrary, and observations with ALMA are not yet accurate enough to allow for a quantitative description of dust settling. Figure 2 is a 2D representation of the disk structure in density and in temperature for the gas and for the dust (multi-grain model with settling).  Physical structure of the multi-grain models incorporating a full dust distribution. T g is the gas temperature, n H is the gas number density, G is the local stellar and interstellar field in the case of the HUV models (see Appendix A), ζ is the dust to gas mass ratio, T a (Eq.1) is the area-weighted dust temperature, and n a is the area-weighted grain number density. A V is the visual extinction counted from disk surface toward the disk midplane.

Radiation fields and dust temperature
In this new approach, we use the dust radiative transfer code to calculate throughout the disk the UV field and the dust temperature for each of the 16 grain sizes, in each cell of the disk (in radius and altitude). The dust temperatures are then used to compute the chemistry with NAUTILUS.
However, the UV field derived above only accounts for dust extinction and scattering. A proper account for self and mutual shielding of molecules and atoms is required to handle the photochemical processes.

Radiation sources
As the central heating source we assume a low mass pre-main sequence star with a mass of M = 0.58 M , which radiates as a black body with an effective temperature of T ,eff = 3900 K and a luminosity L = 0.75 L . Additionally, we consider the interstellar radiation field (ISRF) with a spectral energy distribution (SED) from Draine (1978) between 91.2 nm and 200 nm and the extension of van Dishoeck & Black (1982) for longer wavelengths. Besides the stellar radiation field, TTauri stars exhibit a significant UV excess coming from the inner disk boundary and accretion shocks on the stellar surface. This contribution is highly wavelength-dependent, with spectral lines like Lyman α having substantial intensities (France et al. 2014). As concerns the thermal balance of dust, the frequency dependence is largely irrelevant, and can be absorbed by ensuring that the luminosity used in the model contains the UV excess contribution. However, photoprocesses being wavelength dependent, the details of the UV spectrum shape strongly affect the chemistry. For this reason, we take as stellar input spectrum for our radiative transfer Article number, page 3 of 39 A&A proofs: manuscript no. article  (innermost radius) 1 au R out (outermost radius) 250 au R 0 (reference radius) 100 au ζ (dust to gas mass ratio) 0.01 T mid,0 (observed T mid at R 0 ) 10 K T atm,0 (observed T atm at R 0 ) 50 K σ (stiffness of the vertical T profile) 2 q (radial T profile exponent) 0.4 S c (Schmidt number) 1 α (viscosity coefficient) 0.01 d (grain size distribution exponent) 3.5 a min (smallest grain radius) 5 nm a max (biggest grain radius) 1 mm number of bins 16 ρ grain (material density of dust) 2.5 g.cm −3 M disk (total disk mass) 7.5.10 −3 M and chemical model the sum of the stellar black body spectrum and the UV excess typically found in the in TTauri stars. Further details about the radiation fields are given in Section 3.1.

Radiative transfer simulations
The radiation field is computed using the 3D Monte-Carlo continuum and line radiative transfer code POLARIS (Reissl et al. 2016). As in Heese et al. (2017), we calculate the dust temperature distribution of 16 grain size intervals (see Fig. 5). We assume spherical grains with a size independent composition consisting of a homogeneous mixture of 62.5 % silicate and 37.5 % graphite. Optical properties from Draine & Lee (1984) and Laor & Draine (1993) are used to calculate the scattering and absorption coefficients based on Mie scattering (Mie 1908) using the MIEX code (Wolf & Voshchinnikov 2004). The radiative transfer is solved using 100 wavelengths logarithmically spaced between 50 nm and 2 mm. The spatial grid involves 300 radii, logarithmically spaced by a factor 1.05 between 1 and 300 au, and 181 angles with the same factor between 0 and π. The temperature in the upper optically thin disk layers strongly depends on the grain size (Fig. 5, middle panel). The dust temperature increases with grain size for sizes in the range 0.007 µm to 0.070 µm, then decreases for sizes in range 0.32 µm to 15 µm, and increases again, although only slightly, with larger grain sizes. As explained in Heese et al. (2017), this can be understood by looking at the ratios of the absorption crosssection for wavelengths at which the star radiates (≈ 1 µm) and the absorption cross-section for wavelengths at which the dust emits (≈ 20 µm, Fig. 4). The ratio increases up to grain sizes of 0.07 µm, meaning that the ratio of absorbed to emitted radiation increases, leading to higher grain temperatures. For larger grain sizes, the behavior is reversed until grains become larger than about 15 µm. Using these radiative transfer simulations we can obtain a distribution of dust temperatures T d (a i , r, z) that fully depend on the size and location of the grains. To compare our set  of models with those consisting of a single grain size (see Section 3.1.1), it is also convenient to introduce an area-weighted temperature defined as: . (1) We note that the area-weighted dust density n a (r, z) ( Fig. 2) is defined by an equation of the same form. Figure 6 shows the dust temperature profiles at 100 au used by all our different models. Our disk model is less massive than that used by Heese et al. (2017), leading to smaller opacities, and thus the temperature differences between grains of different sizes is much larger in Figs.5 and 6 than in that previous work. T a Fig. 6: Vertical profiles of temperatures at 100 au. The thick black line is the gas vertical temperature profile in Kelvin. The dotted colored lines are temperatures of each grain population and the thick red line is the area-weighted temperature T a . We see that temperatures are roughly constant for z/H < 1.5 and for z/H > 2.5. Between these two altitudes the temperatures exhibit a strong transition.

Dust extinction, self and mutual shielding
The radiation field is a contribution of both the interstellar radiation field (ISRF) and the stellar radiation field (SRF). Radiation from these two sources encounters molecular or atomic gas and dust grains that affects its propagation through extinction and scattering. The effect of the dust is computed by the 3D Monte-Carlo continuum and line radiative transfer code POLARIS as described in Sec. 2.2.2 ( Fig.3 shows the sampling in the 90 to 400 nm wavelengths range relevant for photoprocesses). Figure  5 shows results at a disk radius of r = 100 au.
In chemical codes, the effect of extinction due to gas, known as self and mutual shielding, is often estimated using empirical attenuation factors as a function of line-of-sight visual extinction (e.g. Lee et al. 1996, for a discussion). The situation for disks is more complex, due to the two sources of UV radiation (stellar neighborhood and ISRF) and varying dust properties. Empirical attenuation factors have been derived for disks (see for example Visser et al. 2009;Heays et al. 2017). They are in general used in a 1+1D approach, where the two sources of radiation are treated independently and scattering is neglected. Furthermore, they depend implicitely on the dust settling that was assumed during their derivation, as well as on the shape of the incident UV field.
For better consistency, we use here a procedure that builds on our knowledge of the dust-attenuated UV field. The extinction produced by the gas requires to know the opacity generated by the molecules through which the radiation travels from outside the disk to the local cell.
We show in Appendix B that a first order approximation for self and mutual shielding is obtained by applying the attenuation of the vertical opacity due to molecules to the radiation field computed with dust only by the POLARIS code: where I d (λ, r, z) is the radiation field given by POLARIS, which explicitly includes the impact of dust scattering, and τ V m is the total opacity due to molecules from the (r, z) point to (r, +∞).

Chemistry simulations
The NAUTILUS gas-grain code ) is used to perform chemistry simulations. NAUTILUS is dedicated to computing chemical evolution in different interstellar environments and uses the rate equation approach (Hasegawa et al. 1992;Hasegawa & Herbst 1993) to calculate the abundance as a function of time. The NAUTILUS gas-grain code is a 3-phase model ) that includes gas-phase chemistry and chemically active grain surfaces and mantles. Thus, aside from the gas-phase chemistry, NAUTILUS covers the physisorption of neutral species on the surface, the diffusion of these species and their reactions, desorption of species from the surface and repopulation of the surface by species from the mantle as the species on the surface evaporate. The latest version of NAUTILUS, the Nautilus Multi Grain Code (NMGC) 1 0D model (Iqbal & Wakelam 2018), can perform simulations using a full grain distribution in size. Each grain size bin is treated independently of the others, and only interact with gas. Compared to the single-grain model (model using a single grain size), a multi-grain model changes chemical rates in two ways. First, accretion rates of species depend on grain population. In a multi-grain model, accretion rates vary on different grain sizes according to their total surface area (Table A.1 shows the relative surface areas as function of grain sizes). Second, dust temperature depends on its size (see Fig. 5). Species have higher reaction rates on grains which are hotter due to higher hopping and desorption rates. In general, grains which are hotter have lower abundances of lighter species such as CN, CH 2 , CO, etc. and more of heavier species such as H 2 O, CH 3 OH, etc. . . (see Iqbal & Wakelam 2018, for details). In our current work, we generalize this multi-grain capability to a 1D situation to represent a protoplanetary disk as described in previous Sections. NMGC computes chemistry in all cells given in cylindrical coordinates (r, z) where we provide as input the gas temperature (Eq. A.5), the local radiation field I d (λ, r, z) and dust temperatures T d (a i , r, z) given by the radiative transfer simulations, and the number density of each dust population n d (a i , r, z) in order to compute the local grain abundances relative to the number density of Hydrogen nuclei n H (r, z): The z direction is treated as in Hersant et al. (2009), but using 64 points.

Photorates
Radiations in the UV range have a critical impact on the disk chemistry since UV photons have the necessary energy to photoionize or photodissociate molecules in the disk. For this reason, characterizing the photorates [s −1 ] at which molecules are dissociated or ionized is of high importance in chemistry models and necessary for a full interpretation of observations. We evaluate the photoprocess rate as in Heays et al. (2017): where σ p (X, λ) is the photoprocess cross-section of species X at wavelength λ and I L (λ, r, z) is the local UV radiation. The index p equals i when it stands for ionization and equal d when it stands for dissociation. We perform our integration over the 91.2 to 400 nm wavelength range. To evaluate the rate at each point in the disk, we thus need the local UV radiation field and the (space-independent) cross-sections. Following Section 2.2.3, the UV field is derived using Eq. 2. The molecular cross-sections used in this study are extracted from the Leiden University website 2 (Heays et al. 2017 The NAUTILUS code initially implemented the H 2 formation through the Langmuir-Hinshelwood (LH) mechanism that considers physisorption on grain surfaces. The LH mechanism is efficient only over a relatively narrow temperature range (5-15 K on flat surfaces Katz et al. 1999;Vidali et al. 2004Vidali et al. , 2005. In the case of multiple grain sizes, dust settling implies that only the small grains remain in the PDR regions of the disk. These small grains illuminated by the UV field from the star can get significantly warmer at the disk surface, and the H 2 formation via the LH mechanism as treated in most astrochemistry codes becomes inefficient, leading to much smaller H 2 formation rate than in the case of a single, larger, equivalent-area, grain size with lower surface temperature. The low formation rate of H 2 leaves a significant amount of atomic Hydrogen that can severely affect the chemical balance in the disk upper layer (see Section 4 for a more detailed discussion).
Observations of H 2 in unshielded regions (Habart et al. 2004(Habart et al. , 2011 where dust can reach high temperatures (≥ 20 K), show that the column density can be significantly higher than what is expected by simple equilibrium rate equation and PDR model predictions. This implies that the formation mechanisms of H 2 can be efficient in a wider domain than predicted by the LH mechanism at equilibrium temperature.
To evaluate the H 2 formation rate, many studies adopt a canonical value for interstellar clouds derived by Jura (1975), which is defined by taking half the rate at which hydrogen atoms stick to the grain surface (e.g. Walsh et al. 2014;Agúndez et al. 2018).
Several other studies proposed more sophisticated mechanisms that theoretically produce comparable H 2 column densities to those derived by observations (e.g. Duley 1996;Cazaux & Tielens 2004;Cuppen & Herbst 2005;Iqbal et al. 2012Iqbal et al. , 2014Thi et al. 2020). We follow here the approach of Bron et al. (2014), hereafter B14, which considers the stochastic fluctuations of both the grain temperature and the H-atom surface population.

Effects of temperature fluctuations
Small grains (≤ 0.01 µm), given their very small heat capacity, have their temperature fluctuating widely even when a relatively small quantity of energy is absorbed. In typical unshielded regions, a small grain can absorb UV photons originated from the interstellar field with a rather high probability, causing a sudden spike in temperature. If the interval of time between UV photon absorptions is sufficiently long, the grain which has undergone transient heating has the time to cool down to temperatures smaller than the equilibrium temperature, giving the opportunity for H 2 to form on the surface. Then, the rapid heating of the grain by the absorption of a UV photon will thermally desorb most species on the surface, releasing H 2 molecules to the gas phase. Bron et al. (2014) have studied the effect of these temperature fluctuations on the formation of H 2 in PDRs and high UV radiation regions using a statistical approach. By numerically solving the master equation (Le Bourlot et al. 2012), they show (Bron et al. 2014, see their Fig.6) that the LH mechanism alone is sufficient to reach the typical H 2 formation rates observed in PDRs when the temperature fluctuations of small grains are accounted for.

Implementing the H 2 formation rates
The investigation of B14 lies in a domain of gas temperature (100 K) and gas density (10 0 to 10 6 cm −3 ) distinct from our study but this difference can easily be corrected as both temperature and gas density only affect the results through the collision rate of atomic Hydrogen with the grains, which scales as n H T g . n H = 10 4 cm 3 n H = 10 6 cm 3 n eq = 10 4 cm 3 n eq = 10 6 cm 3 n eq = 10 7 cm 3 n eq = 10 8 cm 3  Bron et al. (2014) are the blue (10 4 cm −3 ) and purple (10 6 cm −3 ) solid lines while black lines correspond to our analytic fit for four different number densities (10 4 , 10 6 , 10 7 , 10 8 cm −3 ).
Thus, we can estimate the rate for a given density n H and a given temperature T g by extracting or interpolating from B14's data the value corresponding to an equivalent gas density n eq equal to: where s(T g ) is the sticking coefficient at temperature T g as defined by Le Bourlot et al. (2012): with T 2 = 464 K and β = 1.5, which approximately gives the same result as in Sternberg & Dalgarno (1995). We interpolated B14 data as a function of proton density and UV field intensity using an analytic polylogarithmic function. The interpolation results are displayed in Fig.7. The PDR layers of protoplanetary disks are dense, typically above 10 6 − −10 8 cm −3 so that essentially the upper curves in Fig.7 apply to our case.
It should be noted that the B14 calculations assume a powerlaw size distribution from 1 nm to 300 nm composed of only carbonaceous grains. Changing the composition and/or size distribution would affect the precise values of the rates and of the fitted coefficients. However, the asymptotic behaviors of the rates should remain similar. The rate is proportional to the density at low densities, and progressively saturates for larger densities as the fraction of atomic H decreases because the UV flux is insufficient to dissociate H 2 . As a function of the UV field intensity the (unsaturated) rate is expected to decrease as 1/G since the time spent by any grain in the relevant temperature range where the LH formation is efficient is inversely proportional to the UV photon flux.

Detailed description of models
We implement a set of twelve different models (see Table 3). All models share a common set of disk structure and physical parameters (see Fig. 2) as described in Table 1, but differ for a limited number of assumptions on grain sizes and temperatures, UV field or H 2 formation mechanism.
The set is divided into three groups, i.e. the single-grain models, the multi-grain models and an intermediate category of models with single grains locally varying to mimic the grain distribution used in the multi-grain model. We describe their specific characteristics in detail in Section 3.1.1, 3.1.2 and 3.1.3. As mentioned in Section 2.2.1, the excess of UV radiation exhibited by typical TTauri stars, given the strong wavelength-dependence of photoprocesses, inevitably impacts the chemistry of the disk. In order to have a broad view on how our chemical models react to different intensities of UV flux, we adopt two regimes of radiation. One regime uses the spectrum of a relatively high-UV emitting TTauri star and another uses the spectrum of a relatively low UV emitting one, hereafter called HUV and LUV, respectively. For the HUV regime we use the spectrum of TW Hya (Heays et al. 2017) while for the LUV we use the one of DM Tau. Both FUV spectra are provided in the CTTS FUV spectra database 3 (France et al. 2014). The blue solid line in Fig. B.3 exhibits the spectrum of TW Hya as provided by France et al. (2014) before it is attenuated by the disk.
The chemical abundances in the disk are calculated using the time-dependent equation embedded in NAUTILUS. We use atomic initial abundances described in Table B.1 and all the chemical output are extracted after an evolution time of 5.10 6 yrs. We note that using molecular initial abundances would give almost the same final results for disks of this age (see Wakelam et al. 2019). All models incorporate the same chemical network and chemical processes described in Section 2.3. We ignore X-rays because of the large cumulative column densities found around the midplane, the main interest of this paper. Xrays would mostly affect the disk surface layers. We run all models using a standard cosmic-ray ionization rate of 1.3×10 −17 s −1 in the entire disk (see Cleeves et al. 2015, for a detailed discussion on the ionization rate).

Single-grain models
We use single-grain models as a comparison to our more sophisticated multi-grain models described in Section 3.1.2. These models are treated as in Wakelam et al. (2016), in particular with respect to the UV field penetration and photoprocesses (see Hersant et al. 2009). We use the same gas density distribution and same overall gas-to-dust ratio (and hence the same dust mass) as in the multi-grain models. All single-grain models have a single dust size of 0.1 µm. The grains are settled toward the midplane assuming a dust scale height as defined in Equation A.20. However, 0.1 µm grains exhibit little settling, as it can be seen from Model 1: xUV-LH-T g The dust temperature (T d ) profile follows that of the gas (T g ) which is given by the solid black line in Fig.  6), T d is equal to T g everywhere in the disk. x means either H for high or L for low, hence Model 1 is either HUV-LH-T g or LUV-LH-T g , depending on whether the stellar UV field is high (HUV) or low (LUV). The differences are detailed in Section 3.1.
Model 2: xUV-LH-T a Model 2 is similar to Model 1, except for the dust temperature. In Model 2, the grain temperature is given by Eq. 1 and shown in Fig. 6 (red solid line). The temperature T d is hence the area-weighted temperature from the grain size distribution of the multi-grain models and it is called T a . As T a is significantly higher than T g , we expect depletion to occur at lower elevations than in Model 1.
Model 3: xUV-B14-T a In Model 3, the prescription of Bron et al. (2014) as described in Section 2.5 is used to form H 2 . The expected effect is a higher production of H 2 in the PDRs region of the disk than in the other models.

Multi-grain models
Multi-grain models, as opposed to single-grain ones, include the full grain distribution detailed in Table A.1 where each grain population has a specific temperature (Fig. 6), a specific settling factor and is chemically active. For all multi-grain models, we use a size range of [5 mm, 1 mm] with a size exponent d = −3.5 (Eq. A.16) and divide the grain size interval into 16 logarithmically distributed subintervals using Eq. A.19. An overview of the grain size interval with their respective relative dust mass and surface area is given in Table A Model 5: M-xUV-B14 We add the H 2 formation mechanisms prescribed in B14 (see Section 2.5). These models, with grain growth, dust settling, one temperature per grain size and H 2 formation derived from Bron et al. (2014) are the most complex models we compute.

Intermediate models
These models use locally a single grain size mimicking the grain growth and dust settling of the multi-grain models (see Appendix. C.1 for details). At each point, the single grain size, which is chosen in order to obtain the same grain mass and area, are summed up over the grain size distribution of the multi-grain models. The purpose of these models is to better illustrate the impact of a dust temperature spread onto the chemistry, therefore we only run these models at radii 50, 100 and 200 au.
Model 6: Set-HUV-y-T a There are two models. The dust temperature is T a while the H 2 formation (y) either follows the LH mechanism or the B14 approximation. Table 3 summarizes the various models. Table 2 lists the different figures in the text in relation with their respective models. A few figures relevant of the HUV cases are presented in Appendix C, while all figures for the LUV cases appear in Appendix D.  HUV 8,9,11,12,14,15,17,18,19, 20, C.2 multi, LUV D.1, D.4 Set, HUV 8, 12, 15 , C.3, C.4

Results
We present here some general issues about the UV field, the dust temperature and the H 2 formation. A more thorough discussion of the chemistry is done in the next section.

Impact of UV field
A comparison between Tables 4 and D.1 shows that, for a given model, the integrated column densities (or molecular surface densities) have about the same order of magnitude whether the flux is high (HUV) or low (LUV). Unlike in the PDR layer, the midplane chemistry is not dependent on the UV flux regime. The dust opacity is such that the UV flux is sufficiently attenuated in the HUV and LUV models for the chemistry to be similar in both cases. D.3, and D.4 show that the chemistry is nearly the same for all models from the midplane to roughly 1.5 H. We observe more H 2 , CO, CS and CN above 3 scale heights in the LUV models than in the HUV ones because a high UV flux implies more efficient photodissociation. Near the midplane, the UV field is similarly attenuated and the radiation field intensities are about the same in all models. As a consequence, the only parameter that impacts the chemistry is the dust temperature T d and this applies both to the single-grain models and the multi-grain models.

Impact of dust temperature on H 2 formation and chemistry
A first effect of the adopted grain temperature is related to the formation of H 2 and the remaining amount of H, which is very important for the gas-phase chemistry. Adopting the formalism of Bron et al. (2014) to form H 2 produces, as expected, more H 2 . In these cases, Figures C.1 and C.2 show a shift in the H/H 2 transition toward the disk surface which implies less gas-phase H and more H 2 in the upper layers. For single-grain models, the gas temperature we take is relatively low. Table 4 shows that when T d = T g , the remaining amount of H is still reasonably low (intermediate), contrary to the case where T d = T a (HUV-LH-T a ) which exhibits the highest amount of H, together with the multi-grain model M-HUV-LH. Using only the LH mechanism to form H 2 on grains with high temperatures reduces significantly the amount of H 2 above two scale heights. This is particularly true for the M-HUV-LH Model, where the grain temperature depends on the grain size: the reservoir of "cold" grains is not large enough to allow for an efficient H 2 formation on grains in the upper disk layer.
On the contrary, near the midplane, in most models the grain temperature is always low enough for the LH mechanism to efficiently form H 2 and regulate the atomic H abundance.
Notes. HUV stands for high UV field and LUV for low UV field. LH is for the classical traitement of the Langmuir-Hinshelwood mechanism and B14 when the prescription of Bron et al. (2014) is used instead. T g means that the dust temperature equals to that of the gas, T a is the weightedarea dust temperature as defined by Eq. 1. T i is the temperature of i th grain population. M stands for multigrain models and Set for dust models mimicking grain growth and dust settling similar to multi-grain models. As an example, the model M-HUV-B14 is a multi-grain model computed using the B14 method with a high UV field.

Discussion
We first discuss (Sec.4.1) the vertical variations of abundances for CO, CN and CS. We then investigate the C,N and O reservoirs and surface chemistry (Sec.4.2-4.3) before discussing the formation of water and complex organic molecules (Sec.4.4-4.5). We also discuss some possible implications of our results on planet formation and embryos compositions. We focus the discussion that follows (except for Section 4.1) on the HUV case (corresponding to the UV flux received by the TW Hya disk) because we have seen that most of the analysis of HUV models can be considered valid for LUV models. Figures Figure 8 shows the molecular surface densities of a few popular species detected in TTauri disks. Table 4 and D.1 present the column densities at a radius of 100 au for the most abundant species such as CO, CN and CS. These tables reveal that there is no specific trend for a given model, although there can be variations of up to 2.5 orders of magnitude in the column densities. An analysis of Figures 9, C.1 and C.2, which provide the vertical profile of molecular densities at 100 au, is required to understand the origin of these variations. For the column density of CO it shall be noted that 1) CO is formed in the gas phase and 2) is sensitive to photodissociation. The model that produces the most CO is thus LUV-B14-T a because this model combines both a low UV flux, a large column density of H 2 that shields CO from the UV and a high grain temperature that prevents CO from being adsorbed. The second model producing the most CO is M-LUV-B14: the CO column density is slightly lower than in LUV-B14-T a because in multigrain models a fraction of the grains is significantly settled, involving a more effective UV penetration. In that sense, it is not surprising that M-HUV-LH is the model that produces the smallest column density: this model combines a high UV flux, more penetration due to settling and a small production of H 2 which decreases the shielding.
For the column densities of CS and CN, Table 4 shows that using predicted molecular surface densities is not enough to derive general trends. Nevertheless, we can note that M-HUV-B14 is the model that produces the most CS and CN. On the other hand, HUV-B14-T a produces the smallest column density of CS while HUV-LH-T a is the model that produces the smallest column density of CN. This suggests routes of formation which do not simply depend on the grain temperature.
In the case of the single-grain models (see Figs. C.1 and D.3), all models using the same grain temperature (T d = T g or T a ) have the same vertical profile at altitudes z < 1.5H regardless of the flux regime. Similarly, since all multi-grain models use exactly the same grain temperature and size distribution, they all exhibit the same vertical density profile at altitudes z < 1.5H (for a given species) regardless of the incident flux and of the assumed H 2 formation mechanism (see for example, the CO abundance profiles for the two multi-grain models in Figure C.2, middle panels of columns b-c). Accordingly, in this relatively high mass disk, differences in column densities (Tab. 4) result from what happens above z = 1.5H.
Finally, concerning the intermediate model, Set-HUV-B14-T a , Tab. 4 shows that the column densities are rather similar to those found for model HUV-LH-T a , with the exception of CS where the column density is closer to that found for HUV-LH-T g . However, these values are integrated column densities and Figs. C.3 and C.4 show that the vertical profiles are more complex than the comparison based on surface densities may suggest. ] of main molecules at 100 au for HUV models. Last three columns summarize the main properties of the model with respect to the H 2 formation and grains (% of grains, in sites, with a temperature above 20 K between 0 and 1.5H and 1.5 and 2.5H).

Vertical variations
CO Gas-phase CO abundance is both influenced by grain temperature and by abundance of H 2 . For multi-grain models, the flux penetration is larger than in single-grain models, so that CO column densities appear to be globally smaller than in singlegrain models (Fig. 8). In M-HUV-LH, the CO abundance remains low above z ∼ 3H (see Fig. C.2). On the other hand, M-HUV-B14 produces more H 2 allowing CO to survive at higher altitudes. The resulting CO column density is of the same order than in HUV-LH-T a and HUV-B14-T a . Although this has little impact on the total column density, the abundance of CO at z ≤ 1H is much higher in multi-grain models than in singlegrain ones. In single-grain models, the single grain temperature becomes low enough for CO to freeze out efficiently while in Article number, page 11 of 39 A&A proofs: manuscript no. article multi-grain models there is always a fraction of abundant small grains whose temperature remains high enough to prevent CO from being depleted as efficiently as in single-grain models.
CN and Nitrogen bearing Species For all models, the main reactions that create CN in the upper layers at 100 au are bimolecular reactions in the gas-phase: N + C 2 → C + CN (69%) and N + CH → H + CN (22%). Note that CH is mostly destroyed by the reaction H + CH → C + H 2 . Therefore, more atomic H in the upper layers implies less CH available to form CN. The effect is clearly visible for the single-grain models (Fig.C.1) where HUV-B14-T a produces more CN above 3.5H because H 2 is more abundant compared to HUV-LH-T a . The effect is less obvious in the multi-grain models (C.2). Below z = 1.5 H, all the single-grain models exhibit roughly the same abundance, including HUV-LH-T g , which suggests a low sensitivity of CN due to the dust temperature in this range of values.
CS The main formation process of CS in the upper layers is a sequence of a bimolecular reaction and a dissociative recombination reaction: H 2 + CS + → H + HCS + and HCS + + e − → H + CS. The sequence shows that the production of CS is dependent on H 2 . As a consequence, forcing the formation of H 2 leads to the formation of CS, providing that the ionization fraction is the same (as it is the case here).
For single-grain models, the grain temperature is higher in HUV-LH-T a than in HUV-LH-T g . The production of H 2 being more efficient in the latter model at high elevation, the peak of CS in HUV-LH-T g (Fig. C.1) is broader (∼ 2 to 3 H). In the meantime, CS also tends to undergo sticking processes on surfaces higher in the disk, competing with its production in the gas-phase. This explains the dramatic drop of CS density below 2 H. As regards HUV-LH-T a , the production of H 2 starts to be efficient at a lower scale height (∼ 2 versus 3 H). Then, CS abundance dramatically increases at ∼ 1.8 H (because of the steep temperature slope) before rapidly decreasing as temperature is low enough for CS to deplete on the grains, resulting in a very narrow peak. In HUV-B14-T a , the production of H 2 is efficient even in the upper layers and CS starts to form higher in the disk, explaining why the CS abundance is two orders of magnitude larger than in HUV-LH-T g and HUV-LH-T a at 4 H (Fig.C.1).
For multi-grain models, we find a larger CS column density but the location of the peak of abundance is lower in altitude ( Fig. 9 (c)) because of the relatively low visual extinction. Just as for the single-grain models, the surface density is larger in the case of M-HUV-B14 because the formation of H 2 is more efficient and the quantity of H lower. Below z = 1.5 H, the difference in abundance is strictly dependent on T d and as expected, the hot-grain models (HUV-LH-T a and HUV-B14-T a ) produce more CS in the gas-phase than the cold-grain model HUV-LH-T g , whereas all multi-grain models exhibit the same vertical profile.
Concerning show that above the altitude of z ∼ 2H, these models are, as expected, very similar to what is observed for models HUV-LH-T a because they share the same grain temperature for a number of grain sites which is of the same order. Beyond z ∼ 2H, the local gas densities for CO, CN and CS are closer to what is found for the multi-grain models. This is because in the intermediate models the number of grain sites around the midplane is lower than in the case of HUV-LH-T a models and equal to those of multi-grain models. Moreover, the amount of UV should be locally higher than in single grain models because of dust settling. Note that in all cases, at the midplane, the abundance of gas species is always very low, well below observational levels.

Element reservoirs: C, N, O in the cold midplane
The relative importance of reservoirs appearing roughly constant along the distance from the star, we decide to study the reservoirs at 100 au. All species that have been investigated represent at least 0.5 % of the total quantity of carbon, nitrogen or oxygen. As previously seen, the chemistry below z = 1.5 H mostly depends on the grain temperature. Therefore, it turns out to be sufficient to discuss HUV-LH-T g , HUV-LH-T a , M-HUV-LH and Set-HUV-LH-T a only. Figure 12 (top) shows that frozen CO 2 is the main carrier of carbon, both in HUV-LH-T a and M-HUV-LH (and Set-HUV-LH-T a ) with 58.1% and 51.3% of elemental carbon, respectively, while the quantity of CO 2 is negligible in comparison to the other species in HUV-LH-T g . The dust temperature appears too high in these two models to allow for a large quantity of hydrogen to successively hydrogenate the frozen atomic oxygen to form water. Instead, the atomic oxygen easily diffuses on the surface and the relatively high temperature allows to overcome the activation barrier to rapidly form CO 2 by reacting with the frozen CO. In HUV-LH-T g , on the other hand, Fig. 12 (top) shows that the main carrier of carbon is CH 4 with 49.1%. In this colder dust temperature model, CO is initially formed in the gas before being efficiently and rapidly adsorbed onto the grain surfaces, just like in hot-grain models. However, once on the surface, CO is channeled to s-CH 3 OH through hydrogenation sequences. The photodissociation of CH 3 OH then leads to CH 4 formation. Hence, in the midplane, disks with colder grains will produce more CH 4 while warmer disks will produce more CO 2 . The other main carrier is the complex organic species CH 3 OH, which holds about 10% of carbon in M-HUV-LH and HUV-LH-T g . Finally, it is important to mention that CO contributes as a reservoir of C for less than 0.5%, in all models.

Carbon bearing species
Oxygen bearing species CO 2 ice is the main carrier of oxygen in HUV-LH-T a (and Set-HUV-LH-T a ) and M-HUV-LH with 82.0% and 72.5%, respectively (Fig. 12, bottom). The temperature of the grains in these models is too high (T a = 24.6 K) to allow for a large quantity of atomic hydrogen to remain on the grains and hydrogenate the atomic oxygen to form water. As described above, atomic oxygen easily diffuses and also combines with CO to form CO 2 . In HUV-LH-T g , H 2 O is by far the largest carrier of oxygen with 76.4%. Disks with colder grains will produce more H 2 O while disks with warmer grains will produce more CO 2 . CH 3 OH is, again, an important carrier and holds about 10% of oxygen in HUV-LH-T g . As for C, CO holds less than 0.5% of Oxygen, in all models.
Nitrogen bearing species Nitrogen is mostly in the form of N 2 in M-HUV-LH (67.8%) and in HUV-LH-T g (32.4%) (Fig. 12,  middle). Moreover, in almost all models, the gas-phase abundance of N 2 is noticeably larger than the one of the other species. The reason for that is a combined effect of the faster conversion of atomic nitrogen into N 2 in the gas-phase compared to the depletion of N and the high grain temperature that prevents N 2 from depleting onto the surfaces. Indeed, The binding energy of N 2 (E b (N 2 ) = 1100 K on amorphous water ice surface) is smaller than most of the other molecules and slightly smaller than that of CO (Minissale et al. 2016;Wakelam et al. 2017) so that N 2 is retained more efficiently in the gas. We note that, given the low dust temperature, the gas-phase abundance of N 2 in HUV-LH-T g is the smallest of all models. In HUV-LH-T a and in Set-HUV-LH-T a , more N 2 remains in the gas-phase because of the large dust temperatures. On the grains, s-NH 3 is the main carrier of nitrogen with 42.7% of the elemental nitrogen. Therefore N 2 is expected to be the carrier of nitrogen in the disk midplane but remains more in the gas-phase in hot-grain models. The other main carriers of nitrogen are HCN (29.2% in HUV-LH-T g , 12.3% in HUV-LH-T a ) and CH 2 NH 2 (22.9% in HUV-LH-T a , 8.3% in HUV-LH-T g ).

Surface chemistry in the midplane
One important conclusion of Section 4.2 is that disks with colder grains will produce more s-CH 4 and s-H 2 O in the midplane while disks with warmer grains will produce more s-CO 2 . This clearly demonstrates the major role of the grain temperature on the surface chemistry in the plane. We investigate here more deeply the midplane surface chemistry of a few key molecules for the same four models i.e. HUV-LH-T g , HUV-LH-T a , M-HUV-LH, and Set-HUV-LH-T a . Another important conclusion is that our intermediate models, Set-HUV-LH-T a and Set-HUV-B14-T a , are not able to reproduce the complex behavior of the gas-grain chemistry observed with the multi-grain model, since it happens to be closer to the single grain model HUV-LH-T a in terms of main reservoirs in the midplane. Small differences occur for relatively minor constituents (s-CO, s-H 2 CO, s-CH 3 OH, s-HCN for example) because of the different effective dust area in the two cases (Single vs Set), and hence different timescales for conversion on grain surfaces. Figures 13 and 14 present maps of number density for s-H 2 , s-CO, s-CS and s-CN locked on grain surfaces in the case of single-grain and multi-grain models, respectively. Figure 15 presents the surface abundance per total atomic hydrogen for various molecules as a function of the grain radius a and their temperature in the midplane of the disk. The crosses, the square markers, and the round markers stand for the multi-grain model at 30 au, 100 au and 200 au respectively. The triangle pointing upward represents the surface abundance in the single-grain model HUV-LH-T g and the triangle pointing downward represents the abundance in the single-grain model HUV-LH-T a , both at 100 au.
4.3.1. Carbon and oxygen-bearing molecules: CO, CO 2 and CH 4 CO is formed in the gas-phase prior to being accreted on small grains during the first 10 6 yrs of the simulation. However, because of its small binding energy (E b (CO) = 1300 K Wakelam et al. 2017), it is rapidly thermally desorbed from the hot (T g > 20 K) smaller grains and gets accreted by large cold grains (T g < 20 K), where it remains locked until the simulation is stopped. Considering this mechanism, CO is distributed according to the grain population and Fig. 15 (top left panel) shows two regimes of abundance, accounting for the desorption barrier, where CO ice tends to stay locked on grains of size 1 µm at all radii. For the single-grain models, the same mechanism is at play. However, rather than involving a distribution of s-CO on various grain populations, this simply results in a small abundance in the hot-grain model HUV-LH-T a and a high abundance in the coldgrain model HUV-LH-T g .
As for CO 2 , its binding energy is relatively high (E b (CO 2 ) = 2600 K Wakelam et al. 2017). As a consequence the abundance of s-CO 2 depends very little on the dust temperature and simply depends on the dust surface area instead. s-CO 2 is efficiently formed through s-O + s-CO → s-CO 2 (E A = 1000 K) (7) s-OH + s-CO → s-CO 2 + s-H.
Therefore, s-CO 2 is more abundant on the smaller grains because their temperature allows to overcome the activation barrier of reaction 7 and the drop in abundance on large grains at 200 au (and further) originates from the grain temperature getting too low for reaction 7 to be activated. Therefore, the outer disk forms fewer CO 2 ice. In cold-grain models, s-CH 4 is the largest carrier of carbon (see Section 4.2). As expected, Figure 15 (Bottom-left) shows that s-CH 4 is more abundant on cold grains (T g < 20 K). s-CH 4 is formed via the hydrogenation of s-CH 3 . The latter is formed through the following two sequences: s-CO + hν → s-C + s-O (9) s-H 2 + s-C → s-CH 2 (10) s-H 2 + s-CH 2 → s-CH 3 + s-H and s-CH 3 OH + hν → s-OH + s-CH 3 .
In view of the above, s-CH 3 OH is formed on the grains prior to being photodissociated (if the computation domain is sufficiently extended). Then the released atomic carbon (Reaction 9) and s-CH 3 (Reaction 12) lead to s-CH 4 formation. Consequently, the formation of s-CH 4 is more efficient on large cold grains since 1) CO must be adsorbed and 2) hydrogenation is needed to create s-CH 3 OH. s-CH 4 ice is thus found to be more abundant on large grains and in the outer disk.

Sulfur-bearing species: H 2 S
In the gas phase, H 2 S has been identified in dense cloud cores, cometary comae and Phuong et al. (2018) reported the first detection of H 2 S in the cold and dense ring about the TTauri star GG Tau A. The low H 2 S column densities observed in GG Tau A and in dense molecular clouds are assumed to be the result of a strong sulfur depletion (Hudson & Gerakines 2018;Phuong et al. 2018).
We investigate here the formation routes of S-H bonds and more specifically the simple S-bearing molecule Hydrogen sulfide, H 2 S, in the midplane icy mantles of our multi-grain models. Figure 15 (second row, right column) shows the abundance of s-H 2 S as a function of the grain size. There is a strong anticorrelation with grain temperatures. The smallest s-H 2 S abundances are located on the hottest grains (T d ∼ 27 K) while the largest abundances are found on large colder grains (T d 15 K). This anti-correlation is observed at all radii. The abundances on the large grains (> 10µm) are independent of the radius. For the small grains, however, the inner disk exhibits larger abundances than the outer disk because of higher collision rates in the inner regions. The strong anti-correlation with temperature is directly related to the main formation pathways of H 2 S that implies hydrogenation sequences: We note that the channels to s-CH 3 SH also involve successive hydrogenations (of s-CS, mainly). Reaction 13 occurs on grains of all sizes but is more effective on the large cold grains (T d 15 K) and Reaction 14 is only found to be effective on the colder grains.

Nitrogen-bearing species: NH 3
The main nitrogen-bearing icy species in the midplane are s-NH 3 , s-HCN, s-N 2 and s-CN at all radii. In the multi-grain models, as seen in Fig. 15, the abundance of s-NH 3 is inversely proportional to the grain size until ≈ 1µm where the abundance stops decreasing. The abundances of s-NH 3 are roughly the same at 100 and 200 au but not at 30 au.
The main pathway to s-NH 3 originates from the accretion of N at the surface of the grains and its successive hydrogenations to form s-NH 3 (see Aikawa et al. 2015;Eistrup et al. 2018;Ruaud & Gorti 2019). Such hydrogenation reactions are efficient on cold surfaces (T d 15 K), which explains why the abundance of s-NH 3 does not keep dropping when grain temperatures are low. As for the warmer small grain surfaces, their high quantity implies a high collision rate with NH 3 which stays on their surface given its binding energy (E b (NH 3 ) = 5500 K). As the disk evolves, a large fraction of nitrogen-bearing species are trans-formed into complex molecules and other routes become effective such as the destructive hydrogenation of s-NH 2 CO:

Water
The formation of water is key to explain the evolution of protoplanetary disks and formation of comets, both decisive elements to understand the delivery of water to planetary surfaces. that water on grain surfaces cannot exist inside the snowline (T d 150 K), which represents radii smaller than a few to ten au around T Tauri stars. The innermost radius computed by our models is 30 au, which is beyond the radial ice line of water. Therefore, we do not discuss the radial location of the water ice line.

Gas-phase
Three main parameters control the abundance of water in gasphase, the photodesorption by FUV photons and cosmic rays, the dust temperature and the chemical reactions in the gas-phase, in particular: Thus H 2 and OH + are the main precursors for H 2 O. Whether H 2 or OH + is the limiting factor is determined by this other sequence: These two sequences show that H 2 O formation depends on H 2 and OH + but the latter also depends on H 2 . Therefore, this is the H 2 abundance that sets the gas-phase formation of water and accounts for the difference between the models (not cosmic-rays or OH + ). However, the gas-phase formation of H 2 O only represents a small fraction of all H 2 O molecules formed as most H 2 O are formed on grain surfaces. The abundance of gas-phase water is thus mainly governed by both the formation rate of icy water and by the photodesorption rate. Thermal-desorption is ineffective in most parts of the disk because the binding energy of water is high. The formation rate depends on the grain temperature while the photodesorption rate depends on the FUV flux and cosmicrays. The analysis of the results herein (Figs. 16 and 17) shows that the dust temperature plays the major role in the production of gas-phase water. Single-grain models The vertical snowline forms a clear boundary between high and low gas-phase H 2 O abundance at z/r ∼ 0.15 (Fig. 16) and we see that HUV-B14-T a produces the most distinct boundary of all three models. The snowline corresponds roughly to the abrupt drop in the UV flux as shown in Fig. 2 (bottom-right panel). The cold-grain model HUV-LH-T g produces substantially more water in the upper layers than the other single-grain models, because low dust temperature increase the hydrogenation rates on the grain surfaces. The photodesorption rates are approximately equal in the three models (same flux) but since HUV-LH-T g produces a larger water abundance on grain surfaces the final gas-phase abundance is much larger. It results in a column density of water about a factor 10 larger in HUV-LH-T g than in the other two models. On the other hand, the B14's prescription increases only marginally the column density of gas-phase water.
Multi-grain models The vertical snowline is far less apparent (Fig. 17, top row, middle and right column). Both models have approximately the same distribution of abundance. because both have the same dust temperature distribution and same UV. We notice, however, that B14's prescription in M-HUV-B14 produces marginally more water in the upper layer of the disk than in M-HUV-LH, due to the gas-phase formation of water via the two sequences presented above. Fig. 15: Grain surface abundance per total atomic hydrogen and per grain size bin for various molecules at the final stage of integration (i.e. 5 10 6 yrs) as a function of the grain radius a in relation to the grain temperatures. Cross markers stand for the multi-grain model at 30 au, the square markers stand for the multi-grain model at 100 au and the round markers stand for 200 au. The triangles pointing upward represent the surface abundance in the single-grain model HUV-LH-T g (grains of size 0.1 µm) and the triangle pointing downward represents the abundance in the single-grain model HUV-LH-T a , both at 100 au. Abundances being per grain size bin, it is the sum of the abundances in the multi-grain models that should be compared to the abundances in the single grain case.

Water ice
Single-grain models The vertical snowline location is clearly defined (see Fig. 16, bottom row), although we note the presence of a notable band around z = 3 H (z = 2 H) in HUV-LH-T g (resp. HUV-LH-T a ). The reason is the following activated hydrogenation reaction: Indeed, the band illustrates the location in the disk where this reaction is activated. Below the band, the reaction is not activated and water ice is mostly photodesorbed while above the band the photodesorption rates competes with the formation rate and the abundance drops. In HUV-LH-T a (Fig. 16, middle column), Reaction 22 is activated at smaller altitudes in the disk since the dust temperature is globally higher, resulting in a peak in abundance at around 2 scale heights. In contrast, HUV-B14-T a does not exhibit such a peak near z = 2 H since a large fraction of atomic hydrogen on surfaces is used to create molecular hydrogen. This dramatically decreases the quantity of available H atoms necessary to form water.
Multi-grain models The same mechanisms are at work (see middle and right columns of Fig.17). Yet, the dust temperature distribution involves that we find both hot (up to ∼ 95 K) and cold (as low as ∼ 7 K) grains nearly everywhere in the disk. Consequently, Reaction 22 is activated at lower altitudes and the formation of s-H 2 O through hydrogenation can occur at higher altitudes on cold medium-and large-sized grains so that the snowline is located at a higher altitude.
Around the midplane In all models, water is found in the form of ice around the midplane, where the desorption by UV flux and cosmic-rays is limited and where hydrogenation is effective. For single-grain models, the width of the region where water is on the icy mantles is larger in HUV-LH-T g than in the two others because, again, the dust temperature is smaller. As presented in Fig. 15 (bottom right panel), the abundance is roughly inversely proportional to the grain size at 30 au. The reason is the combination of the large binding energy of water and the high density of small grains. At 100 and 200 au, however, the abundance on large grains (size > 0.1 µm) is larger because the grain temperature is low enough to allow for effective hydrogenation.

COMs
Complex organic molecules (COMs), supposedly constituting the crucial bond between simple ISM molecules and prebiotic chemistry, are hard to detect due to their low abundance and weak emission lines originating from their molecular complexity. This is especially true for objects with small spatial extensions like disks. Simple organic molecules have been regularly detected over the last decades such as H 2 CO (Dutrey et al. 1997;Aikawa et al. 2003;Öberg et al. 2011), HC 3 N (Chapillon et al. 2012, CH 3 CN (Öberg et al. 2015) and CH 3 OH (Walsh et al. 2016). We investigate here the effectiveness of our models to synthesize COMs from surfaces or gas-phase when given sufficient time (5 × 10 6 yrs).

Methyl Cyanide (CH 3 CN)
C-N bonded species are very important for the synthesis of Glycine (Goldman et al. 2010). Öberg et al. (2015) have observed methyl cyanide (CH 3 CN) in the disk around the Herbig Ae star MWC 480, a star exhibiting a stronger UV flux and a warmer disk than a TTauri star.
CH 3 CN to HCN abundance ratio To facilitate comparisons, we discuss the abundance of CH 3 CN relative to HCN, as practiced in cometary studies (Mumma & Charnley 2011;Öberg et al. 2015). The model that produces the largest CH 3 CN/HCN ratio is the single-grain model HUV-LH-T g with ∼ 14% (vertically integrated) at 100, 120 and 140 au while the other singlegrain models produce nearly more than ∼ 1% at these radii. The outer regions (200-250 au) produce a larger ratio (up to 9%) than the inner disk regions in the warm grain single-models, whereas the peak is in the region around 100 au in HUV-LH-T g . In multi-grain models, the gas-phase abundance ratio never exceeds 1%. However, it never goes below 0.1%, implying a rather uniform production of CH 3 CN along the radius as compared to the single-grain models. M-HUV-B14 produces slightly more CH 3 CN than M-HUV-LH. The model HUV-LH-T g (with the coldest grains overall) is the one that globally produces the largest CH 3 CN/HCN ratio.
CH 3 CN in gas-phase In the upper regions, the main known pathways to form gas-phase CH 3 CN are the reaction of HCN with CH + 3 and the photodesorption of icy CH 3 CN (thermal desorption is not effective because of the large binding energy of CH 3 CN). Results show that the reaction between HCN and CH + 3 is the most effective route in M-HUV-B14 while only photodesorption is effective in M-HUV-LH. The reason is straightforward since the production of CH + 3 requires the presence of molecular hydrogen: and H 2 is much more abundant in M-HUV-B14 in the upper layers. At 100 au, the abundance of CH + 3 in M-HUV-LH is very low (∼ 10 −20 ) as compared to the abundance in M-HUV-B14 (∼ 10 −10 ). It follows that, as seen in Fig. 18, M-HUV-B14 produces much more CH 3 CN than in M-HUV-LH. Therefore, the presence of H 2 appears to be decisive for the production of gasphase CH 3 CN in the upper regions of the disk. However, this difference between the two models does not significantly impact the total column density of CH 3 CN as most of it is produced near the midplane. Figure 20 (bottom) shows the abundance of s-CH 3 CN ices in the midplane relative to the grain size as a function of the radius. s-CH 3 CN mostly stays on warm small grains while there is a large drop in abundance around grain populations of size ∼ 0.317 µm. The results show that the two main reactions that lead to s-CH 3 CN are the adsorption of gas-phase CH 3 CN onto the surfaces and the following surface reaction:  . We see that the abundance dramatically drops above ∼ 2 scale heights in M-HUV-LH while M-HUV-B14 keeps a rather high abundance at high altitudes. This is due to the available quantity of gas-phase H 2 in M-HUV-B14, which enhances the production of CH + 3 .
small grains that implies large collision rates with 1) gas-phase produced CH 3 CN and 2) with CH 3 that, despite a low binding energy, will briefly stick and quickly react with s-CN (see reaction 24). The second reason is the high binding energy of CH 3 CN (E b (CH 3 CN) = 4680 K Wakelam et al. 2017) and of the precursors (s-CN) so that they stay on small grain surfaces. Moreover, the temperature of small grains is sufficiently high to allow precursors to diffuse on the surface and increase reaction rates. The last reason is that even if s-CH 3 CN ices can be formed on large grain surfaces, photodesorption rates are large enough for CH 3 CN to be desorbed from large grains and seed small grains as the disk evolves.

Methanol and formaldehyde
CH 3 OH in gas-phase For multi-grain models, Fig. 19 shows that the vertical profile of gas-phase CH 3 OH is divided into two layers separated by an abrupt drop in abundance at ∼ 2 H, which corresponds roughly to the dust temperature transition (Fig. 6).
There is a peak in gas-phase CH 3 OH at ∼ 1.8 H in both models. This corresponds to the location where dust temperature is low enough for adsorbed CO to be hydrogenated and where the UV flux is high enough to efficiently photodissociate a significant fraction of newly formed methanol. The lower layer (< 2 H) contains most of CH 3 OH column density. In the upper layer, M-HUV-B14 produces more CH 3 OH than M-HUV-LH.
CH 3 OH to H 2 CO abundance ratio In multi-grain models, the ratio CH 3 OH/H 2 CO increases by a factor of ∼ 2 to 3 in the outermost region (≥ 200 au) between 10 6 and 5 × 10 6 years while it remains stable in the inner region (< 100 au). In M-HUV-B14, after 10 6 years, the averaged ratio in the region 80 -120 au is equal to 1.3. This appears to be consistent with the observed ratio in TW Hya (Walsh et al. 2016;Carney et al. 2019). In M-HUV-LH, the averaged ratio is equal to 2.4. The difference is due to a more efficient pathway to form H 2 CO in the upper disk regions of M-HUV-B14. Gas-phase H 2 CO is efficiently formed through the following sequence: This shows that molecular hydrogen is a key element in the formation process of gas-phase formaldehyde and M-HUV-B14 provides much more H 2 than M-HUV-LH.
Grain size dependence of surface chemistry for CH 3 OH Figure 20 gives a set of maps of ice abundance in the midplane with grain sizes as a function of the computed radii. Top left panel is the grain temperature which allows to analyze the surface abundance in view of the surface temperature. We have same results for both multi-grain models. On icy mantles, both s-H 2 CO and s-CH 3 OH are formed via the same hydrogenation sequence: However, we note that H 2 CO is mainly formed in the gas-phase.
As seen in Fig. 20 (top right), CO tends to be efficiently adsorbed onto cold surfaces (T d 17 K). The map of H 2 CO (bottom left) grossly follows the same distribution although the trend is globally less pronounced. CH 3 OH, on the other hand, exhibits a totally different distribution and tend to concentrate on small grain surfaces (bottom right). The main reason comes from the fact that the binding energy increases as the molecules become more complex (e.g., E b (CO) ∼1200 K, E b (H 2 CO) ∼3200 K, E b (CH 3 OH) ∼5500 K Minissale et al. 2016;Noble et al. 2012;Collings et al. 2004) and that the hydrogenation of s-H 2 CO (which leads to s-CH 3 OH) has an activation barrier and cannot efficiently occur on cold grain surfaces. s-CH 3 CN t=5 × 10 6 yrs 10 12 10 11 10 10 10 9 10 8 ab(s-CH 3 CN) Fig. 20: multi-grain models: Maps of ice abundance relative to H in the midplane, with grain sizes as a function of the radius. Top left panel is the distribution in grain temperature. We see that the hottest region is located at around 60 au from the star on grains of size 0.148 µm.

Impact for planet formation
Planet formation is strongly linked to grain growth and dust disk evolution. Our study reveals some interesting trends which should influence it.
For a T Tauri disk of a few Myr, radial changes of the millimeter spectral index (e.g. Birnstiel et al. 2010) reveal that, after settling, large grains have drifted toward the central star. Typically, grains larger than a few 0.1-1 mm are located inside a radius of 50-80 au. Beyond this, dust disks have a spectral in-dex which is similar to that observed in the ISM, suggesting a similar size distribution, even if settling still occurs with larger grains residing onto the midplane. As a consequence, toward the midplane, one expects to observe different icy layers on grains in the inner (large cold grains) and outer disk (warmer smaller grains). We have seen in the previous section that s-CH 4 and s-H 2 O are preferentially located on large grains while s-CO 2 remains on smaller ones (Fig. 15). This implies also a different chemical coupling between the gas and solid phase. Moreover, we have seen that complex molecules tend to stick and remain on small warmer grains (≤ ∼ 0.1 µm), which are less likely to drift toward the inner disk. The outer disk (r ≥ 80 au) is therefore assumed to support rich chemistry on small grains which are probable precursors to cometary nuclei. On the other hand, as larger grains drift efficiently toward the inner region, their surface temperature is likely to vary during the journey which can allow for the growing chemical complexity to occur until the less refractory molecules become vaporized when sufficiently close to the star. Therefore, rich chemistry of the inner disk region can be assumed to come from the gas-phase which is later accreted by planetary embryos, while rich chemistry of the outer disk region occurs on the small grain surfaces and then seeds the inner region through planet-comet collisions in later phases of the disk.
Grain-grain collisions may modify this picture. The characteristic timescale (the time for a large grain to encounter small grains of equivalent area) around 100 au is of order 10 4 -10 5 years. However, while large grains may accrete the smaller ones, and inherit part of their surface composition while growing, small grains are re-created by disruptive events that may severely affect their surface layer and ice mantles. The outcome of the surface state during such grain-grain collisions remains to be studied.
Nevertheless, in all cases, the grain size distribution, because it results in a spread of dust temperature, allows for a faster depletion of C and O, that are converted to s-CO 2 on the colder grains in the disk midplane, but also a more active surface chemistry leading to formation of COMs on the warmer grains in the upper layer. Such a behavior cannot be reproduced by a single grain temperature chemical model, because the formation processes are highly non-linear, since the mobility and evaporation efficiency depend exponentially on dust temperature.
The range in dust temperatures also spreads out spatially the snowlines that become fuzzy. The spatial extent of these "snow bands" will depend on the specific dust properties and disk mass. The temperature spread is larger for less massive disks.
As disks evolve, the amount of dust (and gas) decreases until they reach the optically thin, debris disk phase. A less massive dust disk necessarily implies a warmer dust disk, which should affect the composition of the icy layers at the grain surface in evolved disks. Recent surveys of 5-40 Myr old disks suggest the existence of an intermediate phase (at least for Herbig Ae stars) between debris and protoplanetary disks (e.g. Kóspál et al. 2013;Péricaud et al. 2017). The so-called hybrid disks share dust properties similar to those observed around debris disks, while the amount of gas is still significant, and a larger than standard gasto-dust-ratio is observed. For instance, in HD141569, Di Folco et al. (2020) observed a CO surface density which is only a factor 10 lower than in a typical protoplanetary HAe disk, whereas the amount of millimeter grains is a factor 50 times smaller than in typical young disks. The reduced amount of grains, and in some cases at least possible changes in the grain size distribution (e.g. Hughes et al. 2018), should affect the grain surface chemistry during the late phases of disk evolution. More simulations of older, gas-rich disks with a small amount of dust and a larger gas-to-dust ratio are needed to quantify the impact on chemistry and planet formation.

Summary
. We presented here a new chemical model of protoplanetary disks that consistently accounts for the impact of grain sizes on the dust temperatures and UV penetration. Representative results were obtained for a typical disk around T Tauri stars, using parameters derived from the Flying Saucer observations . We coupled the chemical code NAUTILUS and radiative transfer code POLARIS to use a consistent distribution of grain sizes and temperatures, accounting for size-dependent dust settling. We implemented in NAUTILUS a new method to account for self and mutual shielding of molecules in photo processes, method that accounts for the frequency dependence of photorates. The outputs of multi-grain models are compared to those of single grain size, using two different assumptions for the single-grain temperature: equal to the (assumed) gas temperature, or to the area-weighted mean temperature of the multigrain models.
-We identify the formation of H 2 through Langmuir-Hinshelwood mechanism as inefficient for a grain size distribution, and the amount of atomic H as critical for the chemistry of the upper layers (2-4 H). The parametric method of Bron et al. (2014), that accounts for temperature fluctuations for small grains, appears more appropriate to produce enough H 2 . -Comparisons between the models reveal differences of more than 1-2 orders of magnitudes for some abundant gas-phase species. Most of these differences arise from above 1.5H, despite the larger UV penetration in multi-grain models. -The location of the H 2 O vertical snowline is partly regulated by the UV penetration but also by the amount of H 2 . In the upper disk layers, the amount of H 2 also regulates the formation of H 2 CO, although this molecule is mostly formed, like CH 3 OH, on grain surfaces. -CH 3 CN, like CH 3 OH to a lesser extent, is preferentially formed on warmer grains that allow for a larger diffusivity of precursors on the grain surface. -At 100 au, in the midplane, disks with colder dust produce more H 2 O and CH 4 while those with warmer dust produce more CO 2 on grain surfaces. -As soon as CO is locked onto large (hence colder) grains (a > 1µm), it remains trapped until the end of the simulated evolution. -The spread of temperatures available with grain size distributions simultaneously provides surface chemistry for cold and warm grains that cannot be mimicked by a unique dust temperature. In particular, it allows depletion of C and O to proceed efficiently (by conversion to CO 2 or CH 4 and H 2 O on colder grains), while the COMs production is boosted on the warmer grains. -Using a single grain size that locally mimicks grain growth and dust settling fails to reproduce the complex gas-grain chemistry resulting from grain temperature depending on grain size. -Another direct consequence of having different dust temperatures at the same location is a spatial spread-out of the snowlines. High angular observations would be required to test this prediction.
Our results clearly show that the chemistry of protoplanetary disks cannot simply be represented by a single grain size and temperature.
Finally, our disk model assumes a relatively high dust mass (i.e. an optically thick dust disk). An optically thinner dust disk would impact the UV penetration and dust temperature, changing the chemistry at play even on the midplane. This remains to be investigated to determine how it can affect the grain surfaces and then the composition of planetesimals and embryos. The scale height H g (r) is given by and using the definition of the temperature law from Eq.A.4 where µ m is the mean molecular weight, and m H the Hydrogen mass. The temperature and velocity being power laws, the scale height H(r) also follows a power law: Vertical temperature gradient: With a vertical temperature gradient, the vertical gas profile deviates from a Gaussian profile. In Nautilus, we compute the gas density for each vertical cell z i at a given r by solving out locally the equation as in Hersant et al. (2009): Appendix A.2: Dust distribution and properties Appendix A.2.1: Dust material We assume that the grains are made of an homogeneous mixture of a population of silicate grains (62.5%) and graphite grains (37.5%). This composition is used both for radiative transfer and chemical modeling. This leads to a material specific mass density for dust grains in all this study equal to ρ m = 2.5 g.cm −3 . The grain temperatures are defined for each grain sizes using radiative transfer simulations as described in Section 2.2.2.
We consider a size dependent grain population. Overall, this grain population is derived assuming no radial drift, but accounting for dust settling as a function of height. The surface density of a grain population of size a is thus given by where σ d,0 (a) is the surface density of grains of size a at reference radius. Eq. A.13 implies that the total dust surface density is: which we associate to the overall (vertically integrated, equivalently disk averaged) dust-to-gas mass ratio ζ through Σ d (r) = ζΣ g (r). The mass fraction of grains of size [a, a + da] is given by where the size distribution in number dn(a) is assumed to follow a power law of exponent d: where C is a normalization constant. C is related to the dust mass by simply integrating over the size as follows: where ρ m is the specific mass density defined in the previous section. Using Eq.A.15 and the definition of ζ yields For radiative transfer and chemistry simulations the grain size distribution is discretized into several logarithmically distributed subranges with relative dust masses of the i th grain size interval: The relative dust mass per bin and the relative area per bin (which also quantifies the number of chemical sites) are given in Table A.1.

Appendix A.2.3: Dust settling
To account for size-dependent dust settling, we assume that the vertical profile for a given grain size has a Gaussian shape: (A.20) Because dust grains do not react to pressure forces in the same way as gas, the dust scale height will differ from that of the gas. Following Boehler et al. (2013), we define the settling factor, that can be related to the dust properties through the stopping time. The stopping time of a dust particle τ s is the characteristic time for a dust particle initially at rest to reach the local gas velocity. We use the dimensionless stopping time T s (r, z) defined as the product of τ s and the Keplerian angular momentum Ω(r) which is a way to compare the stopping time to the dynamical time in the disk. Garaud & Lin (2004) showed that dust friction with gas in protoplanetary disks is well described by the Epstein regime, in which the dimensionless stopping time is given by: The stopping time depends linearly on the grain size. Using Eq. A.9 we can write: The settling factor is given by the approximation of Dong et al. (2015): where α is the α viscosity coefficient, is the dimensionless stopping time in the midplane, and S c is the Schmidt number. The Schmidt number is a dimensionless number defined by the ratio between the turbulent viscosity ν t and the turbulent diffusion D t :   Figure A.2 shows the settling factor at the reference radius as defined in Dong et al. (2015) compared to the asymptotic behavior used by Boehler et al. (2013) and the values from the numerical simulations of Fromang & Nelson (2009). Table A.2 provides the total dust cross-section at 50, 100, and 200 au from the star, defined by the term n d (r, z = 0, a)πa 2 .

Appendix B: Self and Mutual Shielding
The local UV flux in the disk is not a simple sum of the attenuated stellar and interstellar fields because of scattering by dust grains. Scattering can be significant near the disk midplane, where the optical depths are large. However, molecular shield-  ing is most important only in the upper layers, when dust optical depth remains low while molecular optical depth can be significant.
Thus, at first order, to evaluate the impact of self-shielding, we neglect scattering in our analysis by writing in a simple form the local flux I in the cell of coordinates (r, z) as follows: where I IS RF is the incident flux from the ISRF, I * (r, z) is the spatially diluted, unattenuated, incident flux from the central star, τ V is the vertical opacity generated by the disk matter above the local point (r, z) and τ R is the radial opacity generated by the matter between the central star and the local point. As the opacities are a contribution of both dust and gas, we write: where m stands for molecules and d for dust. So the local expression becomes: The local flux weighted by the dust only I d can also be approximated by: POLARIS treats consistently scattering by dust and this simple expression was checked to be a reasonable approximation in the upper layers where molecular shielding is important.
At a given wavelength, the molecular opacity is defined as follows: where X includes all molecules for which we have cross-sections (including H 2 O), n(X, r, z) is the number density [cm −3 ] of species X at coordinates (r, z) and σ(X, λ) is the sum of the absorption, dissociation and ionization cross-sections of species X at wavelength λ. Depending on whether the radiation is from the interstellar medium or from the central star, the line l over which we integrate is either the altitude z above the cell of coordinates (r, z) or the radial distance √ r 2 + z 2 between the central star and the considered cell, respectively. Figure B.2 shows the vertical opacities generated both by the dust (dotted lines) and the molecules (solid lines). The opacities are given for different altitudes at 100 au. H 2 , CO and N 2 are the main contributors to the line opacity.
As shown in Figs. B.2, at all wavelengths where H 2 contributes to the attenuation (λ 115 nm), its opacity is in all parts of the disk always substantially larger than the dust opacity.
While calculating the vertical molecular opacity τ V m is trivial, characterizing the radial molecular opacity τ R m is more difficult as it couples different radii and thus their chemical evolution together. Given the significant computing time involved in the multi-grain model, we show below that we can reasonably simplify the problem through a simple approximation.
Self or mutual shielding is only relevant for high molecular opacities, because frequencies where molecular shielding occurs cover a limited fraction of the UV domain. In regard to Figure B.1, we introduce the ratio between the radial and vertical molecular optical depths as follows: and Eq. B.3 can be rewritten as: Strictly speaking, RV is position and frequency dependent, since the radial and vertical distribution of molecules differ. The geometry shown in Figure B.1 indicates that RV should be ≥ 1 and we approximate Eq. B.7 by: generally the physico-chemistry in the upper disk atmosphere. Studying the disk upper layer is out of the scope of our paper (and would require the use of a PDR code). In fact, our approximation becomes reasonable as soon as dust attenuation is significant, a situation which happens at z/H < 3.5 − 3.8 in our dust disk model (which has a relatively high dust mass), and does not affect our results on the molecular layer and below. In summary, everywhere in the disk, self and mutual shielding are considered by applying the attenuation by the vertical opacity due to molecules to the radiation field computed with dust only by the POLARIS code: where I d (λ, r, z) explicitly includes the impact of dust scattering that is handled by POLARIS. It is important to mention that compared to empirical, precomputed factors as a function of visual extinction and H 2 column densities (or also CO and N 2 for self-shielding), this approach has the advantage to take consistently the dust extinction into account for the specific dust properties and distribution.
Consequently, we use Eq. B.9 to compute the photorates and the vertical molecular opacities are computed using Eq. B.5 in all our simulations. Figure B.3 gives the flux as a function of the wavelength in the midplane and in the upper atmosphere (4 H) at 100 au for the disk whose parameters are given in Table 1. It shows a comparison between the flux as given by the radiative transfer simulation (dotted lines) and the flux as given by Eq. 2 which takes into account the molecular contribution to the opacity (solid line). The attenuation by spectral lines, especially those of H 2 , is clearly visible. The contribution of the molecules to the opacity occurs in a limited range of the UV domain. Although our chemical code (Section 2.3) computes other molecular species that contribute to the UV opacity, their abundance relative to H 2 remains small and limit their contribution.