© S. Gavino et al. 2021

1 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 so-called protostellar phase) small dust grains are essentially dynamically well coupled to the gas. When grain growth occurs (up to at least centimeter-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 CO snowline, which is typically located at a radius about 20 au in a T Tauri disk, the midplane is very dense and cold (≤20K). This area is essentially shielded from the stellar radiation, with 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 The Atacama Large Millimeter/submillimeter Array (ALMA), several high angular resolution and sensitive observationsof disks orbiting either low- or intermediate-mass young stars (T Tauri 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), also revealing that the CO gas is located in the intermediate warm molecular layer. More recently, the ALMA observations of the edge-on T Tauri disk called the Flying Saucer (Dutrey et al. 2017) directly show the molecular layer at intermediate scales (1–2 H) with CS; the denser gas tracer is located at about 1H below the CO emission (1–3 H), which better traces the whole envelope of the molecular layer.

In the last ten years, several astrochemical models have emerged that attempt to incorporate the observed complexity of disk physics to provide more accurate molecular abundance 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. This step includes three material phases: the gas-phase and two phases for the gas-grain chemistry with grains modeled by an icy surface layer and an active mantle (Ruaud & Gorti 2019; Wakelam et al. 2019). Thermochemical models that calculate the density and thermal structures in a self-consistent manner have also been developed (e.g. Hollenbach & Gorti 2009). Recently, Ruaud & Gorti (2019) have also coupled such a thermochemical model (Gorti et al. 2011) to a three-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 on 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, where 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 in which 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 that 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 an important CO depletion in the disk of CQ Tau, whose temperature is well above the CO snowline temperature. These authors suggested that CO may remain trapped onto larger grains that are cold enough to prevent thermal CO desorption. Heese et al. (2017) investigated the dust temperature variations with grain size and position (radius and altitude above the midplane) in a typical T Tauri 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 dependence questions the applicability of elaborate, self-consistent, and thermochemical models that rely on a single-grain size and temperature.

We explore in this work the impact of the variations of the grain temperature with their size onto the chemistry of a representative T Tauri 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 we explore the main consequences of dust grain temperature dependence with size.

For this purpose, we 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 three-phase version, which also takes into account a grain size distribution (Iqbal & Wakelam 2018). We present the two codes and how we parameterize the gas and dust disk in Sect. 2. Section 3 deals with the disk model description and global results. We discuss in more detail the results in Sect. 4 and we state our conclusions in Sect. 5.

2 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 that we independently assume a gas and dust disk structure, a gas temperature, a dust distribution, and settling. Moreover, the coupling between the two codes requires a new method to calculate the dust extinction for NAUTILUS, 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.

2.1 Disk model

Our disk structure is derived from a simple disk model in which a 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^−dda with dn(a), the number of grains 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 that have grain sizes ranging from a_min = 5 nm to a_max = 1 mm. The fractionof mass and area per bin are given in Table A.1.

Because our goal is 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) to avoid chemical effects that could be due to gas alone. We selected values derived from the observations of the Flying Saucer T Tauri disk (see Dutrey et al. 2017), whose edge-on orientation facilitates 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).

Fig. 1 Multi-grain model: coupling scheme between NAUTILUS and POLARIS.

2.2 Radiation fields and dust temperature

In this new approach, we use the dust radiative transfer code to calculate the uv field and dust temperature for each of the 16 grain sizes, in each cell of the disk (in radius and altitude), throughout the disk. The dust temperatures are then used to compute the chemistry with NAUTILUS.

However, the uv field derivedabove only accounts for dust extinction and scattering. A proper accounting for self and mutual shielding of molecules and atoms is required to handle the photo-chemical processes.

2.2.1 Radiation sources

As the central heating source we assume a low-mass, pre-main-sequence star with a mass of M_⋆ = 0.58M_⊙, which radiates as a blackbody with an effective temperature of T_⋆,eff = 3900K 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, T Tauri 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; spectral lines such as Lyman α have 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, because photo processes are 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 and chemical model the sum of the stellar blackbody spectrum and the uv excess typically found in T Tauri stars. Further details about the radiation fields are given in Sect. 3.1.

2.2.2 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. 3). 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. 3, middle panel). The dust temperature increases with grain size for sizes in the range 0.007 μm-0.070 μm, then decreases for sizes in the range 0.32–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 cross section 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 dependon the size and location of the grains. To compare our set of models with those consisting of a single-grain size (see Sect. 3.1.1), it is also convenient to introduce an area-weighted temperature defined as $$T_{\text{a}}(r,z)=\frac{{\sum }_i{a}_i^2{T}_{\text{d}}(a_i,r,z)n_{\text{d}}(a_i,r,z)}{{\sum }_i{a}_i^2{n}_{\text{d}}(a_i,r,z)}.$$(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 5 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. 3 and 5 than in that previous work.

Fig. 2 Physical structure of the multi-grain models incorporating a full dust distribution. The quantity Tg is the gas temperature, nH 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 togas mass ratio, Ta (Eq. (1)) is the area-weighted dust temperature, and na is the area-weighted grain number density. The quantity AV is the visualextinction counted from disk surface toward the disk midplane.

Fig. 3 Results of the radiative transfer simulation in vertical direction at a disk radius of r = 100 au. Left: dust density distribution for the 16 grain size intervals as described in Appendices A.2.2 and A.2.3. Middle: dust temperature of the 16 grain size intervals. Right: stellar (solid lines) and ISRF (dotted lines).

2.2.3 Dust extinction, self and mutual shielding

The radiation field is a contribution of both the 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 Sect. 2.2.2 (Fig. 6 shows the sampling in the 90–400 nm wavelengthsrange relevant for photo processes). Figure 3 shows results at a disk radius of r = 100au.

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 owing to the two sources of uv radiation (stellar neighborhood and ISRF) and varying dust properties. Empirical attenuation factors have been derived for disks (see, e.g., Visser et al. 2009; Heays et al. 2017). These factors are in general used in a 1 + 1D approach, where the two sources of radiation are treated independently and scattering is neglected. Furthermore, they implicitly depend on the dust settling that was assumed during their derivation, as well as on the shape of the incident uv field.

For better consistency, in this work we use a procedure that builds on our knowledge of the dust-attenuated uv field. The extinction produced by the gas requires knowledge of 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, that is, $$I_L(\lambda ,r,z)={\text{e}}^{-{\tau }_{\text{m}}^V}{I}_{\text{d}}(\lambda ,r,z),$$(2)

where I_d(λ, r, z) is the radiation field given by POLARIS, which explicitly includes the impact of dust scattering, and ${\tau }_{\text{m}}^V$ is the total opacity due to molecules from the (r, z) point to (r, +∞).

Fig. 4 Ratio of the absorption cross sections for wavelengths characteristic for stellar radiation (Cabs,1 μm) and dust emission (Cabs,20 μm).

Fig. 5 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 Ta. 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.

Fig. 6 Absorption (blue dashed line) and scattering cross sections (green dotted line) for the first grain size bin (7 nm). The verticalgray lines denote the wavelengths for which the radiation fields are simulated.

2.3 Chemistry simulations

The NAUTILUS gas-grain code (Ruaud et al. 2016) is used to perform chemistry simulations. The NAUTILUS code 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 three-phase model (Ruaud et al. 2016) 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, 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) 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 interacts 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 in 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. 3). Species have higher reaction rates on hotter grains as a result of higher hopping and desorption rates. In general, hotter grains have lower abundances of lighter species such as CN, CH₂, and CO and have higher abundances of heavier species such as H₂O and CH₃OH (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 the cylindrical coordinates (r, z) for which we provide as input the gas temperature (Eq. (A.5)), local radiation field I_d (λ, r, z) and dust temperatures T_d(a_i, r, z), given by the radiative transfer simulations and number density of each dust population n_d (a_i, r, z) to compute the local grain abundances relative to the number density of hydrogen nuclei n_H (r, z): $$A{B}_{\text{d}}(a_i,r,z)=\frac{n_{\text{d}}(a_i,r,z)}{n_{\text{H}}(r,z)}.$$(3)

The z direction is treated as in Hersant et al. (2009), but using 64 points.

2.4 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⁻¹ ] 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 photo process rate as in Heays et al. (2017) as follows: $$k_{\text{p}}(X,r,z)={\displaystyle {\int }_{\lambda }{\sigma }_{\text{p}}}(X,\lambda )I_L(\lambda ,r,z)\text{d}\lambda ,$$(4)

