\section{Introduction}
\label{sec:intro}

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. \revisionA{When} grain growth occurs  (up to at least cm-sized particles)
larger particles, which dynamically decouple from the gas-phase,
settle toward the disk midplane. This vertically changes the
gas-to-dust ratio in the disk. Moreover, grain growth also allows UV to
penetrate deeper in the disk, changing its chemistry.

As a consequence of the complex interplay of physical and chemical
processes, it is now demonstrated that the disk structure can be
separated vertically in three distinct regions: from top to bottom, the
irradiated tenuous disk atmosphere (above $3-4\,H$ or scale heights) where
the gas can be hot and is ionized or in atomic form, the warm denser
molecular layer ($1-3\,H$) where most molecules reside and below $1\,H$, the
midplane. Beyond the CO snowline, which is typically located at a
radius about 20 au in a TTauri disk, the midplane is very dense and
cold ($\leq 20$ $\mathrm{K}$). This area is essentially shielded from the
stellar radiation, with both low ionization and turbulence levels. In
this cold region, most molecules are stuck onto dust grains which
exhibit icy mantles that can still be processed by cosmic rays.

Thanks to ALMA, several high angular resolution and sensitive
observations of disks orbiting either low or intermediate mass young
stars (TTauri and Herbig Ae stars) have quantitatively confirmed this
general scheme. CO observations of HD 163296 have clearly shown the CO
gas-phase depletion in the dark, cold disk midplane \citep{
De-Gregorio-Mionsalvo+etal_2013}, revealing also that the CO gas is
located in the intermediate warm molecular layer.  More recently, the
ALMA observations of the edge-on TTauri disk called the Flying Saucer
\citep{Dutrey+etal_2017} directly show the molecular layer at
intermediate scales ($1-2 H$) with CS, the denser gas tracer being
located at about $1 H$, below the CO emission ($1-3 H$) which better traces
the whole envelope of the molecular layer.

In the last 10 years, several astrochemical models have emerged that
attempt to incorporate the observed complexity of disk physics \revisionA{in order} 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
\citep[e.g.][]{Woitke+etal_2009} and then reactions on the grain
surface \citep[e.g.][]{Semenov+Wiebe_2011, Walsh+etal_2014}.  Another
step was finally achieved with models which include 3 material phases:
the gas-phase and 2 phases for the gas-grain chemistry with grains
modeled by an icy surface layer and an active mantle
\citep{Ruaud+Gorti_2019, Wakelam+etal_2019}. Thermochemical models
which calculate the density and thermal structures in a self-consistent
manner have been also developed \citep[e.g.][]{Hollenbach+Gorti_2009}.
Recently, \citet{Ruaud+Gorti_2019} have also coupled such a
thermochemical model \citep{Gorti+etal_2011} to a 3-phase chemical
model.

An important improvement in chemical models was to introduce more
realistic dust disk structures by taking into account grain growth and
dust settling. Grain growth reduces the dust cross-section and
therefore increases the UV penetration while dust settling locally
changes the gas-to-dust ratio. Grain growth was first introduced in
disk chemistry by  \citet{Aikawa+Nomura_2006}. Dust settling is often
mimicked by adding larger grains in the midplane
\citep[e.g.][]{Wakelam+etal_2019}.

\revisionA{Several authors have already studied the effect of multiple
grain sizes onto chemistry. As an example,  \citet{Acharyya+etal_2011}
used five different grain sizes in their models but the grain
temperature was the same for all grains. Following this work,
\citet{Pauly+Garrod_2016} tested the impact of a simple temperature
distribution, with $T(a) \propto a^{-(1/6)}$, $a$ being the grain size,
while \citet{Ge+etal_2016} explored a small range of temperatures
(14.9--17.9) for the Cold Neutral Medium. More recently, \citet{Iqbal+Wakelam_2018}
also developed a model handling different grain sizes.}

\revisionA{On another hand, \citet{Akimkin+etal_2013} developed
a disk model where they couple the dust grain time evolution (grain
growth and settling) to chemistry. Another improvement was achieved
with the development of grids of thermochemical models allowing
multiwavelength fitting of dust and molecular lines such as DIANA
\citep{Woitke+etal_2019}.
However, these elaborate approaches use a single grain temperature
which is either fixed or derived from a dust radiative transfer
simulation while the temperature of a grain depends on its size
\citep[e.g.][]{Wolf_2003}, larger grains being colder.
\citet{Chapillon+etal_2008} found in the disk of CQ Tau, whose
temperature is well above the CO snowline temperature, an important CO
depletion. They suggested that CO may remain trapped onto larger grains
that are cold enough to prevent thermal CO desorption.
\citet{Heese+etal_2017} have investigated the dust temperature
variations with grain size and position (radius and altitude above the
midplane) in a typical TTauri disk using the 3D dust radiative
transfer code Mol3D \citep{Ober+etal_2015}.  The dust temperature
variations in the molecular layer appear to be significant
enough to affect the disk chemistry. \revisionA{This dependency questions
the applicability of elaborate, self-consistent, thermochemical
models that rely on a single grain size and temperature.} }

\revisionA{We explore here the impact of the variations of the grain temperature
with their size onto the chemistry of a representative TTauri disk. Our goal is not to
build a complete chemical model that incorporates the best of all
previous studies, neither to make grid models that could provide relevant comparisons
with the observations. Instead, we only intend to check the limits of
an assumption that is used in most previous models, and explore the
main consequences of dust grain temperature dependence with size.}

For this purpose, we have coupled the 3D Monte-Carlo continuum
radiative transfer code \textsc{POLARIS}
\citep{Reissl+Wolf+Brauer_2016} to the NAUTILUS Multi Grain Code
\citep{Hersant+etal_2009, Ruaud+etal_2016} using the 3-phase version
which also takes into account a grain size distribution
\citep{Iqbal+Wakelam_2018}. We present the two codes and how we
parametrize the gas and dust disk in Section \ref{sec:mod}. Section
\ref{sec:model_param} deals with the disk model description and global
results. We discuss in more detail the results in Section
\ref{sec:discuss} and we state our conclusions in Section
\ref{sec:concl}.

\section{Model description}
\label{sec:mod}

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

\begin{figure*}
  \centering
  \includegraphics[width=0.7\linewidth]{figures/pipeline3.pdf}
  \caption{Multi-grain model: Coupling scheme between NAUTILUS and POLARIS. }
\label{fig:coupling}
\end{figure*}

\subsection{Disk Model}
\revisionA{Our disk structure is derived from a simple disk model
where truncated power  laws describe most of the quantities: surface density, midplane
temperature, velocity field, with sharp inner and outer radii. A
detailed model description is given in Appendix \ref{app:model}}

\revisionA{We adopt a power-law grain size distribution $dn(a) = C a^{-d} da$ with $dn(a)$, the number of grain
of size $a$ and $C$ a normalization constant that can be derived from the dust mass
(see Appendix \ref{app:dust}). The resulting distribution is
discretized into 16 logarithmically distributed intervals, with grain
sizes ranging from $a_{min} = 5\,$nm to $a_{max}= 1\,$mm. \revisionA{The fraction of mass and area per bin
are given in Table \ref{tab:grains}.}}

\revisionA{Our goal being to evaluate the impact of dust grain temperatures,
we assume a common gas temperature structure for all our simulations (including single and multi-grain models),
in order to avoid chemical effects that could be due to gas only.
We selected values derived from the observations of the Flying Saucer TTauri disk \citep[see][]{Dutrey+etal_2017},
whose edge-on orientation allow a direct measure of the 2D gas
temperature distribution from ALMA observations of the optically thick CO J=2-1 line.
Values are given in Table \ref{tab:param}.}

\revisionA{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 \ref{app:dust}).}
\revisionA{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 $\alpha$ viscosity of 0.01 and a Schmidt number of $S_c = 1$.
Although often assumed in disk modeling, these values are somewhat arbitrary,
\revisionA{and observations with ALMA are not yet accurate enough to allow for a quantitative
description of dust settling.} }

\revisionA{Figure \ref{fig:struct} 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).  }

\begin{table}
\centering
\caption{Overview of the disk model parameters \label{tab:param}}
\begin{tabular}{p{0.68\linewidth}r}
\hline
\noalign{\smallskip}
\multicolumn{2}{c}{\textbf{Disk Parameters}}\\
\noalign{\smallskip}
\hline
 \noalign{\smallskip}
\bm{$T_\mathrm{\star,eff}$} (star temperature)            	& 3900 K \\
\bm{$L_\star$} (star luminosity)       & 0.75 $\mathrm{L_{\odot}}$  \\
\bm{$M_\star$} (star mass)        & 0.58 $\mathrm{M_{\odot}}$ \\
\textbf{ISRF}								  & DRAINE \\

\bm{$\Sigma_{g,0}$} (gas (H + H$_2$) surface density at $R_0$) & 0.335 $\mathrm{g.cm^{-2}}$ \\
\bm{$p$} (surface density exponent) & 1.5 \\
\bm{$R_\mathrm{in}$} (innermost radius)  & 1 $\mathrm{au}$ \\
\bm{$R_\mathrm{out}$} (outermost radius) & 250 $\mathrm{au}$ \\
\bm{$R_0$} (reference radius)& 100 $\mathrm{au}$  \\
\bm{$\zeta$} (dust to gas mass ratio) & 0.01 \\

\bm{$T_\mathrm{mid,0}$} (observed T$_\mathrm{mid}$ at $R_0$)& 10 $\mathrm{K}$\\
\bm{$T_\mathrm{atm,0}$} (observed T$_\mathrm{atm}$ at $R_0$)& 50 $\mathrm{K}$\\
\bm{$\sigma$} (stiffness of the vertical T profile) & 2 \\
\bm{$q$} (radial T profile exponent) & 0.4 \\
\bm{$S_{c}$} (Schmidt number) & 1 \\
\bm{$\alpha$} (viscosity coefficient) & $0.01$ \\

\bm{$d$} (grain size distribution exponent) & $3.5$ \\

\bm{$a_\mathrm{min}$} (smallest grain radius)  & 5 $\mathrm{nm}$ \\
\bm{$a_\mathrm{max}$} (biggest grain radius) & 1 $\mathrm{mm}$ \\
\textbf{number of bins} & 16 \\
\bm{$\rho_\mathrm{grain}$} (material density of dust)  & 2.5 $\mathrm{g.cm^{-3}}$ \\
\bm{$\mathrm{M_{disk}}$} (total disk mass) & 7.5.10$^{-3}$ $\mathrm{M_{\odot}}$\\
\end{tabular}
\end{table}

\begin{figure*}
  \centering

  \includegraphics[width=1\linewidth]{figures/physical4.pdf}
\caption{Physical structure of the multi-grain models incorporating a
full dust distribution. $T_g$ is the gas temperature, $n_H$ is the gas
number density, $G$ is the local stellar and interstellar field in the
case of the HUV models (see Appendix \ref{app:model}), $\zeta$ is the
dust to gas mass ratio, \revisionA{$T_a$ (Eq.\ref{eq:ta}) is the area-weighted
dust temperature, and $n_a$ is the area-weighted grain number density.
A$_V$ is the visual extinction counted from disk surface toward the disk midplane.}
\label{fig:struct}}
\end{figure*}