where σ_p(X, λ) is the photo process 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–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 Sect. 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¹ (Heays et al. 2017). These cross sections are collected from experimental and theoretical studies. Eighty molecules are treated in this way. The photorates are recomputed at each explicit time step of the chemical evolution. For all other gas constituents(atoms or molecules unavailable in Heays et al. 2017), we use approximate shielding factors from Visser et al. (2009).

2.5 Molecular hydrogen formation at the disk surface

2.5.1 The atomic H problem

The NAUTILUS code initially implemented the H₂ formation through the Langmuir-Hinshelwood (LH) mechanism that considers physisorption on grain surfaces. The LH mechanism is only efficient over a relatively narrow temperature range (5–15 K on flat surfaces Katz et al. 1999; Vidali et al. 2004, 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₂ formation via the LH mechanism as treated in most astrochemistry codes becomes inefficient. This leads to much smaller H₂ formation rate than in the case of a single, larger, equivalent-area grain size with lower surface temperature. The low formation rate of H₂ leaves a significant amount of atomic hydrogen that can severely affect the chemical balance in the disk upper layer (see Sect. 4 for a more detailed discussion).

Observations of H₂ in unshielded regions (Habart et al. 2004, 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₂ can be efficient in a wider domain than predicted by the LH mechanism at equilibrium temperature.

To evaluate the H₂ 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₂ column densities to those derived by observations (e.g., Duley 1996; Cazaux & Tielens 2004; Cuppen & Herbst 2005; Iqbal et al. 2012, 2014; Thi et al. 2020). We follow the approach of Bron et al. (2014), hereafter B14, which considers the stochastic fluctuations of the grain temperature and the H-atom surface population.

2.5.2 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, thereby causing a sudden spike in temperature. If the interval of time between uv photon absorptions is sufficiently long, the grain that has undergone transient heating has the time to cool down to temperatures smaller than the equilibrium temperature, thereby giving the opportunity for H₂ to form on the surface. Then, the rapid heating of the grain by the absorption of a uv photon thermally desorbs most species on the surface, releasing H₂ molecules tothe gas phase.

Bron et al. (2014) studied the effect of these temperature fluctuations on the formation of H₂ 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₂ formation rates observed in PDRs when the temperature fluctuations of small grains are accounted for.

2.5.3 Implementing the H₂ formation rates

The investigation of B14 lies in a domain of gas temperature (100 K) and gas density (10⁰ to 10⁶ cm⁻³) distinct from our study but this difference can easily be corrected because temperature and gas density only affect the results through the collision rate of atomic hydrogen with the grains, which scales as $n_{\text{H}}\sqrt{T_{\text{g}}}$. Thus, we can estimate the rate for a given density n_H and a given temperature T_g by extracting or interpolating from B14 data the value corresponding to an equivalent gas density n_eq equal to $$n_{\text{eq}}=n_{\text{H}}\sqrt{\frac{T_{\text{g}}}{100\text{K}}}\frac{s(T_{\text{g}})}{s(100\text{K})},$$(5)

where s(T_g) is the sticking coefficient at temperature T_g as defined by Le Bourlot et al. (2012), $$s(T)={\left(1+{\left(\frac{T}{T_2}\right)}^{\beta }\right)}^{-1},$$(6)

where T₂ = 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 shown in Fig. 7. The PDR layers of protoplanetary disks are dense, typically above 10⁶ − 10⁸ cm⁻³ so that essentially the upper curves in Fig. 7 apply to our case.

It should be noted that the B14 calculations assume a power-law 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 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₂. 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 at which the LH formation is efficient is inversely proportional to the uv photon flux.

Fig. 7 H2 formation rate as a function of the uv field intensity (in units of the Habing interstellar field) for various proton number density. Data adapted from Fig. 6 of Bron et al. (2014) are shown as the blue (104 cm−3) and purple (106 cm−3) solid lines, while black lines correspond to our analytic fit for four different number densities (104, 106, 107, 108 cm−3).

3 Model parameters and results

3.1 Detailed description of models

We implement a set of 12 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₂ formation mechanism.

The set is divided into three groups, that is, the single-grain models, multi-grain models, and an intermediate category of models with single grainslocally varying to mimic the grain distribution used in the multi-grain model. We describe their specific characteristics in detail in Sects. 3.1.1, 3.1.2, and 3.1.3. As mentioned in Sect. 2.2.1, the excess of uv radiation exhibited by typical T Tauri stars, given the strong wavelength-dependence of photo processes, 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 T Tauri star and a second regime uses the spectrum of a relatively low uv emitting T Tauri star, 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 that of DM Tau. Both far-ultraviolet (FUV) spectra are provided in the CTTS FUV spectra database² (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⁶ yr. 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 Sect. 2.3. We ignore X-rays because of the large cumulative column densities found around the midplane, which is the main focus of this paper. X-rays would mostly affect the disk surface layers. We run all models using a standard cosmic-ray ionization rate of 1.3 × 10⁻¹⁷ s⁻¹ in the entire disk (see Cleeves et al. 2015, for a detailed discussion on the ionization rate).

3.1.1 Single-grain models

We use single-grain models as a comparison to our more sophisticated multi-grain models described in Sect. 3.1.2. These models are treated as in Wakelam et al. (2016), in particular with respect to the uv field penetration and photo processes (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 Eq. (A.20). However, 0.1 μm grains exhibit little settling, as shown in Fig. A.1, and such models are virtually equivalent to non-settled dust models (Wakelam et al. 2016). Abundance maps of simple molecules are presented in Figs. 8 and D.2.

Model 1: xUV-LH-T_g

The dust temperature (T_d) profile follows that of the gas (T_g) (given by the solid black line in Fig. 5); T_d is equal to T_g everywhere in the disk. The quantity 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 Sect. 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. 5 (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 Sect. 2.5 is used to form H₂. The expected effect is a higher production of H₂ in the PDRs region of the disk than in the other models.

3.1.2 Multi-grain models

Multi-grain models, as opposed to single-grain models, include the full grain distribution detailed in Table A.1, where each grain population has a specific temperature (Fig. 5), 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.1 while Fig. A.1 shows the scale height of each 16 grain size population as a function of the radius. The abundance maps of simple molecules are presented in Fig. 9. Multi-grain models are composed of two categories.

Model 4: M-xUV-LH

Either called M-HUV-LH or M-LUV-LH depending on the flux, they use the classical LH mechanisms for H₂ formation without considering small grains temperature fluctuation.

Model 5:M-xUV-B14

We add the H₂ formation mechanisms prescribed in B14 (see Sect. 2.5). These models, with grain growth, dust settling, one temperature per grain size, and H₂ formation derived from Bron et al. (2014), are the most complex models we compute.

3.1.3 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 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₂ 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 for the HUV cases are presented in Appendix C, while all figures for the LUV cases appear in Appendix D.

3.2 Results

In this section, we present some general issues about the uv field, dust temperature, and H₂ formation. The next section has a more thorough discussion of the chemistry.

3.2.1 Impact of uv field

A comparison between Tables 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.

Figures C.1, C.2, 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₂, CO, CS, and CN above three scale heights in the LUV models than in the HUV models 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.

3.2.2 Impact of dust temperature on H₂ formation and chemistry

A first effect of the adopted grain temperature is related to the formation of H₂ 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₂ produces, as expected, more H₂. In these cases, Figs. C.1 and C.2 show a shift in the H/H₂ transition toward the disk surface, which implies less gas-phase H and more H₂ 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 in which 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₂ on grains with high temperatures significantly reduces the amount of H₂ 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₂ 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₂ and regulate the atomic H abundance.

4 Discussion

We first discuss (Sect. 4.1) the vertical variations of abundances for CO, CN, and CS. We then investigate the C, N, and O reservoirs and surface chemistry (Sects. 4.2–4.3) before discussing the formation of water and complex organic molecules (Sects. 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 Sect. 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 relevant to the LUV case (similar to what is expected for the DM Tau disk) are presented inAppendix D. All results are presented at the final stage of time evolution (i.e., 5 × 10⁶ yr).

Fig. 8 Density (cm−3) of H2, CO, CS and CN in the gas phase of the single-grain models in HUV regime. Left column is the HUV-LH-Tg model, the middle column is HUV-LH-Ta, and the rightcolumn is HUV-B14-Ta. The black contours represent the dust temperature (Td = Tg in the left column and Td = Ta in the middle and right columns).

4.1 Vertical distribution of abundant (easily observed) molecules

4.1.1 Column densities

Figure 10 shows the molecular surface densities of a few popular species detected in T Tauri disks. Tables 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 ofup to 2.5 orders of magnitude in the column densities. An analysis of Figs. 11, 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 CO is formed in the gas phase and 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₂ 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 multi-grain models a fraction of the grains is significantly settled, involving a more effective uv penetration. In that sense, it is not surprising that the M-HUV-LH model produces the smallest column density: this model combines a high uv flux, more penetration from settling, and a small production of H₂, 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 the M-HUV-B14 model produces the most CS and CN. On the other hand, HUV-B14-T_a produces the smallest column density of CS while the HUV-LH-T_a model produces the smallest column density of CN. This suggests routes of formation that 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.5 H (for a given species) regardless of the incident flux and of the assumed H₂ formation mechanism (see for example, the CO abundance profiles for the two multi-grain models in Fig. C.2, middle panels of columns b-c). Accordingly, in this relatively high-mass disk, differences in column densities (Table 4) result from what happens above z = 1.5 H.

Finally, concerning the intermediate model, Set-HUV-B14-T_a, Table 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.

Fig. 9 Density (cm−3) of H2, CO, CS, andCN in the gas phase of HUV-LH-Tg (left column) and of the multi-grain models in HUV regime. The black contours represent the dust temperature (Td = Tg in the left column and Td = Ta in the middle and right columns).

4.1.2 Vertical variations

CO

Gas-phase CO abundance is influenced by grain temperature and by abundance of H₂. 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 single-grain models (Fig. 10). In M-HUV-LH, the CO abundance remains low above z ~ 3 H (see Fig. C.2). On the other hand, M-HUV-B14 produces more H₂ allowing CO to survive at higher altitudes. The resulting CO column density is of the same order as 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 higherin multi-grain models than in single-grain ones. In single-grain models, the single-grain temperature becomes low enough for CO to freeze out efficiently, while in 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.

Fig. 10 Surface densities of the single-grain models (solid lines) and multi-grain models (dotted lines) as a function of the radius in high uv flux regime. The crosses are the intermediate model for a selection of radii.

Fig. 11 Vertical cumulative surface density (cm−2) of H2, CO, CS andCN at 100 au from the star in the HUV 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₂ → C + CN (69%) and N+ CH → H + CN (22%). We note that CH is mostly destroyed by the reaction H + CH → C + H₂. 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.5 H because H₂ 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₂ + CS⁺ → H + HCS⁺ and HCS⁺ + e⁻ → H + CS. The sequence shows that the production of CS is dependent on H₂. As a consequence, forcing the formation of H₂ 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. Because the production of H₂ is 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₂ starts to be efficient at a lower scale height (~2 versus 3 H). Then, CS abundance dramaticallyincreases 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₂ 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. 11c) 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₂ is more efficient and the quantity of H lower. Belowz = 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 the intermediate models (Set-HUV-B14-T_a and Set-HUV-LH-T_a), Fig. C.3 and Fig. C.4 show that above the altitude of z ~ 2 H, 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 ~ 2 H, 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 higher locally than in single-grain models because of dust settling. We note that in all cases, at the midplane, the abundance of gas species is always very low, well below observational levels.

4.2 Element reservoirs: C, N, O in the cold midplane

Because the relative importance of reservoirs appears roughly constant along the distance from the star, we decided to study the reservoirs at 100 au. All species that were 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.

Carbon bearing species

Figure 12 (top) shows that frozen CO₂ 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₂ 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 us to overcome the activation barrier to rapidly form CO₂ 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₄ 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₃OH through hydrogenation sequences. The photodissociation of CH₃OH then leads to CH₄ formation. Hence, in the midplane, disks with colder grains produce more CH₄ while warmer disks produce more CO₂. The other main carrier is the complex organic species CH₃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.

Oxygen-bearing species

CO₂ 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₂. In HUV-LH-T_g, H₂O is by far the largest carrier of oxygen with 76.4%. Disks with colder grains produce more H₂O while disks with warmer grains produce more CO₂. CH₃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₂ 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₂ is noticeably larger than that of the other species. The reason for that is a combined effect of the faster conversion of atomic nitrogen into N₂ in the gas phase compared to the depletion of N and the high grain temperature that prevents N₂ from depleting onto the surfaces. The binding energy of N₂ (E_b (N₂) = 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), therefore N₂ is retained more efficiently in the gas. We note that, given the low dust temperature, the gas-phase abundance of N₂ 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₂ remains in the gas phase because of the large dust temperatures. On the grains, s-NH₃ is the main carrier of nitrogen with 42.7% of the elemental nitrogen. Therefore, N₂ is expected to be the carrier of nitrogen in the disk midplane but is more present 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₂NH₂ (22.9% in HUV-LH-T_a, 8.3% in HUV-LH-T_g).

4.3 Surface chemistry in the midplane

One important conclusion of Sect. 4.2 is that disks with colder grains produce more s-CH₄ and s-H₂O in the midplane while disks with warmer grains produce more s-CO₂. This clearly demonstrates the major role of the grain temperature on the surface chemistry in the plane. We investigate 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₂CO, s-CH₃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₂, s-CO, s-CS, and s-CN locked on grain surfaces in the case of single-grain and multi-grain models, respectively. Figure 15presents the surface abundance per total atomic hydrogen for various molecules as a function of the grain radius a and their temperaturein the midplane of the disk. The crosses, square markers, and 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₂ and CH₄

CO is formed in the gas phase prior to being accreted on small grains during the first 10⁶ yr of the simulation. However, because of its small binding energy (E_b (CO) = 1300 K; Wakelam et al. 2017), CO 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 cold-grain model HUV-LH-T_g.

As for CO₂, its binding energy is relatively high (E_b(CO₂) = 2600 K; Wakelam et al. 2017). As a consequence the abundance of s-CO₂ depends very little on the dust temperature and simply depends on the dust surface area instead. s-CO₂ is efficiently formed through $$\begin{array}{lll}\hfill & \hfill & \text{s0O+s0CO}\to {\text{s0CO}}_{\text{2}}({\text{E}}_{\text{A}}\text{=}\text{1000}\text{K})\hfill \\ \hfill & \hfill & \text{s0OH+s0CO}\to {\text{s0CO}}_{\text{2}}\text{+s0H}.\hfill \end{array}$$

Therefore, s-CO₂ is more abundant in the smaller grains because their temperature allows these grains to overcome the activation barrier of reaction 7 and the drop in abundance in large grains at 200 au (and further) originates from the grain temperaturegetting too low for reaction 7 to be activated. Therefore, the outer disk forms less CO₂ ice.

In cold-grain models, s-CH₄ is the largest carrier of carbon (see Sect. 4.2). As expected, Fig. 15 (Bottom-left) shows that s-CH₄ is more abundant in cold grains (T_g < 20 K). s-CH₄ is formed via the hydrogenation of s-CH₃. The latter is formed through the following two sequences: $$\begin{array}{lll}\hfill & \hfill & \text{s0CO+h}\nu \to \text{s0C+s0O}\hfill \\ \hfill & \hfill & {\text{s0H}}_{\text{2}}\text{+s0C}\to {\text{s0CH}}_{\text{2}}\hfill \\ \hfill & \hfill & {\text{s0H}}_{\text{2}}{\text{+s0CH}}_{\text{2}}\to {\text{s0CH}}_{\text{3}}\text{+s0H}\hfill \end{array}$$

and $${\text{s0CH}}_{\text{3}}\text{OH+h}\nu \to {\text{s0OH+s0CH}}_{\text{3}}.$$(12)

In view of the above, s-CH₃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₃ (Reaction 12) lead to s-CH₄ formation. Consequently, the formation of s-CH₄ is more efficient on large cold grains since (1) CO must be adsorbed and (2) hydrogenation is needed to create s-CH₃OH. s-CH₄ ice is thus found to be more abundant on large grains and in the outer disk.

4.3.2 Sulfur-bearing species: H₂S

In the gas phase, H₂S has been identified in dense cloud cores, cometary comae and Phuong et al. (2018) reported the first detection of H₂S in the cold and dense ring about the TTauri star GG Tau A. The low H₂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₂S, in the midplane icy mantles of our multi-grain models. Figure 15 (second row, right column) shows the abundance of s-H₂S as a function of the grain size. There is a strong anti-correlation with grain temperatures. The smallest s-H₂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₂S that implies hydrogenation sequences: $$\begin{array}{lll}\hfill & \hfill & \text{s0S}\stackrel{\text{s0H}}{\to }\text{s0HS}\stackrel{\text{s0H}}{\to }{\text{s0H}}_{\text{2}}\text{S}\hfill \\ \hfill & \hfill & {\text{s0CH}}_{\text{3}}\text{SH}\stackrel{\text{s0H}}{\to }{\text{s0CH}}_{\text{3}}{\text{+s0H}}_{\text{2}}\text{S}.\hfill \end{array}$$

We note that the channels to s-CH₃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.

Fig. 12 Species that contain at least 0.5% of elemental carbon (top), nitrogen (middle), and oxygen (bottom) in the midplane at 100 au: percentage of the carbon/nitrogen/oxygen locked in the species as a function of the gas-phase species abundance.

Fig. 13 Density (cm−3) of H2, CO, CS, and CN in the gas phase of the single-grain models in the HUV regime. The left column shows the HUV-LH-Tg model, middle column shows HUV-LH-Ta and right shows HUV-B14-Ta. The black contours represent the dust temperature (Td = Tg in the left column and Td = Ta in the middle and right columns).

Fig. 14 Density (cm−3) of H2, CO, CS and CN in the gas phase of HUV-LH-Tg (left column) and of the multi-grain models in HUV regime. The black contours represent the dust temperature (Td = Tg in the left column and Td = Ta in the middle and right columns).

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 106 yr) as a functionof the grain radius a in relation to the grain temperatures. The crosses stand for the multi-grain model at 30 au, the squares stand for the multi-grain model at 100 au, and the circles stand for 200 au. The triangles pointing upward represent the surface abundance in the single-grain model HUV-LH-Tg (grains of size 0.1 μm) and the triangle pointing downward represents the abundance in the single-grain model HUV-LH-Ta, both at 100 au. Abundances being per grain size bin, the sum of the abundances in the multi-grain models should be compared to the abundances in the single-grain case.

4.3.3 Nitrogen-bearing species: NH₃

The main nitrogen-bearing icy species in the midplane are s-NH₃, s-HCN, s-N₂, and s-CN at all radii. In the multi-grain models, as seen in Fig. 15, the abundance of s-NH₃ is inversely proportional to the grain size until ≈1 μm, where the abundance stops decreasing. The abundances of s-NH₃ are roughly the same at 100 and 200 au but not at 30 au. The main pathway to s-NH₃ originates from the accretion of N at the surface of the grains and its successive hydrogenations to form s-NH₃ (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₃ 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₃, which stays on their surface given its binding energy (E_b (NH₃) = 5500 K). As the disk evolves, a large fraction of nitrogen-bearing species are transformed into complex molecules and other routes become effective such as the destructive hydrogenation of s-NH₂CO: $${\text{s0NH}}_{\text{2}}\text{CO}\stackrel{\text{s0H}}{\to }{\text{s0NH}}_{\text{3}}\text{+s0CO}.$$(15)

4.4 Water

The formation of water is fundamental to explaining the evolution of protoplanetary disks and the formation of comets, which are both decisive elements to understand the delivery of water to planetary surfaces. Spectral lines of H₂O in the gas phase of inner disk regions have now been widely detected in disks around T Tauri stars using near- and mid-infrared (Carr & Najita 2008; Pontoppidan et al. 2010) or far-infrared observations (Riviere-Marichalar et al. 2012). Cold water in the outer disk regions, on the other hand, is harder to detect as a consequence of low excitation lines.

It is assumed (Dominik et al. 2005; Lecar et al. 2006; Podio et al. 2013; Du & Bergin 2014) that water on grain surfaces cannot exist inside the snowline(T_d ≳ 150 K), which represents radii smaller than a few to ten astronomical units 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.

4.4.1 Gas phase

Three main parameters control the abundance of water in the gas phase, the photodesorption by FUV photons and cosmic rays, the dust temperature, and the chemical reactions in the gas phase, in particular, $$\begin{array}{lll}\hfill & \hfill & {\text{H}}_{\text{2}}{\text{+OH}}^{\text{+}}\to {\text{H+H}}_{\text{2}}{\text{O}}^{\text{+}}\hfill \\ \hfill & \hfill & {\text{H}}_{\text{2}}{\text{O}}^{\text{+}}{\text{+H}}_{\text{2}}\to {\text{H+H}}_{\text{3}}{\text{O}}^{\text{+}}\hfill \\ \hfill & \hfill & {\text{H}}_{\text{3}}{\text{O}}^{\text{+}}{\text{+e}}^{\text{-}}\to {\text{H+H}}_{\text{2}}\text{O}\hfill \end{array}$$

Thus H₂ and OH⁺ are the main precursors for H₂O. Whether H₂ or OH⁺ is the limiting factor is determined by another sequence as follows: $$\begin{array}{lll}\hfill & \hfill & {\text{H}}_{\text{2}}\text{+CR}\to {\text{H}}_{\text{2}}^{\text{+}}{\text{+e}}^{\text{-}}\hfill \\ \hfill & \hfill & {\text{H}}_{\text{2}}^{\text{+}}{\text{+H}}_{\text{2}}\to {\text{H}}_{\text{3}}^{\text{+}}\text{+H}\hfill \\ \hfill & \hfill & {\text{H}}_{\text{3}}^{\text{+}}\text{+O}\to {\text{OH}}^{\text{+}}{\text{+H}}_{\text{2}}\hfill \end{array}$$

These two sequences show that H₂O formation depends on H₂ and OH⁺ but the latter also depends on H₂. Therefore, this is the H₂ 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₂O only represents a small fraction of all H₂O molecules formed as most H₂O are formed on grain surfaces. The abundance of gas-phase water is thus mainly governed by 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 cosmic rays. The analysis of the results herein (Figs. 16 and 17) shows that the dust temperature plays a 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₂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 increases 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 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 they both have the same dust temperature distribution and uv. We note, however, that the B14 prescription in M-HUV-B14 produces marginally more water in the upper layer of the disk than in M-HUV-LH owing to the gas-phase formation of water via the two sequences presented above.

4.4.2 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). This is because of the following activated hydrogenation reaction: $$\text{s0H+s0HOOH}\to {\text{s0H}}_{\text{2}}\text{O+s0OH}({\text{E}}_{\text{A}}\text{=}\text{1400}\text{K}).$$(22)

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 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₂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.

Fig. 16 Water abundance in the gas phase (top) and on the grain surfaces (bottom) of the single-grain models. The left column shows the abundances in HUV-LH-Tg, middle in HUV-LH-Ta, and right in HUV-B14-Ta. In the top row, the solid black lines show the dust temperature isocontours; in the bottom row they denote isocontours of scale heights.

Fig. 17 Water abundance in the gas phase (top) and on the grain surfaces (bottom). The left column shows the abundances in HUV-LH-Tg and the middle and right column in the multi-grain models. In the top row, the solid black lines show the dust temperatureisocontours; in the bottom row they denote isocontours of scale heights.

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.

4.5 COMs

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

4.5.1 Methyl cyanide (CH₃CN)

C–N bondedspecies are very important for the synthesis of glycine (Goldman et al. 2010). Öberg et al. (2015) have observed methyl cyanide (CH₃CN) in the disk around the Herbig Ae star MWC 480, a star exhibiting a stronger uv flux and a warmer disk than a T Tauri star.

CH₃CN to HCN abundance ratio

To facilitate comparisons, we discuss the abundance of CH₃CN relative to HCN, as practiced in cometary studies (Mumma & Charnley 2011; Öberg et al. 2015). The model that produces the largest CH₃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 single-grain 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₃CN along the radius as compared to the single-grain models. M-HUV-B14 produces slightly more CH₃CN than M-HUV-LH. The model HUV-LH-T_g (with the coldest grains overall) globally produces the largest CH₃CN/HCN ratio.

CH₃CN in gas phase

In the upper regions, the main known pathways to form gas-phase CH₃CN are the reaction of HCN with CH${}_3^{+}$ and the photodesorption of icy CH₃CN; thermal desorption is not effective because of the large binding energy of CH₃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 for this is straightforward since the production of CH${}_3^{+}$ requires the presence of molecular hydrogen as follows: $${\text{C}}^{\text{+}}\stackrel{{\text{H}}_{\text{2}}}{\to }{\text{CH}}_{\text{2}}^{\text{+}}\stackrel{{\text{H}}_{\text{2}}}{\to }{\text{H+CH}}_{\text{3}}^{\text{+}}$$(23)

and H₂ 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⁻²⁰) as compared to the abundance in M-HUV-B14 (~10⁻¹⁰). It follows that, as seen in Fig. 18, M-HUV-B14 produces much more CH₃CN than in M-HUV-LH. Therefore, the presence of H₂ appears to be decisive for the production of gas-phase CH₃CN in the upperregions of the disk. However, this difference between the two models does not significantly impact the total column density of CH₃CN as most of it is produced near the midplane.

Fig. 18 Altitude above the midplane as a function of the abundance relative to H of CH3CN at 100 au. The blue line shows the abundance produced by model M-HUV-LH; the red line shows the abundance produced by model M-HUV-B14. 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 H2 in M-HUV-B14, which enhances the production of CH${}_3^{+}$.

CH₃CN ices

Figure 19 (bottom) shows the abundance of s-CH₃CN ices in the midplane relative to the grain size as a function of the radius. s-CH₃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₃CN are the adsorption of gas-phase CH₃CN onto the surfaces and the following surface reaction: $${\text{s0CH}}_{\text{3}}\text{+s0CN}\to {\text{s0CH}}_{\text{3}}\text{CN}.$$(24)

s-CH₃CN is more abundant on small grains for three main reasons. The first reason is the large number density of small grains that implies large collision rates with gas-phase produced CH₃CN and with CH₃ that, despite a low binding energy, briefly sticks and quickly reacts with s-CN (see reaction 24). The second reason is the high binding energy of CH₃CN (E_b (CH₃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₃CN ices can be formed on large grain surfaces, photodesorption rates are large enough for CH₃CN to be desorbed from large grains and seed small grains as the disk evolves.

Fig. 19 Multi-grain models: maps of ice abundance relative to H in the midplane, with grain sizes as a function of the radius. The top left panel is the distribution in grain temperature. The hottest region is located at around 60 au from the star on grains of size 0.148 μm.

4.5.2 Methanol and formaldehyde

CH₃OH in gas-phase

For multi-grain models, Fig. 20 shows that the vertical profile of gas-phase CH₃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. 5). There is a peak in gas-phase CH₃OH at ~ 1.8 H in both models. This corresponds to the location at which the dust temperature is low enough for adsorbed CO to be hydrogenated and at which 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₃OH column density. In the upper layer, M-HUV-B14 produces more CH₃OH than M-HUV-LH.

CH₃OH to H₂CO abundance ratio

In multi-grain models, the ratio CH₃OH/H₂CO increases by a factor of ~2–3 in the outermost region (≥ 200 au) between 10⁶ and 5 × 10⁶ yr, while it remains stable in the inner region (<100 au). In M-HUV-B14, after 10⁶ yr, 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₂CO in the upper disk regions of M-HUV-B14. Gas-phase H₂CO is efficiently formed through the following sequence: $$\begin{array}{lll}\hfill & \hfill & {\text{H}}_{\text{2}}{\text{+CH}}_{\text{3}}^{\text{+}}\to {\text{CH}}_{\text{5}}\text{+}\hfill \\ \hfill & \hfill & {\text{CH}}_{\text{5}}^{\text{+}}{\text{+e}}^{\text{-}}\to {\text{H+H+CH}}_{\text{3}}{\text{/H}}_{\text{2}}{\text{+CH}}_{\text{3}}\hfill \\ \hfill & \hfill & {\text{O+CH}}_{\text{3}}\to {\text{H+H}}_{\text{2}}\text{CO}.\hfill \end{array}$$

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₂ than M-HUV-LH.

Fig. 20 Altitude above the midplane as a function of the abundance relative to H of CH3OH at 100 au. The blue line shows the abundance produced by model M-HUV-LH; the red line shows the abundance produced by model M-HUV-B14. The abundance dramatically drops at ~ 2 scale heights in both models, although M-HUV-B14 maintains a larger abundance than M-HUV-LH at high altitudes.

Grain size dependence of surface chemistry for CH₃OH

Figure 19 gives a set of maps of ice abundance in the midplane with grain sizes as a function of the computed radii. The top left panel is the grain temperature that allows us to analyze the surface abundance in view of the surface temperature. We have same results for both multi-grain models. On icy mantles, s-H₂CO and s-CH₃OH are formed via the same hydrogenation sequence as follows: $$\text{s0CO}\stackrel{\text{H}}{\to }\text{s0HCO}\stackrel{\text{H}}{\to }{\text{s0H}}_{\text{2}}\text{CO}\stackrel{\text{H}}{\to }{\text{s0CH}}_{\text{3}}\text{O}\stackrel{\text{H}}{\to }{\text{s0CH}}_{\text{3}}\text{OH}.$$(28)

However, wenote that H₂CO is mainly formed in the gas phase. As seen in Fig. 19 (top right), CO tends to be efficiently adsorbed onto cold surfaces (T_d≾ 17 K). The map of H₂CO (bottom left) grossly follows the same distribution although the trend is globally less pronounced. CH₃OH, on the otherhand, exhibits a totally different distribution and tends to concentrate on small grain surfaces (bottom right). The main reasoncomes from the fact that the binding energy increases as the molecules become more complex (e.g., E_b (CO) ~ 1200 K, E_b (H₂CO) ~ 3200 K, E_b (CH₃OH) ~ 5500 K Minissale et al. 2016; Noble et al. 2012; Collings et al. 2004) and that the hydrogenation of s-H₂CO (which leads to s-CH₃OH) has an activation barrier and cannot efficiently occur on cold grain surfaces.

4.6 Impact for planet formation

Planet formation is strongly linked to grain growth and dust disk evolution. Our study reveals some interesting trends that should influence planetary formation.

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 index that 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, we expect 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₄ and s-H₂O are preferentially located on large grains while s-CO₂ remains on smaller grains (Fig. 15). This also implies 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 on the order of 10⁴ –10⁵ yr. However, while large grains may accrete the smaller grains, and inherit part of their surface composition while growing, small grains are recreated by disruptive events that may severely affect their surface layer and ice mantles. The outcome of the surface state duringsuch 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, which are converted to s-CO₂ 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 nonlinear, since the mobility and evaporation efficiency depend exponentially on dust temperature.

The range in dust temperatures also spatially spreads out the snowlines that become fuzzy. The spatial extent of these “snow bands” depends 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 gas-to-dust-ratio is observed. For instance, in HD141569, Di Folco et al. (2020) observed a CO surface density that 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.

5 Summary

In this work we presented 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 (Dutrey et al. 2017). 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, which 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 multi-grain models.

• We identify the formation of H₂ through the LH 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), which accounts for temperature fluctuations for small grains, appears more appropriate to produce enough H₂:

• Comparisons between the models reveal differences of more than one to two orders of magnitudes for some abundant gas-phase species. Most of these differences arise from above 1.5 H, despite the larger uv penetration in multi-grain models;

• The location of the H₂O vertical snowline is partly regulated by the uv penetration but also by the amount of H₂. In the upper disk layers, the amount of H₂ also regulates the formation of H₂CO, althoughthis molecule is mostly formed, like CH₃OH, on grain surfaces;

• CH₃CN, like CH₃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₂O and CH₄, while those with warmer dust produce more CO₂ 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, this allows depletion of C and O to proceed efficiently (by conversion to CO₂ or CH₄ and H₂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.

Acknowledgements

We thank E. Bron for fruitful discussions about the H₂ formation in PDRs and M. Ruaud for a detailed reading of the article. We also thank an anonymous referee who helped to improve the paper. This work was supported by the “Programme National de Physique Chimie du Milieu Interstellaire” (PCMI) from INSU/CNRS. J.K. acknowledges support from the DFG grants WO 857/13-1, WO 857/15-1, and WO 857/17-1.

Appendix A Formal description of the model

Unless otherwise specified, we use subscript g for gas quantities and subscript d for dust. Subscript 0 is used for reference values.

A.1 Gas disk structure

We describe the gas temperature and disk structure in this section. We neglect gas pressure and disk self-gravity, and thus assume that the disk is in Keplerian rotation around the host star as follows: $$v_g(r)=\sqrt{\frac{G{M}_{*}}{r}}⟺{\Omega }_g(r)=\sqrt{\frac{G{M}_{*}}{r^3}},$$(A.1)

where r is the radius.

The gas surface density Σ_g(r) also follows a power law $${\Sigma }_g(r)={\Sigma }_{g,0}{(\frac{r}{R_0})}^{-p}.$$(A.2)

and, by integration, the gas surface density Σ_g,0 at the reference radius R₀ for p≠2 is related to the gas disk mass by $${\Sigma }_{g,0}=\frac{M_g{R}_0^{-p}\left(2-p\right)}{2\pi \left(R_{\text{out}}^{2-p}-R_{\text{in}}^{2-p}\right)},$$(A.3)

where R_in and R_out are the inner and outer radii of the disk.

The gas temperature is not self-consistently derived from the dust temperature distribution but imposed by analytical laws. The kinetic temperatureT_k in the midplane is given by $$T_k(r,z=0)=T_{\text{mid},0}{(\frac{r}{R_0})}^{-q}.$$(A.4)

Following Dartois et al. (2003), we allow for a warmer disk atmosphere using the following formulation of Williams & Best (2014): $$T_{\text{g}}(r,z)=T_{\text{mid}}(r)+(T_{\text{atm}}(r)-T_{\text{mid}}(r))\sin {(\frac{\pi z}{2z_{\text{atm}}})}^{2\sigma },$$(A.5)

where σ is the stiffness of the vertical temperature profile and z_atm is the altitude at the upper boundary of our disk model (4 scale heights H), meaning that the temperature above 4 H is constant. The atmosphere temperature is also given by a power law with the same exponent as that of the midplane temperature as follows: $$T_{\text{atm}}(r)=T_{\text{atm},0}{(\frac{r}{R_0})}^{-q}.$$(A.6)

Isothermal case

If the gas temperature is vertically isothermal, we can derive the vertical density structure by assuming hydrostatic equilibrium at the midplane temperature. Consequently the gas vertical mass distribution follows a Gaussian profile $${\rho }_g(r,z)={\rho }_{g,mid}(r)\exp \left(-\frac{z^2}{2H_{\text{g}}\text{}{\left(r\right)}^2}\right),$$(A.7)

where ρ_g,mid(r) is the density in the midplane, that is related to the surface density by $${\rho }_{g,mid}(r)=\frac{{\Sigma }_g(r)}{\sqrt{2\pi }{H}_g(r)}.$$(A.8)

Fig. A.1 Scale heights of the grain distribution for the 16 different grain sizes as a function of the radius. The grain averaged size interval ranges from 0.007 to 651 μm. Small grains scale heights (size < 0.317 μm) practically follow the gas scale height. For larger grains the settling factor (Eq. A.21) becomes significant and the dust scale height decreases.

The scale height H_g(r) is given by $$H_g(r)=\frac{C_s(r)}{{\Omega }_g(r)},$$(A.9)

and using the definition of the temperature law from Eq.A.4 $$H_{g,0}=\frac{C_s\left(R_0\right)}{\Omega \left(R_0\right)}=\sqrt{\frac{k_{\text{b}}{T}_{\text{mid},\text{0}}{R}_0^3}{{\mu }_{\text{m}}{m}_{H}G{M}_{*}}},$$(A.10)

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 $$H_g(r)=H_{g,0}{(\frac{r}{R_0})}^{\frac{3}{2}-\frac{q}{2}}.$$(A.11)

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), $$ln({\rho }_g(z_i)=ln({\rho }_g(z_i-1))-({\Omega }^2\frac{\mu m_H}{k_B{T}_g(z_i)}+(ln(T_g(z_i))-ln(T_g(z_{i-i}))).$$(A.12)

A.2 Dust distribution and properties

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 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⁻³. The grain temperatures are defined for each grain sizes using radiative transfer simulations as described in Section 2.2.2.

A.2.2 Size distribution

We considera 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 $${\sigma }_d(r,a)={\sigma }_{d,0}(a){(\frac{r}{R_0})}^{-p},$$(A.13)

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 written as $${\Sigma }_d(r)={\displaystyle {\int }_a{\sigma }_d}(r,a)da,$$(A.14)

which we associate with 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 $$x_{a}da=\frac{m(a)dn(a)}{{\displaystyle {\int }_{a_{min}}^{a_{max}}m}(a)dn(a)}=\frac{m(a)dn(a)}{m_{dust}},$$(A.15)

where the size distribution in number dn(a) is assumed to follow a power law of exponent d, $$dn(a)=C{a}^{-d}da,$$(A.16)

where C is a normalization constant. The quantity C is related to the dust mass by simply integrating over the size as follows: $$m_{dust}={\displaystyle {\int }_{a_{\text{min}}}^{a_{\text{max}}}m}(a)dn(a)={\displaystyle {\int }_{a_{\text{min}}}^{a_{\text{max}}}\frac{4\pi }{3}}{\rho }_{m}C{a}^{3-d}da,$$(A.17)

where ρ_m is the specific mass density defined in the previous section. Using Eq.A.15 and the definition of ζ yields $${\sigma }_d(r,a)=x_a\zeta {\Sigma }_g(r)=\frac{\zeta {\Sigma }_g(r)a^{3-d}}{{\displaystyle {\int }_{a_{\text{min}}}^{a_{\text{max}}}{a}^{3-d}}da}.$$(A.18)

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 as follows: $$\frac{m_{d,i}}{m_d}=\frac{{\displaystyle {\int }_{a_{\text{min},i}}^{a_{\text{max},i}}\frac{4}{3}}\pi {\rho }_{\text{grain}}{a}^{3}n(a)\text{d}a}{{\displaystyle {\int }_{a_{\text{min}}}^{a_{\text{max}}}\frac{4}{3}}\pi {\rho }_{\text{grain}}{a}^{3}n(a)\text{d}a}=\frac{a_{\text{max},i}^{4-d}-a_{\text{min},i}^{4-d}}{a_{\text{max}}^{4-d}-a_{\text{min}}^{4-d}}.$$(A.19)

The relativedust mass per bin and the relative area per bin (which also quantifies the number of chemical sites) are given in Table A.1.

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 $${\rho }_d(r,z,a)=\frac{{\sigma }_d(r,a)}{\sqrt{2\pi }{H}_d(r,a)}\exp \left(-\frac{z^2}{2H_d{(r,a)}^2}\right).$$(A.20)

Because dust grains do not react to pressure forces in the same way as gas, the dust scale height differs from that of the gas. Following Boehler et al. (2013), we define the settling factor, that is the ratio of H_d (r, a) to the gas scale height near the midplane H_g(r), as $$s(a,r)=\frac{H_d(a,r)}{H_g(r)},$$(A.21)

which 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 $$T_s(r,z)={\tau }_s{\Omega }_g(r)=\frac{a{\rho }_m}{{\rho }_g(r,z)C_s}\Omega (r).$$(A.22)

The stopping time depends linearly on the grain size. Using Eq. A.9 we can write $$T_s(r,z)=\frac{a{\rho }_m}{{\rho }_g(r,z)H_g(r)}=\frac{\sqrt{2\pi }a{\rho }_m}{{\Sigma }_g(r)}\exp \left(-\frac{z^2}{2H_g{(r)}^2}\right).$$(A.23)

The settlingfactor is given by the approximation of Dong et al. (2015), that is, $$s(a,r)=\frac{H_d(a,r)}{H_g(r)}=\frac{1}{\sqrt{1+T_{s,mid}\frac{S_c}{\alpha }}},$$(A.24)

where α is the α viscosity coefficient, $$T_{s,mid}=\frac{\sqrt{2\pi }a{\rho }_m}{{\Sigma }_g(r)}$$(A.25)

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 of the turbulent viscosity ν_t to the turbulent diffusion D_t as follows: $$S_c=\frac{{\nu }_t}{D_t}.$$(A.26)

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).

Fig. A.2 Settling factor as a function of the dimensionless stopping time at R0. In green the model given by Boehler et al. (2013), the red dots show the simulation values of Fromang & Nelson (2009), and the blue curve indicates the model of Dong et al. (2015) used in this work.

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².

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 shielding is most important only in the upper layers, when dust optical depth remains low while molecular optical depth can be significant.

Fig. B.1 Elevation above the midplane as a function of the radius. The dotted black lines represent the scale heights of the gas (1 H, 2 H, 3 H, and 4 H). The red lines indicate radiation emitted from the central star and from the ISRF. The coordinates (r, z) represent a point at which both fields cross.

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: $$I(\lambda ,r,z)\approx I_{ISRF}{e}^{-{\tau }^V}+I_{*}(r,z)e^{-{\tau }^R},$$(B.1)

where I_ISRF 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 $$\begin{array}{l}{\tau }^V={\tau }_m^V+{\tau }_d^V\hfill \\ {\tau }^R={\tau }_m^R+{\tau }_{d,}^R\hfill \end{array}$$(B.2)

where m stands for molecules and d for dust. So the local expression becomes $$I(\lambda ,r,z)\approx I_{ISRF}{e}^{-({\tau }_d^V+{\tau }_m^V)}+I_{*}{e}^{-({\tau }_d^R+{\tau }_m^R)}.$$(B.3)

The local flux weighted by the dust only I_d can also be approximated by $$I_d(\lambda ,r,z)\approx I_{ISRF}{e}^{-{\tau }_d^V}+I_{*}{e}^{-{\tau }_d^R}.$$(B.4)

The POLARIS code treats cscattering by dust consistently 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: $${\tau }_m(\lambda ,r,z)={\sum }_X{\displaystyle {\int }_{l=f(r,z)}n}(X,r,z)\sigma (X,\lambda )dl,$$(B.5)

where X includes all molecules for which we have cross sections (including H₂O), n(X, r, z) is the number density [cm⁻³] 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 $\sqrt{r^2+z^2}$ between the central star and the considered cell, respectively.

Figure B.2 shows the vertical opacities generated by the dust (dotted lines) and the molecules (solid lines). The opacities are given for different altitudes at 100 au. H₂, CO and N₂ are the main contributors to the line opacity.

As shown in Figs. B.2, at all wavelengths where H₂ contributes to the attenuation (λ ≲ 115 nm), its opacity is inall parts of the disk always substantially larger than the dust opacity.

While calculating the vertical molecular opacity ${\tau }_m^V$ is trivial, characterizing the radial molecular opacity ${\tau }_m^R$ 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 at which molecular shielding occurs cover a limited fraction of the UV domain. In regard to Figure B.1, we introduce the ratio of the radial to vertical molecular optical depths as follows: $$RV=\frac{{\tau }_m^R}{{\tau }_m^V}$$(B.6)

and Eq. B.3 can be rewritten as $$I(\lambda ,r,z)\approx e^{-{\tau }_m^V}(I_{ISRF}{e}^{-{\tau }_d^V}+I_{*}{e}^{-{\tau }_d^R}{e}^{-{\tau }_m^V(RV-1)}).$$(B.7)

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 $$I(\lambda ,r,z)\approx e^{-{\tau }_m^V}(I_{ISRF}{e}^{-{\tau }_d^V}+I_{*}{e}^{-{\tau }_d^R})\approx e^{-{\tau }_m^V}{I}_d(\lambda ,r,z).$$(B.8)

This approximation may underestimate the molecular shielding when ${\tau }_m^V$ is ≲ 1, while ${\tau }_m^R$ is substantially larger than 1, thus impacting the details of the H₂ formation layer and more 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; this situation 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.

Fig. B.2 Vertical optical depth as a function of wavelengths given at different scale heights Hg (0 H, 2 H, 4 H) of a disk at 100 au away from the star. The solid lines represent molecular opacities; the dotted lines indicate dust opacities.

In summary, self and mutual shielding are considered throughout the disk, by applying the attenuation by the vertical opacity due to molecules to the radiation field computed with dust only by the POLARIS code as follows: $$I(\lambda ,r,z)=e^{-{\tau }_m^V}{I}_d(\lambda ,r,z)$$(B.9)

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, pre-computed factors as a function of visual extinction and H₂ column densities (or also CO and N₂ for self-shielding), this approach has the advantage of consistently taking 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₂, 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₂ remains small and limit their contribution.

Fig. B.3 Flux as a function of the wavelength in the midplane (black lines), at 2 H (red lines), and at 4 H (blue lines) at 100 au. The unattenuated UV spectrum is derived from France et al. (2014) with spectral resampling preserving the equivalent widths of the lines.

Appendix C HUV cases

In Figs.C.1 and C.2 we present the vertical profiles at 100 au of the fractional abundances and number densities of gas-phase H, H₂, CO, CS, and CN in HUV single-grain models and multi-grain models, respectively.

C.1 Intermediate models

The intermediate models have physical conditions identical to that of the multi-grain models, but the dust grain distribution is represented using (locally) a single-grain size, determined so that its effective area and mass are identical to those of dust grains at the same location in multi-grain models. The temperature of this equivalent grain is the area-weighted temperature T_a. These models thus exhibit dust settling as multi-grain models. They also have a larger UV penetration than in single-models, similar, although this value is not identical to the UV penetration in multi-grain models.

For simplicity, we only compute the chemistry for the HUV flux at 50, 100, and 200 au. The results are presented in Figures. C.3, C.4, and have been added in Figures 10, 12, and 15 for discussion.

Fig. C.1 Vertical profiles of H, H2, CO, CS, and CN at 100 au from the star of the HUV single-grain models. The dotted line is the abundance relative to H and the solid line is the density [cm−3 ].

Fig. C.2 Vertical profiles of H, H2, CO, CS, and CN at 100 au from the star of the HUV single-grain model HUV-LH-Tg on the left column, and multi-grain models on middle and right columns. The dotted line indicates the abundance relative to H and the solid line indicates the density [cm−3 ].

Fig. C.3 Vertical profiles of H, H2, CO, CS, and CN at 100 au from the star of the single-grain model HUV-LH-Ta (left), the multi-grain model M-HUV-LH (middle), and the intermediate model Set-HUV-LH-Ta (right). The dotted line represents the abundance relative to H and the solid line indicates the density [cm−3 ].

Fig. C.4 Vertical profiles of H, H2, CO, CS, and CN at 100 au from the star of the single-grain model HUV-B14-Ta (left), the multi-grain model M-HUV-B14 (middle), and the intermediate model Set-HUV-B14-Ta (right). The dotted line represents the abundance relative to H and the solid line represents the density [cm−3 ].

Appendix D LUV cases

In Figure D.1 we present the vertical cumulative column density of H₂, CO, CS, and CN at 100 au from the LUV models. The profiles are very similar to that of the HUV models below z ~ 2H (see Fig. 11). We note a few noticeable differences. The model M-LUV-LH produces a larger column density of CS than M-HUV-LH. As for CN, it is clear that M-HUV-B14 produces a larger column density than M-HUV-LH. In the case of LUV multi-grain models, the difference is far less pronounced and M-LUV-LH produces a slightly larger column density than M-LUV-B14. These different results are mostly due to the lower UV penetration in the case of the LUV models than in the HUV models.

In Figure D.2 we present maps of number density for H₂, CO, CS, and CN in the case of LUV single-grain models. Again, the similarity with HUV models is clear. The profiles around the midplane are identicalbetween LUV and HUV. The impact of the different flux on the upper layers is visible. The peak of CO, CS, and CN are located at higher altitudes and are globally wider in LUV models. We note one major difference for CN between LUV and HUV models. HUV-LH-T_a exhibits a clear spike of CN density around 2H while, conversely, LUV-LH-T_a exhibits a drop at the same location. Figs. D.3 and D.4 show the vertical profiles at 100 au of the fractional abundances and number densities of gas-phase H, H₂, CO, CS, and CN in LUV single-grain models and multi-grain models, respectively.

Fig. D.1 Vertical cumulative surface density [cm−2] of H2, CO, CS, and CN at 100 au from the star of the LUV models.

Fig. D.2 Density [cm−3] of H2, CO, CS, and CN in the gas phase of the single-grain models in LUV regime. The left column shows the results of LUV-LH-Tg, the middle column indicates the results of LUV-LH-Ta, and right column indicates the results of LUV-B14-Ta. The black contours represent the dust temperature (Tdust = Tg in the left column and Tdust = Ta in the middle and right columns).

Fig. D.3 Vertical profiles of H, H2, CO, CS, and CN at 100 au from the star of the LUV single-grain models. The dotted line represents the abundance relative to H and the solid line shows the density [cm−3 ].

Fig. D.4 Vertical profiles of H, H2, CO, CS, and CN at 100 au from the star of the LUV single-grain model HUV-LH-Tg in the left column, and multi-grain models in middle and right columns. The dotted line indicates the abundance relative to H and the solid line indicates the density [cm−3 ].

References

All Tables

All Figures

	Fig. 1 Multi-grain model: coupling scheme between NAUTILUS and POLARIS.
In the text

Fig. 2 Physical structure of the multi-grain models incorporating a full dust distribution. The quantity Tg is the gas temperature, nH 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 togas mass ratio, Ta (Eq. (1)) is the area-weighted dust temperature, and na is the area-weighted grain number density. The quantity AV is the visualextinction counted from disk surface toward the disk midplane.

In the text

	Fig. 3 Results of the radiative transfer simulation in vertical direction at a disk radius of r = 100 au. Left: dust density distribution for the 16 grain size intervals as described in Appendices A.2.2 and A.2.3. Middle: dust temperature of the 16 grain size intervals. Right: stellar (solid lines) and ISRF (dotted lines).
In the text

	Fig. 4 Ratio of the absorption cross sections for wavelengths characteristic for stellar radiation (Cabs,1 μm) and dust emission (Cabs,20 μm).
In the text

Fig. 5 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 Ta. 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.

In the text

	Fig. 6 Absorption (blue dashed line) and scattering cross sections (green dotted line) for the first grain size bin (7 nm). The verticalgray lines denote the wavelengths for which the radiation fields are simulated.
In the text

	Fig. 7 H2 formation rate as a function of the uv field intensity (in units of the Habing interstellar field) for various proton number density. Data adapted from Fig. 6 of Bron et al. (2014) are shown as the blue (104 cm−3) and purple (106 cm−3) solid lines, while black lines correspond to our analytic fit for four different number densities (104, 106, 107, 108 cm−3).
In the text

	Fig. 8 Density (cm−3) of H2, CO, CS and CN in the gas phase of the single-grain models in HUV regime. Left column is the HUV-LH-Tg model, the middle column is HUV-LH-Ta, and the rightcolumn is HUV-B14-Ta. The black contours represent the dust temperature (Td = Tg in the left column and Td = Ta in the middle and right columns).
In the text

	Fig. 9 Density (cm−3) of H2, CO, CS, andCN in the gas phase of HUV-LH-Tg (left column) and of the multi-grain models in HUV regime. The black contours represent the dust temperature (Td = Tg in the left column and Td = Ta in the middle and right columns).
In the text

	Fig. 10 Surface densities of the single-grain models (solid lines) and multi-grain models (dotted lines) as a function of the radius in high uv flux regime. The crosses are the intermediate model for a selection of radii.
In the text

	Fig. 11 Vertical cumulative surface density (cm−2) of H2, CO, CS andCN at 100 au from the star in the HUV models.
In the text

	Fig. 12 Species that contain at least 0.5% of elemental carbon (top), nitrogen (middle), and oxygen (bottom) in the midplane at 100 au: percentage of the carbon/nitrogen/oxygen locked in the species as a function of the gas-phase species abundance.
In the text

	Fig. 13 Density (cm−3) of H2, CO, CS, and CN in the gas phase of the single-grain models in the HUV regime. The left column shows the HUV-LH-Tg model, middle column shows HUV-LH-Ta and right shows HUV-B14-Ta. The black contours represent the dust temperature (Td = Tg in the left column and Td = Ta in the middle and right columns).
In the text

	Fig. 14 Density (cm−3) of H2, CO, CS and CN in the gas phase of HUV-LH-Tg (left column) and of the multi-grain models in HUV regime. The black contours represent the dust temperature (Td = Tg in the left column and Td = Ta in the middle and right columns).
In the text

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 106 yr) as a functionof the grain radius a in relation to the grain temperatures. The crosses stand for the multi-grain model at 30 au, the squares stand for the multi-grain model at 100 au, and the circles stand for 200 au. The triangles pointing upward represent the surface abundance in the single-grain model HUV-LH-Tg (grains of size 0.1 μm) and the triangle pointing downward represents the abundance in the single-grain model HUV-LH-Ta, both at 100 au. Abundances being per grain size bin, the sum of the abundances in the multi-grain models should be compared to the abundances in the single-grain case.

In the text

	Fig. 16 Water abundance in the gas phase (top) and on the grain surfaces (bottom) of the single-grain models. The left column shows the abundances in HUV-LH-Tg, middle in HUV-LH-Ta, and right in HUV-B14-Ta. In the top row, the solid black lines show the dust temperature isocontours; in the bottom row they denote isocontours of scale heights.
In the text

	Fig. 17 Water abundance in the gas phase (top) and on the grain surfaces (bottom). The left column shows the abundances in HUV-LH-Tg and the middle and right column in the multi-grain models. In the top row, the solid black lines show the dust temperatureisocontours; in the bottom row they denote isocontours of scale heights.
In the text

Fig. 18 Altitude above the midplane as a function of the abundance relative to H of CH3CN at 100 au. The blue line shows the abundance produced by model M-HUV-LH; the red line shows the abundance produced by model M-HUV-B14. 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 H2 in M-HUV-B14, which enhances the production of CH${}_3^{+}$.

In the text

	Fig. 19 Multi-grain models: maps of ice abundance relative to H in the midplane, with grain sizes as a function of the radius. The top left panel is the distribution in grain temperature. The hottest region is located at around 60 au from the star on grains of size 0.148 μm.
In the text

	Fig. 20 Altitude above the midplane as a function of the abundance relative to H of CH3OH at 100 au. The blue line shows the abundance produced by model M-HUV-LH; the red line shows the abundance produced by model M-HUV-B14. The abundance dramatically drops at ~ 2 scale heights in both models, although M-HUV-B14 maintains a larger abundance than M-HUV-LH at high altitudes.
In the text

	Fig. A.1 Scale heights of the grain distribution for the 16 different grain sizes as a function of the radius. The grain averaged size interval ranges from 0.007 to 651 μm. Small grains scale heights (size < 0.317 μm) practically follow the gas scale height. For larger grains the settling factor (Eq. A.21) becomes significant and the dust scale height decreases.
In the text

	Fig. A.2 Settling factor as a function of the dimensionless stopping time at R0. In green the model given by Boehler et al. (2013), the red dots show the simulation values of Fromang & Nelson (2009), and the blue curve indicates the model of Dong et al. (2015) used in this work.
In the text

	Fig. B.1 Elevation above the midplane as a function of the radius. The dotted black lines represent the scale heights of the gas (1 H, 2 H, 3 H, and 4 H). The red lines indicate radiation emitted from the central star and from the ISRF. The coordinates (r, z) represent a point at which both fields cross.
In the text

	Fig. B.2 Vertical optical depth as a function of wavelengths given at different scale heights Hg (0 H, 2 H, 4 H) of a disk at 100 au away from the star. The solid lines represent molecular opacities; the dotted lines indicate dust opacities.
In the text

	Fig. B.3 Flux as a function of the wavelength in the midplane (black lines), at 2 H (red lines), and at 4 H (blue lines) at 100 au. The unattenuated UV spectrum is derived from France et al. (2014) with spectral resampling preserving the equivalent widths of the lines.
In the text

	Fig. C.1 Vertical profiles of H, H2, CO, CS, and CN at 100 au from the star of the HUV single-grain models. The dotted line is the abundance relative to H and the solid line is the density [cm−3 ].
In the text

	Fig. C.2 Vertical profiles of H, H2, CO, CS, and CN at 100 au from the star of the HUV single-grain model HUV-LH-Tg on the left column, and multi-grain models on middle and right columns. The dotted line indicates the abundance relative to H and the solid line indicates the density [cm−3 ].
In the text

	Fig. C.3 Vertical profiles of H, H2, CO, CS, and CN at 100 au from the star of the single-grain model HUV-LH-Ta (left), the multi-grain model M-HUV-LH (middle), and the intermediate model Set-HUV-LH-Ta (right). The dotted line represents the abundance relative to H and the solid line indicates the density [cm−3 ].
In the text

	Fig. C.4 Vertical profiles of H, H2, CO, CS, and CN at 100 au from the star of the single-grain model HUV-B14-Ta (left), the multi-grain model M-HUV-B14 (middle), and the intermediate model Set-HUV-B14-Ta (right). The dotted line represents the abundance relative to H and the solid line represents the density [cm−3 ].
In the text

	Fig. D.1 Vertical cumulative surface density [cm−2] of H2, CO, CS, and CN at 100 au from the star of the LUV models.
In the text

	Fig. D.2 Density [cm−3] of H2, CO, CS, and CN in the gas phase of the single-grain models in LUV regime. The left column shows the results of LUV-LH-Tg, the middle column indicates the results of LUV-LH-Ta, and right column indicates the results of LUV-B14-Ta. The black contours represent the dust temperature (Tdust = Tg in the left column and Tdust = Ta in the middle and right columns).
In the text

	Fig. D.3 Vertical profiles of H, H2, CO, CS, and CN at 100 au from the star of the LUV single-grain models. The dotted line represents the abundance relative to H and the solid line shows the density [cm−3 ].
In the text

	Fig. D.4 Vertical profiles of H, H2, CO, CS, and CN at 100 au from the star of the LUV single-grain model HUV-LH-Tg in the left column, and multi-grain models in middle and right columns. The dotted line indicates the abundance relative to H and the solid line indicates the density [cm−3 ].
In the text