\subsection{Radiation fields and dust temperature
\label{sec:rad}
In this new approach, we use the dust radiative transfer code to
calculate  throughout the disk the UV field  and the dust temperature
for each of the 16 grain sizes, in each cell of the disk (in radius and
altitude). The dust temperatures are then used to compute
the  chemistry with NAUTILUS.

\revisionA{However, the UV field derived above only accounts for dust extinction and
scattering. A proper account for self and mutual shielding of molecules
and atoms is required to handle the photo-chemical processes.}

\subsubsection{Radiation sources}
\label{sec:sources}
As the central heating source we assume a low mass pre-main sequence
star with a mass of $M_\star = 0.58$ $\mathrm{M}_\odot$, which radiates as a black
body with an effective temperature of $T_\mathrm{{\star, eff}}$ = 3900
$\mathrm{K}$ and a luminosity $L_\star$ = 0.75 $\mathrm{L_\odot}$.
Additionally, we consider the interstellar radiation field (ISRF) with
a spectral energy distribution (SED) from \citet{Draine_1978} between
91.2 $\mathrm{nm}$ and 200 $\mathrm{nm}$ and the extension of
\citet{Dishoeck+Black_1982} for longer wavelengths. Besides the stellar
radiation field, TTauri stars exhibit a significant UV excess coming
from the inner disk boundary and accretion shocks on the stellar
surface. This contribution is highly wavelength-dependent, with
spectral lines like Lyman $\mathrm{\alpha}$ having substantial
intensities \citep{France+etal_2014}. As concerns the thermal balance of dust, the frequency
dependence is largely irrelevant, and can be absorbed by ensuring that
the luminosity used in the model contains the UV excess contribution.
However, photoprocesses being wavelength dependent, the details of the UV spectrum shape
strongly affect the chemistry. For this reason, we take as stellar
input spectrum for our radiative transfer and chemical model the
sum of the stellar black body spectrum and the UV excess typically
found in the in TTauri stars. Further details about the
radiation fields are given in Section \ref{sec:set}.

\subsubsection{Radiative transfer simulations}\label{sec:RT_sim}

The radiation field is computed using the 3D Monte-Carlo continuum and
line radiative transfer code \textsc{POLARIS}
\citep{Reissl+Wolf+Brauer_2016}.
As in \citet{Heese+etal_2017},
we calculate  the dust temperature distribution of 16 grain size
intervals (see Fig.\,\ref{fig:RTsimulation}).
\revisionA{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 \citet{Draine_1984} and
\citet{Laor+Draine_1993} are used to calculate the scattering
and absorption coefficients based on Mie scattering \citep{Mie_1908}
using the \textsf{MIEX} code \citep{Wolf+Voshchinnikov_2004}.
\revisionA{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 \revisionA{1.05}
between 1 and 300 au, and 181 angles with the same factor between 0 and \revisionA{$\pi$}.
}

The temperature in the upper optically thin disk layers strongly
depends on the grain size (Fig.\,\ref{fig:RTsimulation}, middle panel).
The dust temperature increases with grain size for sizes in the range
0.007 $\mu$m to 0.070 $\mu$m, then decreases for sizes in range 0.32
$\mu$m to 15 $\mu$m, and increases again, although only slightly, with
larger grain sizes. As explained in \citet{Heese+etal_2017},
this can be understood by looking at the ratios of
the absorption cross-section for wavelengths at which the star radiates
($\approx 1$ $\mu$m) and the absorption cross-section for wavelengths at
which the dust emits ($\approx 20$ $\mu$m, Fig.\,\ref{fig:CabsRatio}). The
ratio increases up to grain sizes of 0.07 $\mu$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 $\mu$m. Using these radiative
transfer simulations we can obtain a distribution of dust temperatures
$T_d(a_i, r, z)$ that fully depend on the size and location of
the grains.
To compare our set of models with those consisting of a
single grain size (see Section \ref{sec:single}), it is also convenient
to introduce an area-weighted temperature defined as:

\begin{equation}
\label{eq:ta
	T_{a}(r, z) =  \frac{ \sum_i a_i^2 T_d(a_i, r, z) n_d(a_i, r, z) }{ \sum_i a_i^2 n_d(a_i, r, z) }.
\end{equation}

\noindent We note that the area-weighted dust density $n_\mathrm{a}(r,
z)$ (Fig.\,\ref{fig:struct}) is defined by an equation of the same
form. Figure \ref{fig:weighted_100} shows the dust temperature profiles
at 100 $\mathrm{au}$ used by all our different models.

Our disk model is less massive than that used by \citet{Heese+etal_2017},
leading to smaller opacities, and thus the temperature differences between
grains of different sizes is much larger in Figs.\ref{fig:RTsimulation} and \ref{fig:weighted_100}
than in that previous work.

\begin{figure}
\centering

\includegraphics[width=1.05\linewidth]{figures/C_90-400um.pdf}
\caption{Absorption (\textit{blue dashed line}) and scattering
cross-sections (\textit{green dotted line}) for the first grain size
bin (7 $\mathrm{nm}$). The vertical gray lines mark the wavelengths for
which the radiation fields are simulated.
\label{fig:opticalProperties}}
\end{figure}

\begin{figure}
\centering

\includegraphics[width=0.95\linewidth]{figures/Cabs.pdf}
\caption{Ratio of the absorption cross-sections for wavelengths characteristic for stellar radiation ($\mathrm{C_{abs, 1 \mu m}}$) and dust emission ($\mathrm{C_{abs, 20 \mu m}}$).     \label{fig:CabsRatio}}
\end{figure}

\begin{figure*}
\includegraphics[width=1.0\textwidth]{figures/Polaris_neu2.pdf}
\caption{Results of the radiative transfer simulation in vertical direction at a disk radius of $r$ = 100 $\mathrm{au}$. \textit{Left:} Dust density distribution for the 16 grain size intervals as described in Appendix \ref{app:distrib} and \ref{app:settle}. \textit{Middle:} Dust temperature of the 16 grain size intervals. \textit{Right:} Stellar (solid lines) and interstellar radiation field (dotted lines).
\label{fig:RTsimulation}}
\end{figure*}

\begin{figure}
\includegraphics[width=\linewidth]{figures/weighted_temperature_100.pdf
\caption{\label{fig:weighted_100} Vertical profiles of temperatures at 100 $\mathrm{au}$. The thick black line is the gas vertical temperature profile in Kelvin. The dotted colored lines are temperatures of each grain population and the thick red line is the area-weighted temperature $T_{a}$. We see that temperatures are roughly constant for $z/H$ < 1.5 and for $z/H$ > 2.5. Between these two altitudes the temperatures exhibit a strong transition.}
\end{figure}

\subsubsection{Dust extinction, self and mutual shielding}
\label{sec:UVfield}
The radiation field is a contribution of both the interstellar
radiation field (ISRF) and the stellar radiation field (SRF). Radiation
from these two sources encounters molecular or atomic gas and dust
grains that affects its propagation through extinction and scattering.
The effect of the dust is computed by the 3D Monte-Carlo continuum and
line radiative transfer code POLARIS as described in
Sec.\,\ref{sec:RT_sim} \revisionA{(Fig.\ref{fig:opticalProperties} shows the sampling
in the 90 to 400 $\mathrm{nm}$ wavelengths range relevant for photoprocesses)}.
Figure \ref{fig:RTsimulation} shows results at a
disk radius of $r =100$ $\mathrm{au}$.

In chemical codes, the effect of extinction due to gas, known as self
and mutual shielding, is often estimated using empirical attenuation
factors as a function of line-of-sight visual extinction
\cite[e.g.][for a discussion]{Lee+etal_1996}. \revisionA{The situation for
disks is more complex, due to the two  sources of UV radiation (stellar
neighborhood and ISRF) and varying dust properties. Empirical
attenuation factors have been derived for disks \citep[see for
example][]{Visser+etal_2009,Heays+etal_2017}. They are in general used
in a 1+1D  approach, where the two sources of radiation are treated
independently and scattering is neglected. Furthermore, they  depend
implicitely on the dust settling that was assumed during their
derivation, as well as on the shape of the incident UV field.}

\revisionA{For better consistency}, we use here a procedure that
builds on our knowledge of the dust-attenuated UV field. The extinction
produced by the gas requires to know the opacity generated by the
molecules through which the radiation travels from outside the disk to the
local cell.

\revisionA{We show in Appendix \ref{app:self} 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:}
\begin{equation}
\label{eq:uvfield
	I_L(\lambda, r, z) =  e^{-\tau^V_m} I_d(\lambda, r, z)
\end{equation}
where $I_d(\lambda, r, z)$ is the radiation field given by POLARIS, which
explicitly includes the impact of dust scattering, \revisionA{and $\tau^V_m$ is the
total opacity due to molecules from the $(r,z)$ point to $(r,+\infty$).}

\subsection{Chemistry simulations} \label{sec:chemistry

The NAUTILUS gas-grain code \citep{Ruaud+etal_2016} is used to perform
chemistry simulations. NAUTILUS is dedicated to computing chemical
evolution in different interstellar environments and uses the rate
equation approach \citep{Hasegawa+etal_1992, Hasegawa+Herbst_1993b} to calculate
the abundance as a function of time. The NAUTILUS gas-grain code is a
3-phase model \citep{Ruaud+etal_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, the diffusion of these species and their
reactions, desorption of species from the surface and repopulation of
the surface by species from the mantle as the species on the surface
evaporate. The latest version of NAUTILUS, the Nautilus Multi Grain Code (NMGC)\footnote{Code available here: \url{https://github.com/sgvn/NMGC}} 0D model
\citep{Iqbal+Wakelam_2018}, can perform simulations using a full
grain distribution in size.
\revisionA{Each grain size bin is treated independently of the others, and only interact with gas.}
\revisionA{Compared to the single-grain model (model using a single grain size), a multi-grain model changes
chemical rates in two ways. First, accretion rates of species depend on
grain population. In a multi-grain model, accretion rates vary on
different grain sizes according to their total surface area (Table
\ref{tab:grains} shows the relative surface areas as function of grain sizes).
Second, dust temperature depends on its
size (see Fig.\,\ref{fig:RTsimulation}). Species have higher reaction
rates on grains which are hotter due to higher hopping and desorption
rates. In general, grains which are hotter have lower abundances of
lighter species such as CN, CH$_2$, CO, etc. and more of heavier
species such as H$_2$O, CH$_3$OH, etc\ldots \citep[see][for details]{Iqbal+Wakelam_2018}.
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 \revisionA{cylindrical} coordinates $(r,z)$ where
we provide as input the gas temperature (Eq. \ref{eq:verticalT}),
the local radiation field $I_d(\lambda,r,z)$ and dust temperatures
$T_d(a_i, r, z)$ given by the radiative transfer simulations,
and the number density of each dust population $n_d(a_i, r, z)$ in order to compute the local
grain abundances relative to the number density of Hydrogen nuclei
$n_H(r, z)$:
\begin{equation}
\label{eq:rapport}
	AB_d(a_i, r, z) = \frac{n_d(a_i, r, z)}{n_H(r, z)}.
\end{equation}
\revisionA{The $z$ direction is treated as in \citet{Hersant+etal_2009}, but using 64 points.}

\subsection{Photorates}

Radiations in the UV range have a critical impact on the disk chemistry
since UV photons have the necessary energy to photoionize or
photodissociate molecules in the disk. For this reason, characterizing
the photorates [s$^{-1}$] at which molecules are dissociated or ionized
is of high importance in chemistry models and necessary for a full
interpretation of observations.

We evaluate the photoprocess rate as in \citet{Heays+etal_2017}:
\begin{equation}
\label{eq:rate
	k_{p}(X, r, z) =  \int_\lambda  \sigma_{p}(X, \lambda) I_L(\lambda, r,z) d\lambda
\end{equation}
where $\sigma_{p}(X, \lambda)$ is the photoprocess cross-section of
species $X$ at wavelength $\lambda$ and $I_L(\lambda,r,z)$ is the local
UV radiation. The index $p$ equals $i$ when it stands for ionization
and equal $d$ when it stands for dissociation. We perform our
integration over the 91.2 to 400 $\mathrm{nm}$ wavelength range. To
evaluate the rate at each point in the disk, we thus need the local UV
radiation field and the (space-independent) cross-sections. Following
Section \ref{sec:UVfield}, the UV field is derived using Eq.
\ref{eq:uvfield}.

The molecular cross-sections used in this study are
extracted from the Leiden University
website\footnote{\url{https://home.strw.leidenuniv.nl/~ewine/photo/index.php?file=cross_sections.php}}
\citep{Heays+etal_2017}. They are collected from experimental and
theoretical studies. \revisionA{Eighty molecules are treated in this
way.}
\revisionA{The photorates are recomputed at each explicit time step of the chemical evolution.}

\revisionA{For all other gas constituents \citep[atoms or molecules unavailable in][]{Heays+etal_2017}, we
use approximate shielding factors from \cite{Visser+etal_2009}}.

\subsection{Molecular Hydrogen formation at the disk surface}
\label{sec:h2_formation

\subsubsection{The atomic H problem}
The NAUTILUS code initially implemented the $\mathrm{H_2}$\xspace formation through
the Langmuir-Hinshelwood (LH) mechanism that considers physisorption on
grain surfaces. The LH mechanism is efficient only over a relatively
narrow temperature range
\citep[5-15\,K on flat surfaces][]{Katz+etal_1999,Vidali+etal_2004,Vidali+etal_2005}.

\revisionC{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 $\mathrm{H_2}$\xspace formation via the LH mechanism as treated in most astrochemistry codes} becomes inefficient, leading
to much smaller $\mathrm{H_2}$\xspace formation rate than in the case of a single,
larger, equivalent-area, grain size with lower surface temperature. The low formation rate of $\mathrm{H_2}$\xspace leaves
a significant amount of atomic Hydrogen that can severely affect the
chemical balance in the disk upper layer (see Section \ref{sec:discuss} for a
more detailed discussion).

Observations of $\mathrm{H_2}$\xspace~in unshielded regions \citep{Habart+etal_2004,
Habart+etal_2011} where dust can reach high
temperatures ($\geq20$ 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 $\mathrm{H_2}$\xspace can be efficient in a wider domain than predicted by the LH
mechanism at equilibrium temperature.

\revisionC{To evaluate the $\mathrm{H_2}$\xspace formation rate, many studies adopt a canonical value for interstellar clouds derived by \citet{Jura_1975}, which is defined by taking half the rate at which hydrogen atoms stick to the grain surface  \citep[e.g.][]{Walsh+etal_2014, Agundez+etal_2018}}.

Several \revisionC{other} studies proposed more sophisticated mechanisms that theoretically produce comparable
$\mathrm{H_2}$\xspace column densities to those derived by observations
\citep[e.g.][]{Duley_1996,Cazaux+Tielens_2004,Cuppen+Herbst_2005,Iqbal+etal_2012,
Iqbal+etal_2014,Thi+etal_2020}. We follow here the approach
of \citet{Bron+etal_2014}, hereafter B14, which \revisionC{considers the stochastic fluctuations of both the grain temperature and the H-atom surface population.}

\subsubsection{Effects of temperature fluctuations}

Small grains ($\leq$ 0.01 $\mu$m), given their
very small heat capacity, have their temperature fluctuating widely
even when a relatively small quantity of energy is absorbed. In typical
unshielded regions, a small grain can absorb UV photons originated from
the interstellar field with a rather high probability, causing a sudden
spike in temperature. If the interval of time between UV photon
absorptions is sufficiently long, the grain which has undergone
transient heating has the time to cool down to temperatures smaller than
the equilibrium temperature, giving the opportunity for $\mathrm{H_2}$\xspace to form on
the surface. Then, the rapid heating of the grain by the absorption of
a UV photon will thermally desorb most species on the surface,
releasing $\mathrm{H_2}$\xspace molecules to the gas phase.

\citet{Bron+etal_2014} have studied the effect of these temperature fluctuations on the formation of $\mathrm{H_2}$\xspace in PDRs and
high UV radiation regions using a statistical approach. By numerically
solving the master equation \citep{LeBourlot+etal_2012}, they show
\citep[][see their Fig.6]{Bron+etal_2014} that the LH mechanism alone is
sufficient to reach the typical $\mathrm{H_2}$\xspace formation rates observed  in PDRs
when the temperature fluctuations of small grains are accounted for.

\subsubsection{Implementing the $\mathrm{H_2}$\xspace formation rates}
The investigation of B14 lies in a
domain of gas temperature (100 K) and gas density (10$^0$ to 10$^6$
$\mathrm{cm^{-3}}$) distinct from our study
but this difference can easily be
corrected as both temperature and gas density only affect the results
through the collision rate of atomic Hydrogen with the grains, which
scales as $n_H\sqrt{T_{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's data the value corresponding to an equivalent
gas density $n_{eq}$ equal to:
\begin{equation}
\label{eq:neq
	n_{eq} =  n_H\sqrt{\frac{T_g}{100\,\mathrm{K}}}\frac{s(T_g)}{s(100\,\mathrm{K})}
\end{equation}
where s($T_g$) is the sticking coefficient at temperature $T_g$ as defined by \citet{LeBourlot+etal_2012}:

\begin{equation}
\label{eq:sticking
	s(T) = \left(1+\frac{T}{T_2}^{\beta}\right)^{-1}
\end{equation}
with $T_2= 464$\,K and $\beta=1.5$, which approximately gives the same result as in \citet{Sternberg+Dalgarno_1995}.

We interpolated B14 data as a function of proton density and UV field
intensity using an analytic polylogarithmic function. The interpolation
results are displayed in Fig.\ref{fig:b14}. The PDR layers of
protoplanetary disks are dense, typically above
$10^6 -- 10^8$\,cm$^{-3}$ so that essentially the upper curves in
Fig.\ref{fig:b14} apply to our case.

\begin{figure}
  \centering

  \includegraphics[width=\linewidth]{figures/fig6_fitc.pdf}
\caption{\label{fig:b14}
$\mathrm{H_2}$\xspace~ 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
\citet{Bron+etal_2014} are the blue (10$^{4}$ cm$^{-3}$) and
purple (10$^{6}$ cm$^{-3}$) solid lines while black lines correspond
to our analytic fit
for four different number densities (10$^{4}$, 10$^{6}$, 10$^{7}$, 10$^{8}$ cm$^{-3}$).}
\end{figure}

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 of the fitted
coefficients.
However, the asymptotic behaviors of the rates
should remain similar. The rate is proportional to the density at low
densities, and progressively saturates for larger densities as the
fraction of atomic H decreases because the UV flux is insufficient
to dissociate H$_2$. As a function of the UV field intensity the (unsaturated) rate is expected to decrease as $1/G$ since the
time spent by any grain in the relevant temperature range where the
LH formation is efficient is inversely proportional to the UV photon flux.

\section{Model parameters and results}
\label{sec:model_param}

\subsection{Detailed description of models}\label{sec:set

We implement a set of \revisionC{twelve} different models (see Table \ref{tab:full3}).
All models \revisionA{share a common set of disk structure and physical parameters (see
Fig.\,\ref{fig:struct}) as described in Table \ref{tab:param}, but
differ for a limited number of assumptions on grain sizes and
temperatures,  UV field or H$_2$ formation mechanism.}

The set is divided into \revisionC{three groups, i.e. the single-grain models,
the multi-grain models and an intermediate category of models with single grains locally varying
to mimic the grain distribution used in the multi-grain model.}
 We describe their specific characteristics in
detail in Section \ref{sec:single}, \ref{sec:multi} and  \ref{sec:interm}. As mentioned in
Section \ref{sec:sources}, the excess of UV radiation exhibited by
typical TTauri stars, given the strong wavelength-dependence of
photoprocesses, inevitably impacts the chemistry of the disk. In order
to have a broad view on how our chemical models react to different
intensities of UV flux, we adopt two regimes of radiation. One regime
uses the spectrum of a relatively high-UV emitting TTauri star and
another uses the spectrum of a relatively low UV emitting one,
hereafter called HUV and LUV, respectively.
For the HUV regime we use
the spectrum of TW Hya \citep{Heays+etal_2017} while for the LUV we use
the one of DM Tau. Both FUV spectra are provided in the CTTS FUV
spectra
database\footnote{\url{https://cos.colorado.edu/~kevinf/ctts_fuvfield.html}}
\citep{France+etal_2014}.
The blue solid line in Fig.\,\ref{fig:flux_approx} exhibits the spectrum
of TW Hya as provided by \citet{France+etal_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 \ref{tab:init_ab} and all the chemical
output are extracted after an \revisionC{evolution time of $5.10^6$ yrs. We note that using molecular initial abundances would give almost the same final results for disks of this age  \citep[see][]{Wakelam+etal_2019}}. All
models incorporate the same chemical network and chemical processes
described in Section \ref{sec:chemistry}. \revisionC{ We ignore X-rays because of the large cumulative column densities found around the midplane, the main interest 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$\times 10^{-17}$ s$^{-1}$ in the entire disk \citep[see][for a detailed discussion on the ionization rate]{Cleeves+etal_2015}}.

\subsubsection{Single-grain models}\label{sec:single}

We use single-grain models as a comparison to our more sophisticated
multi-grain models described in Section \ref{sec:multi}.
\revisionA{These models are treated as in \citet{Wakelam+etal_2016}, in
particular with respect to the UV field penetration and photoprocesses \citep[see][]{Hersant+etal_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 $\mathrm{\mu m}$.
The grains are settled toward the midplane assuming a dust scale
height as defined in \autoref{eq:rhod}. However, 0.1 $\mathrm{\mu m}$
grains exhibit little settling, as it can be seen from Fig.\ref{fig:hd}
and  such models are virtually equivalent to non-settled dust models
\citep{Wakelam+etal_2016}.

\paragraph{Model 1: xUV-LH-$\mathrm{T_{g}}$}
The dust temperature ($T_d$) profile follows that of the gas ($T_g$)
which is given by the solid black line in Fig. \ref{fig:weighted_100}),
$T_d$ is equal to $T_g$ everywhere in the disk. $\mathrm{x}$ means either $\mathrm{H}$ for high or $\mathrm{L}$ for low, hence Model 1 is either
HUV-LH-$\mathrm{T_{g}}$\xspace or LUV-LH-$\mathrm{T_{g}}$\xspace, depending on whether the stellar UV field is high (HUV)
or low (LUV). The differences are detailed in Section \ref{sec:set}.

\paragraph{Model 2: xUV-LH-$\mathrm{T_{a}}$}
Model 2 is similar to Model 1, except for the dust temperature. In
Model 2, the grain temperature is given by Eq. \ref{eq:ta} and shown in
Fig.\,\ref{fig:weighted_100} (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.

\paragraph{Model 3: xUV-B14-$\mathrm{T_{a}}$}
In Model 3, the prescription of \citet{Bron+etal_2014} as described in
Section \ref{sec:h2_formation} is used to form $\mathrm{H_2}$\xspace. The expected effect
is a higher production of $\mathrm{H_2}$\xspace in the PDRs region of the disk than in
the other models.

\subsubsection{Multi-grain models}\label{sec:multi}
Multi-grain models, as opposed to single-grain ones, include the full
grain distribution detailed in Table \ref{tab:grains} where each grain
population has a specific temperature (Fig. \ref{fig:weighted_100}), a
specific settling factor and is chemically active. For all multi-grain
models, we use a size range of [5 $\mathrm{mm}$, 1 $\mathrm{mm}$] with
a size exponent $d=-3.5$ (Eq. \ref{eq:da}) and divide the grain size
interval into 16 logarithmically distributed subintervals using Eq.
\ref{eq:interval}. An overview of the grain size interval with their
respective relative dust mass and surface area is given in
\autoref{tab:grains} while Fig.\,\ref{fig:hd} shows the scale height of
each 16 grain size population as a function of the radius. Multi-grain
models are composed of two categories.

\paragraph{Model 4: M-xUV-LH}
Either called M-HUV-LH\xspace or M-LUV-LH\xspace depending on the flux, they use the classical LH mechanisms for $\mathrm{H_2}$\xspace formation without considering small grains temperature fluctuation.

\paragraph{Model 5: M-xUV-B14}
We add the $\mathrm{H_2}$\xspace formation mechanisms prescribed in B14 (see Section \ref{sec:h2_formation}).
These models, with grain growth, dust settling, one temperature per grain size and H$_2$ formation derived from \citet{Bron+etal_2014} are the most complex models we compute.

\subsubsection{Intermediate models}\label{sec:interm}

 These models use locally a single grain size mimicking the grain growth and dust settling of the multi-grain models (see Appendix.\,\ref{app:interm} for details).  At each point, the single grain size, which is chosen in order to obtain
the same grain mass and area, are summed up over the grain size distribution of the multi-grain models.
The purpose of these models is to better illustrate the impact of a dust temperature spread onto the
chemistry, therefore we only run these models at radii 50, 100 and 200 au.

 \paragraph{Model 6: Set-HUV-y-$\mathrm{T_{a}}$}
There are two models. The dust temperature is T$_a$ while the H$_2$ formation (y) either follows the LH mechanism or the B14 approximation.
\\\\
\revisionA{Table \ref{tab:full3} summarizes the various models.}
\revisionA{Table \ref{tab:fig} lists the different figures in the text in relation with their respective models.
A few figures relevant of the HUV cases are presented in Appendix \ref{app:huv},
while all figures for the LUV cases appear in Appendix \ref{app:luv}.}

\begin{table}
\caption{List of figures. \label{tab:fig}}
\centering
\begin{tabular}{ l c }

   \hline
   \noalign{\smallskip}
  \textbf{Models} &  \textbf{Figures} \\
  \noalign{\smallskip}
 \hline
 \noalign{\smallskip}
   \textbf{single, HUV} & \ref{fig:compare_high}, \ref{fig:100cumul_huv},  \ref{fig:reservoirs}, \ref{fig:s-s-maps-huv},  \ref{fig:ab_surface},  \ref{fig:s-water_maps}, \ref{fig:s-100profile_high}    \\
   \noalign{\smallskip}
   \hline
   \noalign{\smallskip}
    \textbf{single, LUV} & \ref{fig:100cumul}, \ref{fig:s-maps-luv}, \ref{fig:s-100profile_low}  \\
    \noalign{\smallskip}
   \hline
   \noalign{\smallskip}
   \hline
   \noalign{\smallskip}
   \textbf{multi, HUV} &  \ref{fig:compare_high}, \ref{fig:100cumul_huv}, \ref{fig:m-maps},   \ref{fig:reservoirs}, \ref{fig:m-s-maps-huv}, \ref{fig:ab_surface}, \ref{fig:m-water_maps}, \ref{fig:methyl_compare}, \ref{fig:CH3OH_compare}, \ref{fig:surfdensmap}, \ref{fig:m-100profile_high}    \\
   \noalign{\smallskip}
      \hline
      \noalign{\smallskip}
   \textbf{multi, LUV} & \ref{fig:100cumul}, \ref{fig:m-100profile_low}    \\
      \hline
      \noalign{\smallskip}
   \textbf{Set, HUV} &\ref{fig:compare_high}, \ref{fig:reservoirs}, \ref{fig:ab_surface} , \ref{fig:interm-100profile_LH}, \ref{fig:interm-100profile_B14}   \\
\end{tabular}
\end{table}

\subsection{Results}

\revisionA{We present here some general issues
about the UV field, the dust temperature and the H$_2$ formation. A more thorough discussion of}
the chemistry is done in the next section.

\begin{table*}
\centering
 \caption{Description of models}
 \begin{tabular}{l  c  c  c  c  c  c}
   \hline
    & \textbf{bins} & \textbf{sizes} (\bm{$\mathrm{\mu m}$}) & \bm{$\mathrm{T_{d}}$} & \textbf{HUV} & \textbf{LUV} & \textbf{H$_2$ formation} \bigstrut \\
   \hline

   \textbf{HUV-LH-$\bm{\mathrm{T_{g}}}$} & 1 & 0.1 & $T_{g}$ & \checkmark &  & LH \bigstrut \\
   \hline
   \textbf{HUV-LH-$\bm{\mathrm{T_{a}}}$}  & 1 & 0.1 & $T_{a}$ & \checkmark &  & LH  \bigstrut \\
   \hline
    \textbf{HUV-B14-$\bm{\mathrm{T_{a}}}$} & 1 & 0.1 & $T_{a}$ & \checkmark &  & B14  \bigstrut \\
   \hline

    \textbf{LUV-LH-$\bm{\mathrm{T_{g}}}$} & 1 & 0.1 & $T_{g}$ &  & \checkmark & LH  \bigstrut \\
   \hline
    \textbf{LUV-LH-$\bm{\mathrm{T_{a}}}$} & 1 & 0.1 & $T_{a}$ &  & \checkmark & LH  \bigstrut  \\
   \hline
     \textbf{LUV-B14-$\bm{\mathrm{T_{a}}}$} & 1 & 0.1 & $T_{a}$ &  & \checkmark & B14 \bigstrut \\
   \hline
   \hline
   \textbf{M-HUV-LH} & 16 & [$5.10^{-3} - 10^{3}$] & $T_{i}$ $ _{\{i = 1 ... 16\}}$ & \checkmark &  & LH  \bigstrut \\
      \hline
   \textbf{M-HUV-B14} & 16 & [$5.10^{-3} - 10^{3}$] & $T_{i}$ $ _{\{i = 1 ... 16\}}$ & \checkmark &  & B14  \bigstrut \\
      \hline
   \textbf{M-LUV-LH} & 16 & [$5.10^{-3} - 10^{3}$] & $T_{i}$ $ _{\{i = 1 ... 16\}}$ &  & \checkmark & LH  \bigstrut \\
      \hline
   \textbf{M-LUV-B14} & 16 & [$5.10^{-3} - 10^{3}$] & $T_{i}$ $ _{\{i = 1 ... 16\}}$ &  & \checkmark & B14  \bigstrut \\
      \hline
     \textbf{Set-HUV-LH-$\bm{\mathrm{T_{a}}}$} & 1 & See appendix C  &  $T_{a}$ &\checkmark & & LH  \bigstrut \\
      \hline
     \textbf{Set-HUV-B14-$\bm{\mathrm{T_{a}}}$} & 1 & See Appendix C & $T_{a}$  & \checkmark &  & B14  \bigstrut \\
     \hline
    \end{tabular}
 \tablefoot{HUV stands for high UV
 field and LUV for low UV field.
 LH is for the classical traitement of the Langmuir-Hinshelwood
 mechanism and B14 when the prescription of \citet{Bron+etal_2014} is
 used instead.  $T_g$ means that the dust temperature equals to that of
 the gas, $T_a$ is the weighted-area dust temperature as defined by
 Eq.\,\ref{eq:ta}. $T_i$ is the temperature of i$^\mathrm{th}$
 grain population. M stands for multigrain models and Set for dust models mimicking grain growth and dust settling similar to
 multi-grain models. As an example, the
 model M-HUV-B14\xspace is a multi-grain model computed using the B14 method with a
 high UV field.  \label{tab:full3}}

\end{table*}
\begin{figure*}
\begin{subfigure}{.48\linewidth}
  \centering

  \includegraphics[width=1\linewidth]{figures/compare/surfacedensities/Hb.pdf}
\end{subfigure}
\begin{subfigure}{.48\linewidth}
  \centering

  \includegraphics[width=1\linewidth]{figures/compare/surfacedensities/H2b.pdf}
\end{subfigure}

\begin{subfigure}{.48\linewidth}
  \centering

  \includegraphics[width=1\linewidth]{figures/compare/surfacedensities/COb.pdf}
\end{subfigure}
\begin{subfigure}{.48\linewidth}
  \centering

  \includegraphics[width=1\linewidth]{figures/compare/surfacedensities/CSb.pdf}
\end{subfigure}

\begin{subfigure}{.48\linewidth}
  \centering

  \includegraphics[width=1\linewidth]{figures/compare/surfacedensities/CNb.pdf}
\end{subfigure}
\begin{subfigure}{.48\linewidth}
  \centering

  \includegraphics[width=1\linewidth]{figures/compare/surfacedensities/H2COb.pdf}
\end{subfigure}
\caption{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. \revisionC{Crosses are the intermediate model for a selection of radii}.}
\label{fig:compare_high}
\end{figure*}

 \subsubsection{Impact of UV field}

\revisionA{A comparison between Tables
\ref{tab:huv_sigma} and \ref{tab:luv_sigma} shows that, for a given
model, the integrated column densities (or molecular surface densities)
have about the same order of magnitude whether the flux is high (HUV)
or low (LUV). Unlike in the PDR layer, the midplane chemistry is not
dependent on the UV flux regime. The dust opacity is such that the UV
flux is sufficiently attenuated in the HUV and LUV models for the
chemistry to be similar in both cases.}

\revisionA{Figures \ref{fig:s-100profile_high},
\ref{fig:m-100profile_high}, \ref{fig:s-100profile_low}, and
\ref{fig:m-100profile_low}  show that the chemistry is nearly the same
for all models from the midplane to roughly $1.5\,H$.  We observe more
$\mathrm{H_2}$\xspace, CO, CS and CN above 3 scale heights in the LUV models than in the
HUV ones because a high UV flux implies more efficient
photodissociation. Near the midplane, the UV field is similarly
attenuated and the \revisionA{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.}

\subsubsection{Impact of  dust temperature on H$_2$ formation and chemistry}

\revisionA{A first effect of the adopted grain temperature is related to
the formation of $\mathrm{H_2}$\xspace and the remaining amount of H, which is very
important for the gas-phase chemistry. Adopting the formalism of
\citet{Bron+etal_2014} to form $\mathrm{H_2}$\xspace produces, as expected, more $\mathrm{H_2}$\xspace. In
these cases, Figures \ref{fig:s-100profile_high}  and
\ref{fig:m-100profile_high} show a shift in the H/$\mathrm{H_2}$\xspace transition toward
the disk surface which implies less gas-phase H and more $\mathrm{H_2}$\xspace in the
upper layers.}

\revisionA{For single-grain models, the gas temperature we take is
relatively low. Table \ref{tab:huv_sigma} shows that when $T_{d} =
T_g$, the remaining amount of H is still reasonably low (intermediate),
contrary to the case where $T_{d} = T_{a}$ (HUV-LH-$\mathrm{T_{a}}$\xspace) which exhibits the
highest amount of H, together with the multi-grain model M-HUV-LH\xspace.
Using only the LH mechanism to form H$_2$ on grains with high temperatures reduces
significantly the amount of $\mathrm{H_2}$\xspace above two scale heights.
This is
particularly true for the M-HUV-LH\xspace Model, where the grain temperature
depends on the grain size: the reservoir of  ``cold'' grains is not large
enough to allow for an efficient $\mathrm{H_2}$\xspace formation on grains in the upper
disk layer. }

\revisionA{On the contrary, near} the midplane, in most models
the grain temperature is always low enough \revisionA{for the LH mechanism
to efficiently form $\mathrm{H_2}$\xspace  and regulate the atomic H abundance.}

\section{Discussion}
\label{sec:discuss}

\revisionA{
We first discuss (Sec.\ref{subsec:gas-phase})
the vertical variations of abundances for CO, CN and
 CS. We then investigate the C,N and O reservoirs  and surface
 chemistry (Sec.\ref{subsec:reservoirs}-\ref{subsec:surface}) before discussing the formation of water and complex organic
 molecules (Sec.\ref{subsec:water}-\ref{subsec:coms}). We also discuss some possible implications of our results
 on planet formation and embryos compositions.  We
 focus the discussion that follows (except for Section\,\ref{subsec:gas-phase})
 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 in Appendix \ref{app:luv}. All results are presented at the final stage of
time evolution i.e. 5.10$^6$ yrs.}

\subsection{Vertical \revisionA{distribution} of abundant (easily observed) molecules} \label{subsec:gas-phase}

\subsubsection{Column densities}
\revisionA{Figure \ref{fig:compare_high} shows the molecular surface densities of
a few popular species detected in TTauri disks.}
\revisionA{Table \ref{tab:huv_sigma} and \ref{tab:luv_sigma} 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, \revisionA{although there can be variations of up to 2.5 orders
of magnitude in the column densities.}}
\revisionA{An analysis of Figures
\ref{fig:100cumul_huv}, \ref{fig:s-100profile_high} and
\ref{fig:m-100profile_high}, which provide the vertical profile of
molecular densities at 100 au, is required to understand the origin
of these variations.}

\revisionA{For the column density of CO it shall be noted that 1) CO is
formed in the gas phase and 2) is sensitive to photodissociation. The
model that produces the most CO is \revisionA{thus LUV-B14-$\mathrm{T_{a}}$\xspace because this} model combines both a low
UV flux, a large column density of H$_2$ that shields CO from the UV
and a high grain temperature that prevents CO from being adsorbed. The
second model \revisionA{producing} the most CO is M-LUV-B14\xspace:  the CO
column density is slightly lower than in LUV-B14-$\mathrm{T_{a}}$\xspace \revisionA{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
M-HUV-LH\xspace is the model that produces the smallest column density:  this
model combines a high UV flux, more penetration due to settling and a
small production of H$_2$ which decreases the shielding.}

\revisionA{For the column densities of CS and CN, Table
\ref{tab:huv_sigma} shows that using predicted molecular surface
densities is not enough to derive general trends.
\revisionA{Nevertheless,} we can note that M-HUV-B14\xspace is the model
that produces the most CS and CN. On the other hand, HUV-B14-$\mathrm{T_{a}}$\xspace produces the
smallest column density of CS while HUV-LH-$\mathrm{T_{a}}$\xspace is the model that produces
the smallest column density of CN.  This suggests routes of formation
which do not simply depend on the grain temperature. }

\revisionA{In the case of the single-grain models (see
Figs.\,\ref{fig:s-100profile_high} and \ref{fig:s-100profile_low}), all
models using the same grain temperature ($T_d = T_g$ or $T_a$)
have
the same vertical profile at altitudes $z < 1.5 H$  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 \revisionA{of the assumed $\mathrm{H_2}$\xspace formation
mechanism}
\revisionA{(see for example, the CO abundance profiles for the two multi-grain models
in Figure \ref{fig:m-100profile_high}, middle panels of columns b-c).
Accordingly, in this relatively high mass disk, differences in column
densities (Tab.\,\ref{tab:huv_sigma}) result from what happens above $z = 1.5 H$.}}

\vspace{0.2cm}

\revisionC{Finally, concerning the intermediate model, Set-HUV-B14-T$\bm{\mathrm{_a}}$, Tab.\,\ref{tab:huv_sigma} shows that
the column densities are rather similar to those found for  model HUV-LH-$\mathrm{T_{a}}$\xspace, with the exception of CS where the column density is closer to
that found for HUV-LH-$\mathrm{T_{g}}$\xspace.  However, these values are integrated column densities and Figs.\,\ref{fig:interm-100profile_LH} and \ref{fig:interm-100profile_B14}
show that the vertical profiles are more complex than the comparison based on
surface densities may suggest.}

  \begin{table*}
\centering
\caption{Column densities [cm$^{-2}$] of main molecules at 100 au for HUV models. Last three columns summarize the main properties of the model with respect to the $\mathrm{H_2}$\xspace formation and grains (\% of grains, in sites, with a temperature above 20 K between 0 and $1.5 H$ and $1.5$ and $2.5 H$). \label{tab:huv_sigma}}
\begin{tabular}{l c c c c c c c c}
\hline
\noalign{\smallskip}
& \textbf{H} & \textbf{H$_2$} & \textbf{CO} & \textbf{CS} & \textbf{CN} & $\mathrm{H_2}$\xspace formation &  0 to $1.5 H$ &  1.5 to $2.5 H$\\
\noalign{\smallskip}
\hline	\hline
\noalign{\smallskip}
\textbf{HUV-LH-T$\bm{\mathrm{_g}}$}  &	 $1.48.10^{21}$ & $8.19.10^{22}$ & $3.24.10^{17}$ & $4.40.10^{12}$ & $3.79.10^{13}$ & "cold" grains &	0\% & 25\% \\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
\textbf{HUV-LH-T$\bm{\mathrm{_a}}$}&        $9.80.10^{21}$ & $7.78.10^{22}$ & $1.01.10^{18}$ & $3.38.10^{12}$ & $3.08.10^{12}$ & more H &100\% & 100\% \\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
\textbf{HUV-B14-T$\bm{\mathrm{_a}}$} &     $3.20.10^{20}$ & $8.25.10^{22}$ & $1.46.10^{18}$ & $1.35.10^{11}$ & $1.34.10^{13}$ & less H & 100\% & 100\%	\\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
\textbf{M-HUV-LH}  &	$8.26.10^{21}$ & $7.86.10^{22}$ & $7.46.10^{16}$ & $6.69.10^{13}$ & $4.62.10^{13}$ & more H &  96\% & 99\%\\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
\textbf{M-HUV-B14}&        $1.31.10^{20}$ & $8.26.10^{22}$ & $9.02.10^{17}$ & $1.92.10^{14}$ & $3.49.10^{14}$ & less H & 96\% & 99\% \\
\noalign{\smallskip}
\hline
\textbf{Set-HUV-B14-T$\bm{\mathrm{_a}}$}&        {\bf $3.65.10^{20}$} & {\bf $8.25.10^{22}$} & {\bf $1.58.10^{18}$} & {\bf $3.96.10^{12}$} &{\bf $3.96.10^{13}$} &  less H & 100\%  & 100\%  \\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
\end{tabular}
\end{table*}
\begin{figure*}
\begin{subfigure}{.48\linewidth}
  \centering

  \includegraphics[width=0.90\linewidth]{figures/compare/cumul/H2100.pdf}
\end{subfigure}
\begin{subfigure}{.48\linewidth}
  \centering

  \includegraphics[width=0.90\linewidth]{figures/compare/cumul/CO100.pdf}
\end{subfigure}\\
\begin{subfigure}{.48\linewidth}
  \centering

  \includegraphics[width=0.90\linewidth]{figures/compare/cumul/CS100.pdf}
\end{subfigure}
\begin{subfigure}{.48\linewidth}
  \centering

  \includegraphics[width=0.90\linewidth]{figures/compare/cumul/CN100.pdf}
\end{subfigure}

\caption{Vertical cumulative surface density [$\mathrm{cm^{-2}}$] of $\mathrm{H_2}$\xspace, CO, CS and CN  at 100 au from the star in the HUV models.}
\label{fig:100cumul_huv}
\end{figure*}

\begin{figure*}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/H2.pdf}

\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Teff/maps/H2.pdf}

\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_B14_Teff/maps/H2.pdf}

\end{subfigure}

\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/CO.pdf}

\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Teff/maps/CO.pdf}

\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_B14_Teff/maps/CO.pdf}

\end{subfigure}

\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/CS.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Teff/maps/CS.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_B14_Teff/maps/CS.pdf}

\end{subfigure}

\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/CN.pdf}
   \subcaption{HUV-LH-$\mathrm{T_{g}}$\xspace}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Teff/maps/CN.pdf}
   \subcaption{HUV-LH-$\mathrm{T_{a}}$\xspace}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_B14_Teff/maps/CN.pdf}
   \subcaption{HUV-B14-$\mathrm{T_{a}}$\xspace}
\end{subfigure}
\caption{Density [cm$^{-3}$] of H$_2$, CO, CS and CN in the gas-phase of the single-grain models in HUV regime. Left column is the HUV-LH-$\mathrm{T_{g}}$\xspace model, middle one is HUV-LH-$\mathrm{T_{a}}$\xspace and right one is HUV-B14-$\mathrm{T_{a}}$\xspace.  Black contours represent the dust temperature (T$_\mathrm{d}$ = T$_\mathrm{g}$ in the left column and T$_\mathrm{d}$ = T$_\mathrm{a}$ in the middle and right columns).}
\label{fig:s-maps-huv}
\end{figure*}

\begin{figure*}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/H2.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_LH/maps/H2.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_B14/maps/H2.pdf}
\end{subfigure}
\label{fig:m-mapsH2}

\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/CO.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_LH/maps/CO.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_B14/maps/CO.pdf}
\end{subfigure}
\label{fig:m-mapsCO}

\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/CS.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_LH/maps/CS.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_B14/maps/CS.pdf}
\end{subfigure}
\label{fig:m-mapsCS}

\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/CN.pdf}
   \subcaption{HUV-LH-$\mathrm{T_{g}}$\xspace}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_LH/maps/CN.pdf}
   \subcaption{M-HUV-LH\xspace}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_B14/maps/CN.pdf}
   \subcaption{M-HUV-B14\xspace}
\end{subfigure}
\label{fig:m-mapsCN}
\caption{Density [cm$^{-3}$] of H$_2$, CO, CS and CN in the gas-phase of HUV-LH-$\mathrm{T_{g}}$\xspace (left column) and of the multi-grain models in HUV regime. Black contours represent the dust temperature (T$_\mathrm{d}$ = T$_\mathrm{g}$ in the left column and T$_\mathrm{d}$ = T$_\mathrm{a}$ in the middle and right columns).}
\label{fig:m-maps}
\end{figure*}

\subsubsection{Vertical variations}

 \paragraph{CO}
\revisionA{ Gas-phase CO abundance is both influenced by grain temperature
and by  abundance of $\mathrm{H_2}$\xspace. }

\revisionA{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.\,\ref{fig:compare_high}). In
M-HUV-LH\xspace, the CO abundance remains low above $z \sim 3 H$ (see
Fig.\,\ref{fig:m-100profile_high}). On the other hand, M-HUV-B14\xspace produces
more $\mathrm{H_2}$\xspace allowing CO to survive at higher altitudes. The resulting CO
column density is of the same order than in HUV-LH-$\mathrm{T_{a}}$\xspace and HUV-B14-$\mathrm{T_{a}}$\xspace. \revisionA{Although
this has little impact on the total column density, the abundance of CO at
$z \leq 1 H$ is much higher}
in 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.}

\paragraph{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 +
$\mathrm{C_2}$  $\rightarrow$ C + CN (69\%) and N + CH $\rightarrow$ H
+ CN (22\%). Note that CH is mostly destroyed by the reaction H + CH
$\rightarrow$ C + $\mathrm{H_2}$. Therefore, more \revisionA{atomic} H in
the upper layers implies less CH available to form CN. The effect is
clearly visible for the single-grain models
(Fig.\ref{fig:s-100profile_high}) where HUV-B14-$\mathrm{T_{a}}$\xspace produces more CN above
3.5H because $\mathrm{H_2}$ is more abundant compared to HUV-LH-$\mathrm{T_{a}}$\xspace. The
effect is less obvious in the multi-grain models
(\ref{fig:m-100profile_high}). Below $z= 1.5\,H$, all the single-grain
models exhibit roughly the same abundance, including HUV-LH-$\mathrm{T_{g}}$\xspace, which
suggests a low sensitivity of CN due to the dust temperature in this range
of values.

 \paragraph{CS}
The main formation process of CS in the upper
layers is a sequence of a bimolecular reaction and a dissociative
recombination reaction: $\mathrm{H_2}$\xspace + $\mathrm{CS^+}$ $\rightarrow$ H +
$\mathrm{HCS^+}$ and $\mathrm{HCS^+}$ + $\mathrm{e^-}$ $\rightarrow$ H + CS.
The sequence shows that the production of CS is dependent on $\mathrm{H_2}$\xspace.
As a consequence, forcing the formation of $\mathrm{H_2}$\xspace 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-$\mathrm{T_{a}}$\xspace
than in HUV-LH-$\mathrm{T_{g}}$\xspace. The production of $\mathrm{H_2}$\xspace being more efficient in the latter
model at high elevation, the peak of CS in HUV-LH-$\mathrm{T_{g}}$\xspace  (Fig.\,\ref{fig:s-100profile_high})
is broader ($\sim 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-$\mathrm{T_{a}}$\xspace, the production of $\mathrm{H_2}$\xspace starts to be efficient at a lower scale
height ($\sim$ 2 versus $3\,H$). Then, CS abundance dramatically
increases at $\sim 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-$\mathrm{T_{a}}$\xspace, the
production of $\mathrm{H_2}$\xspace 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-$\mathrm{T_{g}}$\xspace and HUV-LH-$\mathrm{T_{a}}$\xspace at $4\,H$
(Fig.\ref{fig:s-100profile_high}).

\revisionA{For multi-grain models, we find a \revisionA{larger} CS column
density but the location of the peak of abundance is lower in altitude
(Fig. \ref{fig:100cumul_huv} (c)) because of the relatively low visual
extinction. Just as for the single-grain models, the surface density is
larger in the case of M-HUV-B14\xspace because the formation of $\mathrm{H_2}$ is
more efficient and the quantity of H lower. Below $z = 1.5\,H$, the
difference in abundance is strictly dependent on $T_d$ and as expected,
the hot-grain models (HUV-LH-$\mathrm{T_{a}}$\xspace and HUV-B14-$\mathrm{T_{a}}$\xspace) produce more CS in the gas-phase
than the cold-grain model HUV-LH-$\mathrm{T_{g}}$\xspace, whereas all multi-grain models exhibit
the same vertical profile.}

\revisionC{Concerning the intermediate models (Set-HUV-B14-T$\bm{\mathrm{_a}}$ and Set-HUV-LH-T$\bm{\mathrm{_a}}$), Fig.\,\ref{fig:interm-100profile_LH}  and
Fig.\,\ref{fig:interm-100profile_B14} show that above the altitude of $z \sim 2H$, these models are, as expected, very similar to what is observed for models HUV-LH-$\mathrm{T_{a}}$\xspace because they share the same grain temperature for a number of grain sites which is of the same order. Beyond $z \sim 2H$,
the local gas densities for CO, CN and CS are closer to what is found for the multi-grain models. This is because in the intermediate models the number of grain sites around the midplane is lower than in the case of HUV-LH-$\mathrm{T_{a}}$\xspace models and equal to those of multi-grain models.  Moreover, the amount of UV should be locally higher than in single grain models because of dust settling. Note that in all cases, at the midplane, the abundance of gas species is always very low, well below observational levels.}

 \subsection{\revisionA{Element} reservoirs: C, N, O in the cold midplane}  \label{subsec:reservoirs}

 \begin{figure}
\begin{subfigure}{1\linewidth}
  \centering

 \includegraphics[width=1\linewidth]{figures/compare/reservoirs/carbon2.pdf}

 \label{fig:reservoirs_C}
\end{subfigure}
\begin{subfigure}{1\linewidth}
  \centering

  \includegraphics[width=1\linewidth]{figures/compare/reservoirs/nitrogen2.pdf}

  \label{fig:reservoirs_N}
\end{subfigure}
\begin{subfigure}{1\linewidth}
  \centering

 \includegraphics[width=1\linewidth]{figures/compare/reservoirs/oxygen2.pdf}

 \label{fig:reservoirs_O}
\end{subfigure}
\caption{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.}
 \label{fig:reservoirs}
\end{figure}

The relative importance of reservoirs appearing roughly constant along
the distance from the star, we decide to study the reservoirs at 100
au. All species that have been investigated represent at least 0.5 \%
of the total quantity of carbon, nitrogen or oxygen. As previously
seen, the chemistry below $z = 1.5\,H$ \revisionA{mostly depends on the
grain temperature}. Therefore, it turns out to be sufficient to discuss
HUV-LH-$\mathrm{T_{g}}$\xspace, HUV-LH-$\mathrm{T_{a}}$\xspace, M-HUV-LH\xspace \, \revisionC{and Set-HUV-LH-T$\bm{\mathrm{_a}}$} only.

\revisionA{\paragraph{Carbon bearing species}}
Figure \ref{fig:reservoirs} (top) shows that frozen CO$_2$ is the main
carrier of carbon, both in HUV-LH-$\mathrm{T_{a}}$\xspace and M-HUV-LH\xspace~ \revisionC{(and Set-HUV-LH-T$\bm{\mathrm{_a}}$)}
with 58.1\% and 51.3\% of elemental carbon, respectively, while the quantity of CO$_2$ is
negligible in comparison to the other species in HUV-LH-$\mathrm{T_{g}}$\xspace. The dust temperature
appears too high in these two models to allow for a large quantity of
hydrogen to successively hydrogenate the frozen atomic oxygen to form
water. Instead, the atomic oxygen easily diffuses on the surface and
the relatively high temperature allows to overcome the activation
barrier to rapidly form CO$_2$ by \revisionA{reacting with the frozen CO}. In
HUV-LH-$\mathrm{T_{g}}$\xspace, on the other hand, Fig.\,\ref{fig:reservoirs} (top)
shows that the main carrier of carbon is CH$_4$ with 49.1\%.  In this
colder dust temperature model, CO is initially formed in the gas before
being efficiently and rapidly adsorbed onto the grain surfaces, just
like in hot-grain models. However, once on the surface, CO is channeled
to s-CH$_3$OH through hydrogenation sequences. The photodissociation of
CH$_3$OH then leads to CH$_4$ formation.  Hence, in the midplane,
disks with colder grains will produce more CH$_4$ while warmer disks will produce
more CO$_2$. \revisionA{The other main carrier is the complex organic species
CH$_3$OH, which holds about 10\% of carbon in M-HUV-LH\xspace and HUV-LH-$\mathrm{T_{g}}$\xspace.} \revisionC{Finally, it is important to mention that CO contributes as a reservoir of C for less than  $0.5$\%, in all models.}

\revisionA{\paragraph{Oxygen bearing species}}  CO$_2$ ice is the main
carrier of oxygen in HUV-LH-$\mathrm{T_{a}}$\xspace~\revisionC{(and Set-HUV-LH-T$\bm{\mathrm{_a}}$)} and M-HUV-LH\xspace  with 82.0\% and 72.5\%,
respectively (Fig. \ref{fig:reservoirs}, 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 \revisionA{the atomic oxygen to form water. As described above,
atomic oxygen easily diffuses and also combines with CO to form CO$_2$.}
In HUV-LH-$\mathrm{T_{g}}$\xspace, H$_2$O is by far the largest carrier of oxygen with
76.4\%. Disks with colder grains will
produce more H$_2$O while disks with warmer grains will produce
more CO$_2$. \revisionA{CH$_3$OH is, again, an important carrier and holds about 10\% of oxygen in HUV-LH-$\mathrm{T_{g}}$\xspace.}
\revisionC{As for C, CO holds less than $0.5$\% of Oxygen, in all models.}

\paragraph{Nitrogen bearing species}  Nitrogen is mostly in the form of N$_2$ in M-HUV-LH\xspace (67.8\%) and in HUV-LH-$\mathrm{T_{g}}$\xspace (32.4\%) (Fig. \ref{fig:reservoirs}, middle). Moreover, in almost all models, the gas-phase abundance of N$_2$ is noticeably larger than the one of the other species. The reason for that is a combined effect of the faster conversion of atomic nitrogen into N$_2$ in the gas-phase compared to the depletion of N and the high grain temperature that prevents N$_2$ from depleting onto the surfaces. Indeed, The binding energy of N$_2$ ($E_{\mathrm{b}}(\mathrm{N}_2) = 1100$\,K on amorphous water ice surface) is smaller than most of the other molecules and slightly smaller than that of CO \citep{Minissale+etal_2016, Wakelam+etal_2017} so that N$_2$ is retained more efficiently in the gas. We note that, given the low dust temperature, the gas-phase abundance of N$_2$ in HUV-LH-$\mathrm{T_{g}}$\xspace is the smallest of all models. In HUV-LH-$\mathrm{T_{a}}$\xspace\  \revisionC{and in Set-HUV-LH-T$\bm{\mathrm{_a}}$}, \revisionC{more N$_2$ remains in the gas-phase because of the large dust temperatures.} On the grains, s-NH$_3$ is the main carrier of nitrogen with 42.7\% of the elemental nitrogen. Therefore N$_2$ is expected to be the carrier of nitrogen in the disk midplane but remains more in the gas-phase in hot-grain models. The other main carriers of nitrogen are HCN (29.2\% in HUV-LH-$\mathrm{T_{g}}$\xspace, 12.3\% in HUV-LH-$\mathrm{T_{a}}$\xspace) and CH$_2$NH$_2$ (22.9\% in HUV-LH-$\mathrm{T_{a}}$\xspace, 8.3\% in HUV-LH-$\mathrm{T_{g}}$\xspace).

\subsection{Surface chemistry in the midplane} \label{subsec:surface

\revisionA{One important conclusion of Section \ref{subsec:reservoirs} is
that disks with colder grains will produce more s-CH$_4$ and s-H$_2$O in the midplane
while disks with warmer grains will produce more s-CO$_2$. This clearly demonstrates
the major role of the grain temperature on the surface chemistry in the plane. We
investigate here more deeply the midplane surface chemistry of a few key molecules
for the same four models i.e. HUV-LH-$\mathrm{T_{g}}$\xspace, HUV-LH-$\mathrm{T_{a}}$\xspace, M-HUV-LH\xspace}\revisionC{, and Set-HUV-LH-T$_\mathrm{a}$}.
\revisionC{Another important conclusion is that our intermediate models, Set-HUV-LH-T$\bm{\mathrm{_a}}$ and Set-HUV-B14-T$\bm{\mathrm{_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-$\mathrm{T_{a}}$\xspace in terms of main reservoirs in the midplane. Small differences occur for relatively
minor constituents (s-CO, s-H$_2$CO, s-CH$_3$OH, s-HCN for example) because of the
different effective dust area in the two cases (Single vs Set), and hence different timescales for conversion on grain surfaces.}

\revisionA{Figures \ref{fig:s-s-maps-huv} and \ref{fig:m-s-maps-huv}
present maps of number density for s-H$_2$, s-CO, s-CS and s-CN locked on grain
surfaces in the case of single-grain and multi-grain models,
respectively. Figure \ref{fig:ab_surface} presents the surface abundance
per total atomic hydrogen for various molecules as a function of the grain radius $a$
and their temperature in the midplane of the disk. The crosses, the
square markers, and the round markers stand for the multi-grain model
at 30 au, 100 au and 200 au respectively. The triangle pointing upward
represents the surface abundance in the single-grain model HUV-LH-$\mathrm{T_{g}}$\xspace and
the triangle pointing downward represents the abundance in the
single-grain model HUV-LH-$\mathrm{T_{a}}$\xspace, both at 100 au. }

\subsubsection{Carbon and oxygen-bearing molecules: CO, CO$_2$ and CH$_4$}
CO is formed in the gas-phase prior to being accreted on small grains
during the first 10$^6$ yrs of the simulation. However, because of its
small binding energy \citep[$E_{\mathrm{b}}$(CO) = 1300 K][]{Wakelam+etal_2017},
it is rapidly thermally desorbed from the hot ($T_g > 20$\,K)
smaller grains and gets accreted by large cold grains ($T_g < 20$\,K), \revisionA{where}
it remains locked until the simulation is
stopped. Considering this mechanism, CO is distributed according to the
grain population and Fig.\,\ref{fig:ab_surface} (top left panel) shows
two regimes of abundance, accounting for the desorption barrier, where CO ice tends to stay locked on grains of size $\gtrsim$
1 $\mu$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-$\mathrm{T_{a}}$\xspace and a high abundance in the cold-grain model
HUV-LH-$\mathrm{T_{g}}$\xspace.

As for CO$_2$, its binding energy is relatively high
\citep[$E_{\mathrm{b}}(\mathrm{CO}_2) = 2600$ K][]{Wakelam+etal_2017}.
As a consequence the abundance of s-CO$_2$ depends very little on the
dust temperature {and simply depends on the dust surface area instead}. s-CO$_2$ is efficiently formed
through
\begin{gather}
\label{eq:CO2b}
\mathrm{s\mhyphen O + s\mhyphen CO \rightarrow s\mhyphen CO_2 \: (E_A = 1000\,K)}\\
\label{eq:CO2d}
\mathrm{s\mhyphen OH + s\mhyphen CO \rightarrow s\mhyphen CO_2 + s\mhyphen H}.
\end{gather}
\noindent Therefore, s-CO$_2$ is more abundant on the smaller grains because their
temperature allows to overcome the activation barrier of reaction
\ref{eq:CO2b} and the drop in abundance on large grains at 200 au (and
further) originates from the grain temperature getting too low for
reaction \ref{eq:CO2b} to be activated. Therefore, the outer disk forms
fewer CO$_2$ ice.

In cold-grain models, s-CH$_4$ is the largest carrier of carbon (see
Section \ref{subsec:reservoirs}). As expected, Figure \ref{fig:ab_surface}
(Bottom-left) shows that s-CH$_4$ is more abundant on cold grains ($T_g < 20$\,K).
s-CH$_4$ is formed via the hydrogenation of s-CH$_3$. The latter is
formed through the following two sequences:
\begin{gather}
\label{eqn:CH4a}
\mathrm{s\mhyphen CO + h\nu \rightarrow s\mhyphen C + s\mhyphen O} \\
\mathrm{s\mhyphen H_2 + s\mhyphen C \rightarrow s\mhyphen CH_2} \\
\mathrm{s\mhyphen H_2 + s\mhyphen CH_2 \rightarrow s\mhyphen CH_3 + s\mhyphen H}
\end{gather}
\noindent and
\begin{gather}
\label{eqn:CH4b}
\mathrm{s\mhyphen CH_3OH + h\nu \rightarrow s\mhyphen OH + s\mhyphen CH_3}.
\end{gather}

\noindent In view of the above, s-CH$_3$OH is formed on the grains prior
to being photodissociated (if the computation domain is sufficiently
extended). Then the released atomic carbon (Reaction \ref{eqn:CH4a})
and s-CH$_3$ (Reaction \ref{eqn:CH4b}) lead to s-CH$_4$ formation.
Consequently, the formation of s-CH$_4$ is more efficient on large cold
grains since 1) CO must be adsorbed and 2) hydrogenation is needed to
create s-CH$_3$OH. s-CH$_4$ ice is thus found to be more abundant on large
grains and in the outer disk.

  \subsubsection{Sulfur-bearing species: H$_2$S}
In the gas phase, H$_2$S has been identified in dense cloud cores,
cometary comae and \citet{Phuong+etal_2018} reported the first
detection of  H$_2$S in the cold and dense ring about the TTauri star
GG Tau A. The low H$_2$S column densities observed in GG Tau A and in
dense molecular clouds are assumed to be the result of a strong sulfur
depletion \citep{Hudson+Gerakines_2018, Phuong+etal_2018}.

We investigate here the formation routes of S-H bonds and more
specifically the simple S-bearing molecule Hydrogen sulfide, H$_2$S, in
the midplane icy mantles of our multi-grain models.
Figure \ref{fig:ab_surface} (second row, right column) shows the abundance
of s-H$_2$S as a function of the grain size. There is a strong
anti-correlation with grain temperatures. The smallest s-H$_2$S
abundances are located on the hottest grains ($T_d \sim 27$\,K) while the largest
abundances are found on large colder grains ($T_d \lessapprox 15$\,K). This anti-correlation is
observed at all radii. The abundances on the large grains ($> 10 \mathrm{\mu m}$)
are independent of the radius. For the small grains,
however, the inner disk exhibits larger abundances than the outer disk
because of higher collision rates in the inner regions. The strong
anti-correlation with temperature is directly related to the main
formation pathways of H$_2$S that implies hydrogenation sequences:

\begin{gather}
\label{eqn:H2Sa}
\mathrm{s\mhyphen S \xrightarrow{s\mhyphen H} s\mhyphen HS \xrightarrow{s\mhyphen H} s\mhyphen H_2S} \\
\label{eqn:H2Sb}
\mathrm{s\mhyphen CH_3SH \xrightarrow{s\mhyphen H} s\mhyphen CH_3 + s\mhyphen H_2S}.
\end{gather}

\noindent We note that the channels to s-CH$_3$SH also involve
successive hydrogenations (of s-CS, mainly). Reaction \ref{eqn:H2Sa}
occurs on grains of all sizes but is more effective on the large cold
grains ($T_d \lessapprox 15$\,K) and Reaction \ref{eqn:H2Sb} is only found to be effective on the
colder grains.

\subsubsection{Nitrogen-bearing species: NH$_3$}

The main nitrogen-bearing icy species in the midplane are s-NH$_3$, s-HCN,
s-N$_2$ and s-CN at all radii. In the multi-grain models, as seen in
Fig.\,\ref{fig:ab_surface}, the abundance of s-NH$_3$ is inversely
proportional to the grain size until $\approx 1 \mathrm{\mu m}$ where the
abundance stops decreasing. The abundances of s-NH$_3$ are roughly
the same at 100 and 200 au but not at 30 au.

The main pathway to s-NH$_3$ originates from the accretion of N at the
surface of the grains and its successive hydrogenations to form s-NH$_3$
\citep[see][]{Aikawa+etal_2015, Eistrup+etal_2018, Ruaud+Gorti_2019}.
Such hydrogenation reactions are efficient on cold surfaces ($T_d \lessapprox 15$\,K), which
explains why the abundance of s-NH$_3$ does not keep dropping
when grain temperatures are low. As for the  warmer small grain surfaces,
their high quantity implies a high collision rate with NH$_3$ which
stays on their surface given its binding energy
(E$_{\mathrm{b}}$(NH$_3$) = 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$_2$CO:

\begin{equation}
\label{eq:nh3
	  \mathrm{s\mhyphen NH_2CO \xrightarrow{s\mhyphen H} s\mhyphen NH_3 + s\mhyphen CO}.
\end{equation}

\begin{figure*}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/JKH2.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Teff/maps/JKH2.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_B14_Teff/maps/JKH2.pdf}
\end{subfigure}

\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/JKCO.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Teff/maps/JKCO.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_B14_Teff/maps/JKCO.pdf}
\end{subfigure}

\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/JKCS.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Teff/maps/JKCS.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_B14_Teff/maps/JKCS.pdf}
\end{subfigure}

\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/JKCN.pdf}
   \subcaption{HUV-LH-$\mathrm{T_{g}}$\xspace}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Teff/maps/JKCN.pdf}
   \subcaption{HUV-LH-$\mathrm{T_{a}}$\xspace}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_B14_Teff/maps/JKCN.pdf}
   \subcaption{HUV-B14-$\mathrm{T_{a}}$\xspace}
\end{subfigure}
\caption{Density [cm$^{-3}$] of H$_2$, CO, CS and CN in the gas-phase of
the single-grain models in HUV regime. Left column is the HUV-LH-$\mathrm{T_{g}}$\xspace model,
middle one is HUV-LH-$\mathrm{T_{a}}$\xspace and right one is HUV-B14-$\mathrm{T_{a}}$\xspace. Black contours represent the dust temperature (T$_\mathrm{d}$ = T$_\mathrm{g}$ in the left column and T$_\mathrm{d}$ = T$_\mathrm{a}$ in the middle and right columns).}
\label{fig:s-s-maps-huv}
\end{figure*}

\begin{figure*}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/JKH2.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_LH/maps/JKH2.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_B14/maps/JKH2.pdf}
\end{subfigure}
\label{fig:m-surf-mapsH2}

\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/JKCO.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_LH/maps/JKCO.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_B14/maps/JKCO.pdf}
\end{subfigure}
\label{fig:m-surf-mapsCO}

\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/JKCS.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_LH/maps/JKCS.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_B14/maps/JKCS.pdf}
\end{subfigure}
\label{fig:m-surf-mapsCS}

\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/JKCN.pdf}
   \subcaption{HUV-LH-$\mathrm{T_{g}}$\xspace}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_LH/maps/JKCN.pdf}
   \subcaption{M-HUV-LH\xspace}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_B14/maps/JKCN.pdf}
   \subcaption{M-HUV-B14\xspace}
\end{subfigure}
\label{fig:m-surf-mapsCN}
\caption{Density [cm$^{-3}$] of H$_2$, CO, CS and CN in the gas-phase of HUV-LH-$\mathrm{T_{g}}$\xspace (left column) and of the multi-grain models in HUV regime. Black contours represent the dust temperature (T$_\mathrm{d}$ = T$_\mathrm{g}$ in the left column and T$_\mathrm{d}$ = T$_\mathrm{a}$ in the middle and right columns).}
\label{fig:m-s-maps-huv}
\end{figure*}

\begin{figure*}
  \centering

  \includegraphics[width=1.0\linewidth]{figures/compare/temp_surface/full2.pdf}
\caption{Grain surface abundance per total atomic hydrogen \revisionC{and per grain size bin} for various
molecules at the final stage of integration (i.e. 5\,10$^6$ yrs) as a
function of the grain radius $a$ in relation to the grain temperatures.
Cross markers stand for the multi-grain model at 30 au, the square
markers stand for the multi-grain model at 100 au and the round markers
stand for 200 au. The triangles pointing upward represent the surface
abundance in the single-grain model HUV-LH-$\mathrm{T_{g}}$\xspace (grains of size 0.1
$\mathrm{\mu m}$) and the triangle pointing downward represents the
abundance in the single-grain model HUV-LH-$\mathrm{T_{a}}$\xspace, both at 100 au.  \revisionC{Abundances
being per grain size bin, it is the sum of the abundances in the multi-grain models
that should be compared to the abundances in the single grain case.}}
\label{fig:ab_surface}
\end{figure*}

\subsection{Water} \label{subsec:water}

The formation of water is key to explain the evolution of
protoplanetary disks and formation of comets, both decisive elements to
understand the delivery of water to planetary surfaces. Spectral lines
of H$_2$O in the gas-phase of inner disk regions have now been widely
detected in disks around TTauri stars using near- and mid-infrared
\citep{Carr+Najita_2008, Pontoppidan+etal_2010} or far-infrared
observations \citep{Riviere-Marichalar+etal_2012}. Cold water in the
outer disk regions, on the other hand, is harder to detect due to low
excitation lines.

It is assumed \citep{Dominik+etal_2005, lecar+etal_2006,
Podio+etal_2013, Du+Bergin_2014} that water on grain surfaces cannot
exist inside the snowline ($T_{d} \gtrsim 150$\,K), which represents radii
smaller than a few to ten au around T\,Tauri stars. The innermost radius
computed by our models is 30 au, which is beyond the radial ice line of
water. Therefore, we do not discuss the radial location of the water
ice line.

\subsubsection{Gas-phase}
\revisionA{Three main parameters control  the abundance of water in
gas-phase, the photodesorption by FUV photons and cosmic rays, the
dust temperature and the chemical reactions in the gas-phase, in
particular: }

\revisionA{ \begin{gather}
\label{eqn:h2od}
\mathrm{H_2 + OH^+ \rightarrow H + H_2O^+} \\
\mathrm{H_2O^+ + H_2 \rightarrow H + H_3O^+} \\
\mathrm{H_3O^+ + e^-  \rightarrow H + H_2O}
\end{gather}}

\revisionA{Thus $\mathrm{H_2}$\xspace and OH$^+$ are the main precursors for H$_2$O. Whether H$_2$ or OH$^+$ is the limiting factor is determined by this other sequence:}

\revisionA{ \begin{gather}
\label{eqn:oh}
\mathrm{H_2 + CR \rightarrow H_2^+ + e^-}\\
\label{eqn:oh2}
\mathrm{H_2^+ + H_2 \rightarrow H_3^+ + H} \\
\label{eqn:oh3}
\mathrm{H_3^+ + O  \rightarrow OH^+ + H_2}.
\end{gather}}

\revisionA{\noindent These two sequences show that H$_2$O formation
depends on $\mathrm{H_2}$\xspace and OH$^+$ but the latter also depends on $\mathrm{H_2}$\xspace.
Therefore, this is the $\mathrm{H_2}$\xspace abundance that sets the gas-phase formation
of water and accounts for the difference between the models (not
cosmic-rays or OH$^+$). However, the gas-phase formation of H$_2$O
only represents a small fraction of all H$_2$O molecules formed as most
H$_2$O are formed on grain surfaces. The abundance of gas-phase water
is thus mainly governed by both the formation rate of icy water and by
the photodesorption rate. Thermal-desorption is ineffective in most
parts of the disk because the binding energy of water is high. The
formation rate depends on the grain temperature while the
photodesorption rate depends on the FUV flux and cosmic-rays. The
analysis of the results herein (Figs.\,\ref{fig:s-water_maps} and
\ref{fig:m-water_maps}) shows that the dust temperature plays the major
role in the production of gas-phase water.}

\revisionA{\paragraph{Single-grain models} The vertical snowline forms a
clear boundary between high and low gas-phase H$_2$O abundance at $z/r
\sim 0.15$ (Fig.\,\ref{fig:s-water_maps}) and we see that HUV-B14-$\mathrm{T_{a}}$\xspace 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.\,\ref{fig:struct} (bottom-right panel). The cold-grain model HUV-LH-$\mathrm{T_{g}}$\xspace
produces substantially more water in the upper layers than the other
single-grain models,
because low dust temperature
increase the hydrogenation rates on the grain surfaces. The
photodesorption rates are approximately equal in the three models
(same flux) but since HUV-LH-$\mathrm{T_{g}}$\xspace 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-$\mathrm{T_{g}}$\xspace than in the
other two models. On the other hand, the B14's prescription increases
only marginally the column density of gas-phase water.}

\revisionA{ \paragraph{Multi-grain models} The vertical snowline is far
less apparent (Fig.\,\ref{fig:m-water_maps}, top row, middle and right
column). Both models have approximately the same distribution of
abundance. because both have the same dust temperature distribution and
same UV. We notice, however, that B14's prescription in M-HUV-B14\xspace produces
marginally more water in the upper layer of the disk than in
M-HUV-LH\xspace, due to the gas-phase formation of water via the two sequences
presented above.}

\subsubsection{Water ice}

\revisionA{ \paragraph{Single-grain models} The vertical snowline location
is clearly defined (see Fig.\,\ref{fig:s-water_maps}, bottom row),
although we note the presence of a notable band around $z = 3\,H$ ($z =
2\,H$) in HUV-LH-$\mathrm{T_{g}}$\xspace  (resp. HUV-LH-$\mathrm{T_{a}}$\xspace). The reason is the following activated
hydrogenation reaction: }

\begin{gather}
\label{eqn:h2oc}
\mathrm{s\mhyphen H + s\mhyphen HOOH \rightarrow s\mhyphen H_2O + s\mhyphen OH \: (E_A = 1400\,K)}.
\end{gather}

\noindent Indeed, the band illustrates the location in the disk where
this reaction is activated. Below the band, the reaction is not
activated and water ice is mostly photodesorbed while above the band
the photodesorption rates competes with the formation rate and the
abundance drops. In HUV-LH-$\mathrm{T_{a}}$\xspace (Fig. \ref{fig:s-water_maps}, middle
column), Reaction \ref{eqn:h2oc} 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-$\mathrm{T_{a}}$\xspace 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.

\revisionA{ \paragraph{Multi-grain models} The same mechanisms are at work
(see middle and right columns of Fig.\ref{fig:m-water_maps}). Yet, the
dust temperature distribution involves that we find both hot (up to $\sim$ 95 K) and cold (as low as $\sim$ 7 K)
grains nearly everywhere in the disk. Consequently, Reaction
\ref{eqn:h2oc} is activated at lower altitudes and the formation of
s-H$_2$O through hydrogenation can occur at higher altitudes on cold
medium- and large-sized grains so that the snowline is located at a
higher altitude. }

\revisionA{\paragraph{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-$\mathrm{T_{g}}$\xspace than in the two others because, again,
the dust temperature is smaller. }
As presented in Fig.\,\ref{fig:ab_surface} (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 $\mu m$) is larger because the grain
temperature is low enough to allow for effective hydrogenation.

\begin{figure*}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/H2O_ab.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Teff/maps/H2O_ab.pdf}
 \end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_B14_Teff/maps/H2O_ab.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/sH2O_ab.pdf}
  \subcaption{HUV-LH-$\mathrm{T_{g}}$\xspace}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Teff/maps/sH2O_ab.pdf}
  \subcaption{HUV-LH-$\mathrm{T_{a}}$\xspace}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_B14_Teff/maps/sH2O_ab.pdf}
  \subcaption{HUV-B14-$\mathrm{T_{a}}$\xspace}
\end{subfigure}

\caption{Water abundance in the gas-phase (top) and on the grain surfaces (bottom) of the single-grain models. Left column shows the abundances in HUV-LH-$\mathrm{T_{g}}$\xspace, middle in HUV-LH-$\mathrm{T_{a}}$\xspace and right one in HUV-B14-$\mathrm{T_{a}}$\xspace. In the top row, the solid black lines show the dust temperature isocontours, in the bottom row they denote isocontours of scale heights.}
\label{fig:s-water_maps}
\end{figure*}

\begin{figure*}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/H2O_ab.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_LH/maps/H2O_ab.pdf}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_B14/maps/H2O_ab.pdf}
 \end{subfigure}

 \begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/SINGLE/HUV_HL_Tg/maps/sH2O_ab.pdf}
  \subcaption{HUV-LH-$\mathrm{T_{g}}$\xspace}
\end{subfigure}
 \begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_LH/maps/sH2O_ab.pdf}
  \subcaption{M-HUV-LH\xspace}
\end{subfigure}
\begin{subfigure}{.33\linewidth}
  \centering

  \includegraphics[width=1.05\linewidth]{figures/MULTI/HUV_B14/maps/sH2O_ab.pdf}
  \subcaption{M-HUV-B14\xspace}
\end{subfigure}

\caption{Water abundance in the gas-phase (top) and on the grain surfaces (bottom). Left column shows the abundances in HUV-LH-$\mathrm{T_{g}}$\xspace and middle and right column in the multi-grain models. In the top row, the solid black lines show the dust temperature isocontours, in the bottom row they denote isocontours of scale heights.}
\label{fig:m-water_maps}
\end{figure*}

  \subsection{COMs}  \label{subsec:coms}

Complex organic molecules (COMs), supposedly constituting the crucial bond between simple ISM molecules and prebiotic chemistry, are hard to detect due to their low abundance and weak emission lines originating from their molecular complexity. This is especially true for objects with small spatial extensions like disks. Simple organic molecules have been regularly detected over the last decades such as H$_2$CO \citep{Dutrey+etal_1997, Aikawa+etal_2003, Oberg+etal_2011}, HC$_3$N \citep{Chapillon+etal_2012}, CH$_3$CN \citep{Oberg+etal_2015} and  CH$_3$OH \citep{Walsh+etal_2016}. We investigate here  the effectiveness of our models to synthesize COMs from surfaces or gas-phase when given sufficient time ($5\times10^6$ yrs).

 \begin{figure}
\includegraphics[width=\linewidth]{figures/MULTI/COMS/ch3cn_compare.pdf}
\caption{Altitude above the midplane as a function of the abundance relative to H of CH$_3$CN at 100 au. Blue line shows the abundance produced by model M-HUV-LH\xspace, red line shows the abundance produced by model M-HUV-B14\xspace. We see that the abundance dramatically drops above $\sim$ 2 scale heights in M-HUV-LH\xspace while M-HUV-B14\xspace keeps a rather high abundance at high altitudes. This is due to the available quantity of gas-phase H$_2$ in M-HUV-B14\xspace, which enhances the production of CH$_3^+$.
\label{fig:methyl_compare}}
\end{figure}

  \subsubsection{Methyl Cyanide (CH$_3$CN)}

\revisionA{C-N bonded species are very important for the synthesis of
Glycine \citep{Goldman+etal_2010}. \citet{Oberg+etal_2015} have
observed methyl cyanide (CH$_3$CN) in the disk around the Herbig Ae star
MWC 480, a star exhibiting a stronger UV flux and a warmer disk than a
TTauri star. }

\revisionA{\paragraph{CH$_3$CN to HCN abundance ratio}
To facilitate comparisons,  we discuss the abundance of CH$_3$CN
relative to HCN, as practiced in cometary studies
\citep{Mumma+Charnley_2011, Oberg+etal_2015}. }The model that produces
the largest CH$_3$CN/HCN ratio is the single-grain model HUV-LH-$\mathrm{T_{g}}$\xspace with
$\sim$ 14\% (vertically integrated) at 100, 120 and 140 au while the
other single-grain models produce nearly more than $\sim$ 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-$\mathrm{T_{g}}$\xspace. In
multi-grain models, the gas-phase abundance ratio never exceeds 1\%.
However, it never goes below 0.1\%, implying a rather uniform
production of CH$_3$CN along the radius as compared to the single-grain
models.  M-HUV-B14\xspace produces slightly more CH$_3$CN than M-HUV-LH\xspace. \revisionA{ The
model HUV-LH-$\mathrm{T_{g}}$\xspace (with the coldest grains overall) is the one that globally
produces the largest CH$_3$CN/HCN ratio.}

\revisionA{\paragraph{CH$_3$CN in gas-phase}
In the upper regions,} the main known pathways to form gas-phase
CH$_3$CN are the reaction of HCN with CH$_3^+$ and the photodesorption
of icy CH$_3$CN (thermal desorption is not effective because of the
large binding energy of CH$_3$CN). Results show that the reaction
between HCN and CH$_3^+$  is the most effective route in M-HUV-B14\xspace while
only photodesorption is effective in M-HUV-LH\xspace. The reason is
straightforward since the production of CH$_3^+$ requires the presence
of molecular hydrogen:

\begin{gather}
\label{eqn:ch3+}
\mathrm{C^+ \xrightarrow{H_2} CH_2^+ \xrightarrow{H_2} H + CH_3^+}
\end{gather}

\noindent and H$_2$ is much more abundant in M-HUV-B14\xspace in the upper layers.
At 100 au, the abundance of CH$_3^+$ in M-HUV-LH\xspace is very low ($\sim
10^{-20}$) as compared to the abundance in M-HUV-B14\xspace ($\sim 10^{-10}$). It
follows that, as seen in Fig.\,\ref{fig:methyl_compare}, M-HUV-B14\xspace produces
much more CH$_3$CN than in M-HUV-LH\xspace. Therefore, the presence of H$_2$
appears to be decisive for the production of gas-phase CH$_3$CN in the
upper regions of the disk. However, this difference between the two
models does not significantly impact the total column density of
CH$_3$CN as most of it is produced near the midplane.

\revisionA{\paragraph{CH$_3$CN ices} Figure \ref{fig:surfdensmap}}
(bottom) shows the abundance of s-CH$_3$CN ices in the midplane relative
to the grain size as a function of the radius. s-CH$_3$CN mostly stays on
warm small grains while there is a large drop in abundance around grain
populations of size $\sim$ 0.317 $\mathrm{\mu m}$. The results show
that the two main reactions that lead to s-CH$_3$CN are the
adsorption of gas-phase CH$_3$CN onto the surfaces and the following
surface reaction:

\begin{gather}
\label{eqn:ch3cn}
\mathrm{s\mhyphen CH_3 + s\mhyphen CN \rightarrow s\mhyphen CH_3CN}.
\end{gather}

\noindent \revisionA{s-CH$_3$CN is more abundant on small grains because of
three main reasons.} The first reason is the large number density of
small grains that implies large collision rates with 1) gas-phase
produced CH$_3$CN and 2) with CH$_3$ that, despite a low binding
energy, will briefly stick and quickly react with s-CN (see reaction
\ref{eqn:ch3cn}). The second reason is the high binding energy of
CH$_3$CN \citep[E$_\mathrm{b}$(CH$_3$CN) = 4680 K][]{Wakelam+etal_2017}
and of the precursors (s-CN) so that they stay on small grain surfaces.
Moreover, the temperature of small grains is sufficiently high to allow
precursors to diffuse on the surface and increase reaction rates. The
last reason is that even if s-CH$_3$CN ices can be formed on large grain
surfaces, photodesorption rates are large enough for CH$_3$CN to be
desorbed from large grains and seed small grains as the disk evolves.

\subsubsection{Methanol and formaldehyde}

\begin{figure}
\includegraphics[width=\linewidth]{figures/MULTI/COMS/ch3oh_compare.pdf}
\caption{Altitude above the midplane as a function of the abundance relative to H of CH$_3$OH at 100 au. Blue line shows the abundance produced by model M-HUV-LH\xspace, red line shows the abundance produced by model M-HUV-B14\xspace. We see that the abundance dramatically drops at $\sim$ 2 scale heights in both models, although M-HUV-B14\xspace maintains a larger abundance than M-HUV-LH\xspace at high altitudes.
\label{fig:CH3OH_compare}}
\end{figure}

\paragraph{CH$_3$OH in gas-phase} For multi-grain models,
Fig.\,\ref{fig:CH3OH_compare} shows that the vertical profile of
gas-phase CH$_3$OH is divided into two layers separated by an abrupt
drop in abundance at $\sim 2\,H$, which corresponds roughly to the dust
temperature transition (Fig. \ref{fig:weighted_100}). There is a peak
in gas-phase CH$_3$OH at $\sim 1.8\,H$ in both models. This corresponds
to the location where dust temperature is low enough for adsorbed CO to
be hydrogenated and where the UV flux is high enough to efficiently
photodissociate a significant fraction of newly formed methanol. The
lower layer ($ < 2\,H$) contains most of CH$_3$OH column density. In the
upper layer, M-HUV-B14\xspace produces more CH$_3$OH than M-HUV-LH\xspace.

\paragraph{CH$_3$OH to H$_2$CO abundance ratio} In multi-grain models,
the ratio CH$_3$OH/H$_2$CO increases by a factor of $\sim$ 2 to 3 in
the outermost region ($\geq$ 200 au) between $10^6$ and $5\times10^6$
years while it remains stable in the inner region (< 100 au). In M-HUV-B14\xspace,
after $10^6$ years, the averaged ratio in the region 80 - 120 au is
equal to 1.3. This appears to be consistent with the observed ratio in
TW Hya \citep{Walsh+etal_2016, Carney+etal_2019}. In M-HUV-LH\xspace, the
averaged ratio is equal to 2.4. The difference is due to a more
efficient pathway to form H$_2$CO in the upper disk regions of M-HUV-B14\xspace.
Gas-phase H$_2$CO is efficiently formed through the following sequence:
\begin{gather}
\label{eqn:h2co}
\mathrm{H_2 +CH_3^+ \rightarrow CH_5+} \\
\mathrm{CH_5^+ + e^- \rightarrow H + H + CH_3/H_2 + CH_3} \\
\mathrm{O + CH_3 \rightarrow H + H_2CO}.
\end{gather}

\noindent
This shows that molecular hydrogen is a key element in the formation
process of gas-phase formaldehyde and M-HUV-B14\xspace provides much more H$_2$
than M-HUV-LH\xspace.

\paragraph{Grain size dependence of surface chemistry for  CH$_3$OH}

\revisionA{Figure \ref{fig:surfdensmap} gives a set of maps of ice abundance
in the midplane with grain sizes as a function of the computed radii.
Top left panel is the grain temperature which allows to analyze the
surface abundance in view of the surface temperature. We have same
results for both multi-grain models. On icy mantles, both s-H$_2$CO and
s-CH$_3$OH are formed via the same hydrogenation sequence:}

\revisionA{
\begin{gather}
\label{eqn:hydro}
\mathrm{s\mhyphen CO \xrightarrow{H} s\mhyphen HCO \xrightarrow{H} s\mhyphen H_2CO \xrightarrow{H} s\mhyphen CH_3O \xrightarrow{H} s\mhyphen CH_3OH}.
\end{gather}
}
\noindent However, we note that H$_2$CO is mainly formed in the
gas-phase. As seen in Fig.\,\ref{fig:surfdensmap} (top right), CO tends
to be efficiently adsorbed onto cold surfaces (T$_d$ $\lessapprox$ 17 K). The map of H$_2$CO
(bottom left) grossly follows the same distribution although the trend
is globally less pronounced. CH$_3$OH, on the other hand, exhibits a
totally different distribution and tend to concentrate on small grain
surfaces (bottom right). The main reason comes from the fact that the
binding energy increases as the molecules become more complex
\citep[e.g., $E_{\mathrm{b}}$(CO) $\sim$1200 K,
$E_{\mathrm{b}}$(H$_2$CO) $\sim$3200 K, $E_{\mathrm{b}}$(CH$_3$OH)
$\sim$5500 K][]{Minissale+etal_2016, Noble+etal_2012,
Collings+etal_2004} and that the hydrogenation of s-H$_2$CO (which leads
to s-CH$_3$OH) has an activation barrier and cannot efficiently occur on
cold grain surfaces.

\begin{figure*}
\begin{subfigure}{.49\linewidth}
  \centering
  \includegraphics[width=1.0\linewidth]{figures/MULTI/COMS/temperature.pdf}

\end{subfigure}
\begin{subfigure}{.49\linewidth}
  \centering

  \includegraphics[width=1.0\linewidth]{figures/MULTI/COMS/CO.pdf}
\end{subfigure}
\begin{subfigure}{.49\linewidth}
  \centering

  \includegraphics[width=1.0\linewidth]{figures/MULTI/COMS/H2CO.pdf}
\end{subfigure}
\begin{subfigure}{.49\linewidth}
  \centering

  \includegraphics[width=1.0\linewidth]{figures/MULTI/COMS/CH3OH.pdf}
\end{subfigure}
\begin{subfigure}{.49\linewidth}
  \centering

  \includegraphics[width=1.0\linewidth]{figures/MULTI/COMS/CH3CN.pdf}
\end{subfigure}

\caption{multi-grain models: Maps of ice abundance relative to H in the
midplane, with grain sizes as a function of the radius. Top left panel
is the distribution in grain temperature. We see that the hottest
region is located at around 60 au from the star on grains of size 0.148
$\mathrm{\mu m}$.}
\label{fig:surfdensmap}
\end{figure*}

\subsection{Impact for planet formation}

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

For a T\,Tauri disk of a few Myr, radial changes of the millimeter
spectral index \citep[e.g.][]{Birnstiel+etal_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 which
is similar to that observed in the ISM, suggesting a similar size
distribution, even if settling still occurs with larger grains residing
onto the midplane. As a consequence, toward the midplane, one
expects to observe different icy layers on grains in the inner (large
cold grains) and outer disk (warmer smaller grains). We have seen in
the previous section that  s-CH$_4$ and s-H$_2$O are preferentially located
on large grains while s-CO$_2$ remains on smaller ones (Fig.\,\ref{fig:ab_surface}).
This implies also a different chemical coupling
between the gas and solid phase. Moreover, we have seen that complex
molecules tend to stick and remain on small warmer grains ($\leq$ $\sim$
0.1 $\mathrm{\mu m}$), which are less likely to drift toward the inner
disk. The outer disk (r $\geq$ 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.

\revisionA{Grain-grain collisions may modify this picture.  The characteristic timescale (the time for a large grain to encounter
small grains of equivalent area) around 100 au is of order $10^4$--$10^5$ years. However, while large grains may accrete the smaller ones, and
inherit part of their surface composition while growing, small grains are re-created by disruptive events that may severely affect
their surface layer and ice mantles. The outcome of the surface state during such grain-grain collisions remains to be studied.}

\revisionA{Nevertheless, in all cases, the grain size distribution, because it results in a spread of dust temperature, allows for a
faster depletion of C and O, that are converted to s-CO$_2$ on the colder grains in the disk midplane, but also a more
active surface chemistry leading to formation of COMs on the warmer grains in the upper layer. Such a behavior cannot
be reproduced by a single grain temperature chemical model, because the formation processes are highly non-linear,
since the mobility and evaporation efficiency depend exponentially on dust temperature. }

\revisionA{
The range in  dust temperatures also spreads out spatially the snowlines
that become fuzzy. The spatial extent of these ``snow bands'' will depend
on the specific dust properties and disk mass. The temperature spread is
larger for less massive disks.} %

As disks evolve, the amount of dust (and gas) decreases until they
reach the optically thin, debris disk phase. A less massive dust disk
necessarily implies a warmer dust disk, which should affect the
composition of the icy layers at the grain surface in evolved disks.
Recent surveys of 5-40\,Myr old disks suggest the existence of an
intermediate phase (at least for Herbig Ae stars) between debris and
protoplanetary disks \citep[e.g.][]{Kospal+etal_2013,
Pericaud+etal_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, \citet{Difolco+etal_2020} observed
a CO surface density which is only a factor 10 lower than in a typical
protoplanetary HAe disk, whereas the amount of millimeter grains is a
factor 50 times smaller than in typical young disks. The reduced amount
of grains, and in some cases at least possible changes in the grain
size distribution \citep[e.g.][]{hughes+etal_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.

\section{Summary}.
\label{sec:concl}
We presented here a new chemical model of protoplanetary disks that
consistently accounts for the impact of grain sizes on the dust
temperatures and UV penetration. Representative results were obtained
for a typical disk around T Tauri stars, using parameters derived from
the Flying Saucer observations \citep{Dutrey+etal_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, method
that accounts for the frequency dependence of photorates.
\revisionA{
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. }

\begin{itemize}
\item We identify the formation of H$_2$ through Langmuir-Hinshelwood
mechanism as inefficient for a grain size distribution, and the amount
of atomic H as critical for the chemistry of the upper layers (2-4\,$H$).
The parametric method of \cite{Bron+etal_2014}, that accounts for
temperature fluctuations for small grains, appears more appropriate to
produce enough H$_2$.
\item \revisionA{Comparisons between the models reveal differences of more than 1-2 orders of magnitudes}
for some abundant gas-phase species. Most of these differences arise
from above $1.5 H$, despite the larger UV penetration in multi-grain
models.
\item The location of the H$_2$O vertical
snowline is partly regulated by the UV penetration but also by the
amount of H$_2$. In the upper disk layers, the amount of H$_2$ also
regulates the formation of H$_2$CO, although this molecule is mostly
formed, like CH$_3$OH, on grain surfaces.
\item CH$_3$CN, like CH$_3$OH to a lesser extent, is preferentially
formed on warmer grains that allow for a larger diffusivity of
precursors on the grain surface.
\item At 100\,au, in the midplane, disks with colder dust produce more
H$_2$O and CH$_4$ while those with warmer dust  produce more CO$_2$ on
grain surfaces.
\item As soon as CO is locked onto large (hence colder) grains
($a >1\mu$m), it remains trapped until the end of the simulated evolution.
\item \revisionA{The spread of temperatures available with grain size distributions
simultaneously provides surface chemistry for cold and warm grains
that cannot be mimicked by a unique dust temperature. In particular,
it allows depletion of C and O to proceed efficiently (by conversion
to CO$_2$ or CH$_4$ and H$_2$O on colder grains), while the COMs
production is boosted on the warmer grains. }
\item 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.
\item \revisionA{Another direct consequence of having different dust temperatures
at the same location is a spatial spread-out of the snowlines.} \revisionC{High angular observations
would be required to test this prediction.}

\end{itemize}

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.

\begin{acknowledgements}
We thank E.Bron for fruitful discussions about the H$_2$ 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.
\end{acknowledgements}