# H2 ro-vibrational excitation in protoplanetary disks and its effects on the
chemistry

Maxime Ruaud NASA Ames Research Center, Moffett Blvd., Mountain View, CA
94035, USA  Carl Sagan Center, SETI Institute, Mountain View, CA 94035, USA

(Received XX; Revised YY; Accepted ZZ)

###### Abstract

The effect of H2 ro-vibrational excitation on the chemistry of protoplanetary
disks is studied using a framework that solves for the disk physical and
chemical structure and includes a detailed calculation of H2 level
populations. Chemistry with ro-vibrationally excited H2 is found to be
important for the formation of several commonly observed species in disks and
this work demonstrates the need to accurately treat PDR chemistry in disks if
we are to make inferences on the chemical state of the disk during planet
formation epochs. This is found to be even more critical for molecules like
C2H, CN or HCN that are commonly used to infer changes in the elemental disk
C/O and N/O ratios, with implications for planetesimal formation and the
composition of exoplanet atmospheres. Computed vertical column densities with
the full H2 population calculation are increased by $\sim 1-2$ orders of
magnitude for molecules such as CN, HCN/HNC compared to calculations with no
treatment of excited H2. For the commonly used pseudo-level approximation, the
computed columns of these molecules are overestimated by a factor of $\sim
3-5$ when compared to the full model. We further note that the computed
abundance for these molecules strongly depends on the strength of the FUV
photons at energies that pump H2 (i.e. 11-13.6 eV), which is not well
constrained in disks, and that rate constants as a function of H2 ro-
vibrational levels for the key reaction N + H${}_{2}\rightarrow$ NH are needed
for a more accurate assessment of CN/HCN chemistry but are currently
unavailable.

Protoplanetary disks (1300), Astrochemistry (75), Chemical abundances (224),
Star formation (1569)

††journal: ApJ

##  1 Introduction

Chemistry involving ro-vibrationally excited H2 has long known to be
significant for several key reactions in photon dominated regions (PDRs) of
the interstellar medium (e.g. Agúndez et al., 2010; Nagy et al., 2013). In
these regions, the absorption of FUV photons in Lyman and Werner bands
efficiently pumps H2 in its excited electronic states. Approximately 10% of
the absorbed photons lead to H2 dissociation, while 90% radiatively de-excite
to bound ro-vibrational levels of the ground electronic state. These excited
ro-vibrational levels decay by means of quadrupole transitions which result in
the emission of infrared photons (e.g. Black & Dalgarno, 1976). The energy in
these excited ro-vibrational levels is sufficient to overcome the reaction
energy barrier of several reactions. Reactions of excited H2 with C+ and S+
for instance are responsible for the formation of CH+ and SH+ ions in PDRs
(e.g. Agúndez et al., 2010; Nagy et al., 2013). The efficiency of these
reactions critically depends on the levels population of H2. For C+ for
instance, the reaction with H2 in ground vibrational state is endothermic by
$\sim 4500$ K but becomes barrierless when reacting with H2 in $v=1$ (Zanchet
et al., 2013b). In the case of S+, the reaction is endothermic by $\sim 10000$
K and H2 needs to be in $v\gtrsim 2$ to yield significant reactivity (Zanchet
et al., 2019).

Chemistry with ro-vibrationally excited H2 in disks is often approximated by
the so called pseudo-level approximation, in which all the excited vibrational
levels of H2 are grouped together into a single vibrational excited state
H${}_{2}^{*}$, corresponding to a pseudo-level $v^{*}=6$ and having a weighted
averaged energy of $\sim 30000$ K (see Tielens & Hollenbach, 1985, and
references therein). The use of this approximation has been shown to lead to
satisfactory results when accounting for processes such as H2 self-shielding
and photodissociation as compared to detailed H2 excitation calculations (see
Fig. 1 & 2 of Draine & Bertoldi, 1996, for instance). However, the use of this
approximation to solve for chemistry with excited H2 may lead to unrealistic
results as the rates critically depends on the fractional populations in each
ro-vibrational state. In particular, rate coefficients of species reacting
with H${}_{2}^{*}$ are often approximated using the rate for reaction with H2
by reducing the exponential factor by the corresponding energy of the level
(i.e. $\sim 30000$K, see Tielens & Hollenbach, 1985). Recently, Cazzoletti et
al. (2018) used this approximation to explain the bright CN emission around
the well studied disk TW Hya. They proposed that CN originates from the PDR
surface by assuming a barrierless reaction of N with H${}_{2}^{*}$, although
there is a barrier of 12600K for it to react with ground state H2. Visser et
al. (2018) used the same assumption to conclude that most of the observed HCN
emission also comes from the PDR layer. Kamp et al. (2017) follow a more
conservative approach where the rate for reaction with H${}_{2}^{*}$ is
obtained by subtracting the energy corresponding to the first vibrational
excited level $v=1$ for which $E_{vJ}\sim 6000$K. Full calculations accounting
for the quantum ro-vibrational state of H2 have been limited to few reactions
such as H2 with C+, He+, S+, O, OH and CN, with important implications for the
formation of molecules such as CH+ and SH+ (Agúndez et al., 2010, 2018). The
goal of this article is to explore the effects of the full vibrational level
calculations of H2 on the chemistry of protoplanetary disks. Modeling results
are presented in Section 2 and Section 3 contains our conclusions.

##  2 Modeling results

In what follows, we consider the fiducial disk model of Ruaud & Gorti (2019)
to study the effects of H2 excitation on the chemistry. The disk is assumed to
have a mass $M_{\text{disk}}=5\times 10^{-3}M_{\odot}$ and a dust-to-gas mass
ratio $\Sigma_{d}/\Sigma_{g}=0.01$, and is constant with radius but it varies
with height owing to dust settling (see Appendix A of Ruaud & Gorti, 2019, for
more details). Stellar parameters are assumed to be $M_{*}=1M_{\odot}$,
$R_{*}=2.45R_{\odot}$ and an X-ray luminosity of $10^{30}$ erg s-1. The
stellar spectrum is the TW Hya spectrum taken from the Leiden database (Heays
et al., 2017).

###  2.1 H2 ro-vibrational excitation and level populations

The model includes 100 ro-vibrational levels of H2 in the ground electronic
state. The level populations of H2 are computed by including the effects of
(1) radiative transitions, (2) collisional processes with He, H, H+ and H2,
(3) FUV radiative pumping (4) X-ray excitation (5) excitation at H2 formation
and (6) the ortho/para conversion on grains. More details on the
implementation of each of these process is given in the appendix. Even though
X-ray pumping is taken into account, we find that it does not play a
significant role in the non-thermal excitation of H2 when compared to FUV
pumping. This result is similar to earlier modeling results of Nomura et al.
(2007) where they study the relative roles of X-ray and UV excitation of H2.
Therefore, in what follows, we do not discuss the effects of X-ray pumping on
chemistry any further.

Figure 1 shows the model calculated vertical column densities of all H2 levels
included in the model, divided by statistical weight, at $r=1$, 10 and 100 au
from the star, and expected LTE values (dashed lines). For $r=1$ and 10 au,
the lowest rotational levels of the ground vibrational level (i.e. $v=0$ and
$J\leq 8$) are at LTE and higher energy levels are mainly populated by FUV
photon pumping. For $r=100$ au, only levels up to $J=3$ are thermally
populated and departure from LTE induced by FUV photon pumping occurs at lower
energy than for $r=1$ and 10 au. This is due to the the lower temperatures and
densities of the outer disk as compared to the inner disk (i.e. radial
temperature gradient). Rotational distributions in vibrationally excited
states are characterized by rotational temperatures typical of fluorescent
emission and ranging from $T_{\text{rot}}\sim$1000K and $\sim$2500K (e.g.
Draine & Bertoldi, 1996). Curvature observed in LTE values is explained by the
vertical temperature gradient. Regions closer to the disk midplane, which
dominate the contribution to the column of the lowest energy H2 levels, are
colder than regions located in the disk atmosphere. As a result, the LTE curve
tends to flatten as the energy of the level increases.

Figure 1:  Vertically integrated H2 column density $N(v,J)$ divided by the
level degeneracy $g(J)$ as a function of the energy of the level at $r=1$, 10
and 100 au from the star. Dashed lines are at LTE.

Figure 2 shows the computed abundance maps of H2 in different $(v,J)$ states.
As seen in this figure (as well as Fig. 1), most of the H2 in the disk is in
the ($v=0,J=0$) and ($v=0,J=1$) ground vibrational states and the abundance of
H2 in higher ro-vibrational states tends to decrease as the energy levels of
the states increase. The abundance of the H2 ($v=1,J=0$) state for instance
reaches a peak abundance of few $10^{-6}$ with respect to the density
n${}_{\text{H}}$, while higher energy states have abundances $\lesssim
10^{-7}$. As will be discussed below, chemistry with excited H2 critically
depends on the fractional populations in each of these ro-vibrational states.

Figure 2:  Computed abundances of H2 in different ro-vibrational levels as a
function of the disk radius $r$ and normalized height $z/r$. Solid line
contours show the $A_{V}$ to the star.

###  2.2 Chemistry involving ro-vibrationally excitated H2

The effects of the full ro-vibrational level calculations of H2 on the
chemistry is explored for species having reaction barriers with ground state
H2. In the chemical network used here (i.e. the KIDA network, Wakelam et al.,
2015), there are $\sim 80$ reactions, including important reactions with
species such as C+, S+, N+, C, O, N, and OH that mediate pathways to many
other chemical species. Rate constants have been measured and/or calculated
for some reactions with species such as C+, S+ and O. For
$\text{C}^{+}(^{2}P)+\text{H}_{2}\rightarrow\text{CH}^{+}+\text{H}$, Gerlich
et al. (1987) experimentally derived rate constants for the reaction with
H${}_{2}(v=0,J=0-7)$ and for which we use analytical fits given in Agúndez et
al. (2010). For higher energy levels we use results from quantum calculations
carried out by Zanchet et al. (2013b) and for which fits are available for
H${}_{2}(v=1,J=0)$ and H${}_{2}(v=2,J=0)$ (see Table 3 of Zanchet et al.,
2013a). For $\text{S}(^{4}S)+\text{H}_{2}\rightarrow\text{SH}^{+}+\text{H}$,
we use rate constants given in Zanchet et al. (2013a) for reaction with
H${}_{2}(v=0,J=0)$ and $(v=1,J=0)$ states, and those given in Zanchet et al.
(2019) for reaction with H${}_{2}(v=2-5,J=0)$ states. In the case of
$\text{O}(^{3}P)+\text{H}_{2}\rightarrow\text{OH}+\text{H}$, we use fits to
quantum calculations performed by (Sultanov & Balakrishnan, 2005) and for
which expressions are given in Agúndez et al. (2010) for reaction with
$\text{H}_{2}(v=0-3,J=0)$. For these three reactions, since rate constants are
available for $J=0$ only (i.e. except for $\text{C}^{+}(^{2}P)+\text{H}_{2}$
for which we have rates with H${}_{2}(v=0,J=0-7)$), rate constants of higher
$J$’s are assumed to be equal to the rate of the corresponding vibrational
level in $J=0$. Moreover, for these reactions, rates for higher energy levels
than the last available state are assumed to be equal to the rate in this
state. For the remaining reactions having a barrier for reaction with H2,
reaction rates are approximated by

$k_{vJ}=\alpha(T/300)^{\beta}\exp(-\max[\gamma-
E_{vJ},0]/T)~{}[\text{cm}^{3}\text{s}^{-1}]$ (1)

where $\alpha$, $\beta$ and $\gamma$ are the usual rate constant parameters.
This expression assumes that the reaction proceeds without a barrier when the
energy of the level overcomes the energy barrier of the reaction. Note that
this assumption could be prone to error as rate constant enhancements are
usually specific to each system and difficult to predict (see also Agúndez et
al., 2010). One such example is $\text{O}(^{3}P)+\text{H}_{2}$, which has an
activation barrier of $\sim 4000$K, but for which an activation barrier
persists even when H2 is in the $v=3$ state, for which $E_{v}\sim 17000$K
(Sultanov & Balakrishnan, 2005; Agúndez et al., 2010).

Three types of models are considered in the following: (1) a reference model
in which the effect of the chemistry with ro-vibrationally excited H2 is
neglected111It should however be noted that reaction parameters that are
determined experimentally have a temperature dependence that implicitly
consider the excitation of H2., (2) a model in which all the excited
vibrational levels of H2 are grouped together into a single vibrational
excited state having a weighted averaged energy of 30160 K corresponding to
H${}_{2}(v=6)$ (the so-called pseudo-level approximation), and (3) a model
where the full H${}_{2}(v,J)$ population calculation is taken into account
(see Table 1). For model M2, we follow the treatment described in Tielens &
Hollenbach (1985) in which H${}_{2}^{*}$ is formed by FUV pumping and
destroyed by radiative decay and collisional de-excitation by H and H2. The
FUV pumping rate is assumed to be 9 times the H2 photodissociation rate and
the radiative decay rate is $2\times 10^{-7}$ s-1. For collisional de-
excitation by H and H2 we use rates given in Equation (A14) of Tielens &
Hollenbach (1985). These rates have recently been updated by Visser et al.
(2018) based on more recent calculations. A comparison with the rates used in
Visser et al. (2018) however shows that the impact of these updated rates is
relatively small. This is due to the fact that radiative decay dominates
H${}_{2}^{*}$ destruction in most of the disk. The only location where the use
of these updated rates leads to differences is the inner disk (i.e. $r\lesssim
20$ au) where higher densities result in collisional de-excitation dominating
radiative decay. Finally, reaction rates with H${}_{2}^{*}$ in this model are
computed using Eq. 1 and assuming $E_{v^{*}}=30160$K. A comparison between the
computed abundance map of H${}_{2}^{*}$ calculated with the pseudo-level
approximation and for all levels with $E_{vJ}\geq 30000$K of the full model is
shown in Fig. 3.

Figure 3:  Computed abundance of excited H2 as a function of the disk radius $r$ and normalized height $z/r$. The top figure shows result obtained in the case of the pseudo-level approximation and the bottom figure shows results obtained with the full model and for levels with $E_{vJ}\geq 30000$K. Table 1:  Model description. | Model description
---|---
M1 | Reference model in which the chemistry
| of excited H2 is neglected
M2 | Chemistry of excited H2 assuming
| the pseudo-level approximation
M3 | Chemistry of excited H2 with
| full H2 population calculation

In the following, the chemistry in regions with $z/r<0.2$, which trace regions
where grain surface chemistry is important, is not presented. Readers
interested in the chemistry of these regions can refer to papers such as
Semenov & Wiebe (2011), Walsh et al. (2014) and Ruaud & Gorti (2019), for
detailed models of gas-grain chemistry using large surface reaction networks.
Figs. 4 and 5 show the computed abundance maps of CH+, SH+, OH, H2O, C2H, CN,
HCN, HNC and N2H+ computed in each of these models and show that the chemistry
with ro-vibrationally excited H2 is important in the surface layers of the
disk (i.e. in regions with $z/r>0.2$, where H2 ro-vibrational levels are
efficiently populated by the absorption of FUV photons, see Fig. 2). The use
of the pseudo-level approximation as compared to the full H2 population shows
important differences for species such as CN, HCN and HNC for which the
computed abundances in the PDR layer are overestimated by approximately one
order of magnitude in the pseudo-level approximation model (M2) as compared to
the full H2 population model (M3). These differences are also seen in Fig. 6
which shows vertical column density ratios for models M1 and M2 with respect
to the full calculation M3 of some abundant gas-phase species. Details on some
important reaction pathways are discussed below.

Figure 4:  Computed abundance map of CH+, SH+ and OH as a function of the disk
radius $r$ and normalized height $z/r$ for model M1, M2 and M3 (see Table 1).

Figure 5:  Same as Fig. 4 but for H2O, C2H, CN, HCN, HNC and N2H+. Figure 6:
Vertical column density ratios of a selection of gas-phase molecules as a
function of the distance from the star. Solid lines correspond to ratios
between column densities computed in model M1 (no chemistry with excited H2)
and M3 (full model) while dashed lines show ratios between results obtained
from model M2 (pseudo-level approximation) and M3. The gray shaded area
indicate a factor of 5 deviation compared to results obtained from the full
model.

####  2.2.1 Carbon chemistry

The inclusion of C+ \+ H${}_{2}(v,J)$ has an important impact on the formation
of CH+ as compared to the reference model M1. At $r>20$ au, CH+ abundance is
increased by almost two orders of magnitude in the M2 and M3 models as
compared to M1, while the M3 model has $\sim 3$ times higher CH+ than the M2
model. In the full calculation model, CH+ is mainly formed by reaction with H2
in the $(v=0,J=8)$, $(v=1,J=0)$ and $(v=1,J=1)$ state, for which
$E_{vJ}\gtrsim 4640$ K. Enhanced formation of CH+ in the irradiated layer is
in agreement with Agúndez et al. (2018), even though our computed abundances
are lower by a factor of $\sim 10$. This could be the result of different disk
physical parameters and/or the use of a different stellar spectrum. As seen in
Figs. 5 and 6, the enhanced formation of CH+ only has a moderate impact on the
formation of C2H, a bright molecule at millimeter wavelengths in disks for
which most models fail to reproduce observed emission without invoking
enhanced C/O ratios (e.g. Bergin et al., 2016; Cleeves et al., 2018; Miotello
et al., 2019). In all three models, C2H formation starts with CH + C+ which
forms C${}_{2}^{+}$. C${}_{2}^{+}$ then reacts with H2 to form successively
C2H+ and C2H${}_{2}^{+}$, which then recombine to form C2H. The key molecule
in this sequence is CH whose formation is primarily initiated by the radiative
association of C+ and H2 to form CH${}_{2}^{+}$. CH${}_{2}^{+}$ then reacts
with H2 to form CH${}_{3}^{+}$, which then dissociatively recombine to form
CH. In model M2 and M3, the reaction of C with excited H2 also contributes to
the formation of CH. In model M3, enhanced formation of CH+ plays a role on
the formation of C2H at r$\lesssim 40$ au, where C2H abundance is enhanced by
a factor of $\sim 5$ when compared to model M1. In this region, the enhanced
formation of CH+ contributes to the formation of CH through CH+ \+ H2 which
forms CH${}_{2}^{+}$ and then CH through the sequence described above. In
model M2, the decreased abundance of C2H at $30\lesssim r\lesssim 100$ au is
found to be linked to the enhanced production of CN in this region which
efficiently competes with C2H formation. It should however be noted that the
enhanced formation of C2H due to chemistry with excited H2 may not be enough
to explain observed C2H emission using ISM-like C/O ratios and C/O$\gtrsim 1$
may still be required. Increased vertical column densities of at least an
order of magnitude seems to be needed in order to reproduce C2H observations
as compared to what is predicted with standards C/O ratios (e.g. Bergin et
al., 2016; Cleeves et al., 2018; Miotello et al., 2019). It is however still
unclear whether this change in disks C/O ratios is induced by the differential
sequestration of oxygen and carbon at the surface of large grains in the disk
miplane (e.g. Krijt et al., 2018, 2020) or by photodestruction of carbon
grains in disks atmosphere (e.g. Anderson et al., 2017; Bosman et al., 2021).

####  2.2.2 Oxygen chemistry

The inclusion of O + H${}_{2}(v,J)$ has only a moderate impact on the computed
abundances of OH and water in the M3 model as compared to M1 (see Fig. 4 and
5). In both models, OH and water formation are mainly driven by X-rays through
the production of H+ and H${}_{3}^{+}$, and their successive reaction with
atomic oxygen. As seen in these figures, M2 overestimates slightly the
abundance of OH and water (i.e. by a factor of $\sim 5$) in a thin layer
located at $z/r\sim 0.3-0.4$ and extending from $r=20$ and 100 au as compared
to M3. The higher abundance of OH and H2O in this layer are found to be linked
to the enhanced production of CN in this model and in particular a slightly
lower abundance of C+ which is one of the main destruction channels of OH and
water in this layer. The chemistry of other molecules such as NO, HNO and NO+
are also found to be affected by chemistry with ro-vibrationally excited H2
(see Fig. 6). For these molecules, this is due to the enhanced formation of NH
and NH2 (see Section 2.2.4) and their subsequent reaction with O.

####  2.2.3 Sulfur chemistry

Even though the abundance of SH+ is low throughout the disk, the inclusion of
S+ \+ H${}_{2}(v,J)$ has an important impact on its formation (see Fig. 4). As
seen in this figure, models M2 and M3 give relatively different results with
SH+ being overestimated in the vertical directions of model M2 and slightly
underestimated in the inner part of the disk (i.e. at $r\lesssim 100$ au).
Quantum calculations show that H2 needs to be in $v\geq 2$ states for the
reaction to proceed without a barrier (Zanchet et al., 2019) and its
efficiency thus depends on the populations in the H${}_{2}(v\geq 2)$ states.
In model M3, H${}_{2}(v\geq 2)$ states are more efficiently populated in the
inner disk (the abundance of H${}_{2}(v\geq 2)$ is on the order of few
$10^{-8}-10^{-7}$ up to $r\sim 100$ au, see Fig. 2). At greater distance from
the star, the populations in these states decrease significantly hence
lowering the formation efficiency of SH+. Except for SH+, sulfur chemistry is
not particularly affected by chemistry with excited H2.

####  2.2.4 Nitrogen chemistry

Nitrogen chemistry is strongly affected by the chemistry with ro-vibrationally
excited H2. In particular, the abundance of species such as CN, HCN and HNC is
increased by a factor of $\sim 100$ in model M3 as compared to the reference
model M1 (see Fig. 5). For model M2, the abundance of these species is
increased by almost a factor of $\sim 10$ as compared to model M3 at the peak
abundance (and a factor of 3 to 5 in terms of the computed vertical column
densities, see Fig. 6). In both model M2 and M3, the formation of CN and
HCN/HNC starts from N + H${}_{2}(v,J)$ to form NH for which the reaction
barrier is 12650 K. NH then reacts with C+ to form CN+. CN is then formed by
charge exchange with H while HCN/HNC are formed through reaction of CN+ with
H2 to form successively HCN+ and HCNH+, which then recombines with electrons
to form HCN and HNC. The overestimation of CN and HCN/HNC in model M2 as
compared to model M3 comes from the fact that the fractional population of H2
in $(v>2)$ states, for which $E(v,J)>11600$ K, are relatively low as compared
to what is predicted for H${}_{2}^{*}$ (see Fig. 2 and 3). As seen in Fig. 5
and 6, the increased formation of CN and HCN also impacts N2H+ for which the
computed vertical column densities are increased by a factor of at least 3 in
models M2 and M3 as compared to model M1. The increased production of NH is
also found to impact the formation of molecules such as NH2, NO, HNO, NO+,
CN+, HCNH+, CNC+, NS+, SiN and SiN+ (see Fig. 6). Recently, Cazzoletti et al.
(2018) used the pseudo-level approximation to explain CN emission around the
well studied disk TW Hya as originating from the PDR surface and being formed
by the same reaction of N with excited H2, where they assumed no barrier (i.e.
similar to our model M2). Visser et al. (2018) used the same assumption to
conclude that most of the observed HCN emission also comes from the PDR layer.
While the impact of this reaction on the formation of CN and HCN is confirmed
by models shown here, it is however found that the use of the pseudo-level
approximation lead to a factor of 3 to 5 overestimation in their computed
vertical column densities (see Fig. 6).

###  2.3 Sensitivity to the stellar FUV photons flux

In this section we study the effect of the FUV flux on the modeling results.
As seen previously, the chemistry of excited H2 sensitively depends on the
fractional populations in each ro-vibrational state of H2 which in turn depend
on FUV pumping. We focus on the 11 to 13.6eV photons range which is poorly
constrained in disks (e.g. France et al., 2014) but yet fundamental for H2
excitation and disk photochemistry. In particular, it is important to note
that most H2 fluorescent lines only depend on the photon flux in this energy
range. This is because FUV pumping tends to be dominated by electronic
transitions from the lower vibrational levels of the ground state
$X^{1}\Sigma_{g}^{+}$, (e.g. electronic Lyman and Werner transitions from
H${}_{2}(v=0,J<8)$ all lie at $\lambda<1150$Å). In the following, we consider
two variations of model M3 in which we vary the FUV flux in the 11 to 13.6 eV
range by a factor of 10. It should however be noted that the disk physical
structure is kept the same and that only the chemistry is recomputed in these
models. Fig. 7 compares the effect of varying the FUV flux on the computed
vertical column density of a selection of gas-phase molecules which are
strongly affected by a change in the FUV flux. In particular, the computed
vertical column densities of CH+, SH+, CN and HCN vary by a factor of $\sim
10$ when the FUV flux in the 11-13.6 eV range varies by a factor of 10. C2H is
found to be relatively insensitive to the FUV flux, except in a small region
located at $30\lesssim r\lesssim 100$ au where its computed vertical column
density varies by at most a factor of $\sim 5$ when compared to the reference
model (i.e. $G_{11-13.6}=1$). Increased CN, HCN and C2H column densities with
FUV flux could explain some of the positive correlations found in the emission
of these molecules (e.g. Bergner et al., 2019).

Figure 7:  Computed vertical column densities of CH+, SH+, N2H+, C2H, CN and
HCN as a function of the disk radius $r$ for model M3 and where the FUV flux
in the 11 to 13.6eV photons range was varied by a factor of 10.
$G_{11-13.6}=1$ used for the reference model corresponds to a photon flux of
$4\times 10^{13}$ photons cm-2 s-1 in the 11-13.6 eV range and at 1 au from
the star.

##  3 Conclusion

The effect of H2 ro-vibrational excitation on the chemistry of protoplanetary
disks is studied using a framework that solves for the disk physical and
chemical structure and includes a detailed calculation of H2 level
populations. Chemistry with ro-vibrationally excited H2 is found to be
efficient in the irradiated surface layer of disks where H2 is efficiently
pumped by FUV photons and its efficiency sensitively depends on the fractional
populations in each ro-vibrational state of H2. Chemistry with excited H2 is
found to play an important role on the formation of several commonly observed
species in disks and in particular CN, HCN and HNC, for which the computed
vertical column densities are increased by 1 to 2 orders of magnitude for the
full H2 population calculation compared to calculations with no treatment of
ro-vibrationally excited H2, and by factors of $\sim 3-5$ when compared to the
commonly used pseudo-level approximation where all excited H2 is considered to
be at a single energy. Chemistry with ro-vibrationally excited H2 also
enhances the formation of molecules such as SH+, CH+, N2H+, NH, NH2, NO, HNO,
NO+, CN+, HCNH+, CNC+, NS+, SiN and SiN+. For most of these molecules, except
CH+ and SH+, the key reaction is N + H${}_{2}\rightarrow$ NH which has a
barrier of 12650 K with ground state H2. Even though few studies exist on this
particular system (e.g Pascual et al., 2002), rate constants as a function of
H2 ro-vibrational levels are missing. Here it was assumed that the reaction
proceeds without a barrier when the energy of the level overcomes the energy
barrier of the reaction (see Eq. 1). However this could be a strong over
simplification of the problem and rate constants should therefore be obtained
for a more accurate assessment of CN/HCN chemistry. This is even more critical
that molecules such as C2H, CN or HCN are now commonly used to infer changes
in the elemental disk C/O and N/O ratios (e.g. Bergin et al., 2016; Cleeves et
al., 2018; Miotello et al., 2019), with implications for planetesimal
formation and the composition of exoplanet atmospheres. It is also shown that
the formation of these molecules is greatly affected by a change in the FUV
flux in the 11 to 13.6eV energy range. In particular we find a positive
correlation between the increase in the FUV flux and the computed vertical
column densities of molecules such as CN, HCN and to a lesser extent C2H.
However, the strength of the FUV photon flux at energies that pump H2 in
protoplanetary disks (i.e. mostly 11-13.6 eV) is not well constrained and this
further impacts our ability to make inferences on the chemical state of the
disk during planet formation epochs. With the upcoming launch of JWST, it will
soon be possible to probe excited H2 which will help to further constrain
disks chemistry and potential molecular tracers of planet formation processes.
Finally, it is found that computed vertical column densities from the common
pseudo-level approximation for H2 chemistry for most species agree within a
factor $\lesssim 5$ with the more detailed calculations presented here.

I thank O. Roncero and A. Zanchet for sharing their computational results on
state-to-state chemistry and for fruitful discussions around these data. I
acknowledge support by the National Aeronautics and Space Administration
through the NASA Astrobiology Institute under Cooperative Agreement Notice
NNH13ZDA017C issued through the Science Mission Directorate.

##  Appendix A Radiative transitions

Due to its homonuclearity and symmetry, H2 has no permanent dipole moment. As
a result, electric dipole ro-vibrational transitions of H2 are forbidden and
only weak quadrupole transitions are allowed. For ro-vibrational levels within
the ground electronic state we use energy levels compiled in Dabrowski (1984)
and quadrupole radiative transition probabilities from Wolniewicz et al.
(1998). For radiative transitions of the Lyman and Werner band systems (i.e.
$B^{1}\Sigma_{u}\rightarrow X^{1}\Sigma^{+}_{g}$ and $C^{1}\Pi_{u}\rightarrow
X^{1}\Sigma^{+}_{g}$), we use energy levels and transition probabilities
calculated by Abgrall et al. (1993a, b, 2000).

##  Appendix B Collisional processes

Collision excitation are included for collisions with H (Wrathmall & Flower,
2007), He (Flower et al., 1998), H2 (Flower & Roueff, 1998) and H+ (Gerlich,
1990). For collisions with He and H2 we use fit to collisional de-excitation
rates provided by Le Bourlot et al. (1999) in the range $100\lesssim
T_{\mathrm{gas}}\lesssim 6000$ K. Collisional excitation rates are then
computed using the detail balance equation

$q(v^{\prime}J^{\prime}\rightarrow
vJ)=\frac{g_{J}}{g_{J}^{\prime}}q(vJ\rightarrow
v^{\prime}J^{\prime})\exp\Bigg{(}-\frac{E_{vJ}-E_{v^{\prime}J^{\prime}}}{T}\Bigg{)}$
(B1)

where $g_{J}$ is the statistical weight of the level $J$ and $E_{vJ}$ its
energy. Statistical weights are given by $g_{J}=(2J+1)$ for para-H2 and
$g_{J}=3(2J+1)$ for ortho-H2.

##  Appendix C FUV radiative pumping

The FUV pumping rate from a ro-vibrational level $(v,J)$ of the ground
electronic sate to level $(v^{*},J^{*})$ of an upper electronic state is given
by (e.g van Dishoeck & Black, 1986; Sternberg & Dalgarno, 1989)

$P(vJ\rightarrow v^{*}J^{*})=\int I_{\nu}\sigma_{\nu}(vJ\rightarrow
v^{*}J^{*})\exp\big{(}-N_{vJ}\sigma_{\nu}(vJ\rightarrow
v^{*}J^{*})\big{)}d\nu$ (C1)

where $\nu$ is the transition frequency, $I_{\nu}$ is the FUV photon intensity
(in photons.cm-2.s-1.Hz-1), $\sigma_{\nu}$ is the photon absorption cross
section of the transition, and $N_{vJ}$ is the column density of H2 molecule
in level $(v,J)$. Assuming that the FUV intensity is constant over the line
profile and defining the equivalent width of the absorption line

$W(vJ\rightarrow
v^{*}J^{*})=\int\Big{[}1-\exp\big{(}-N_{vJ}\sigma_{\nu}(vJ\rightarrow
v^{*}J^{*})\big{)}\Big{]}d\nu$ (C2)

as well as the line self-shielding function

$\theta(vJ\rightarrow v^{*}J^{*})=\Bigg{(}\frac{\pi
e^{2}}{m_{e}c}f_{\text{osc}}(vJ\rightarrow
v^{*}J^{*})\Bigg{)}^{-1}\frac{dW(vJ\rightarrow v^{*}J^{*})}{dN_{vJ}}$ (C3)

equation C1 may be rewritten

$P(vJ\rightarrow v^{*}J^{*})=\frac{\pi
e^{2}}{m_{e}c}I_{\nu}f_{\text{osc}}(vJ\rightarrow
v^{*}J^{*})\theta(vJ\rightarrow v^{*}J^{*})$ (C4)

where $f_{\text{osc}}(vJ\rightarrow v^{*}J^{*})$ is the oscillator strength of
the transition. The line self-shielding function is calculated following the
analytical treatment of Federman et al. (1979).

##  Appendix D Excitation by X-rays

Fast electrons created by X-rays ionization can excite H2 ro-vibrational
levels via collisions. This can occur either by direct collisional excitation
within the ground vibrational level or by cascade following excitation of the
electronic states (e.g. Gredel & Dalgarno, 1995; Tiné et al., 1997). Tiné et
al. (1997) provide entry probabilities for excitation rates by X-rays for H2
in the ground vibrational state to H${}_{2}(v=0-14,J=0-11)$ and use tabulated
entry probabilities obtained for a fractional ionization of $10^{-4}$.

##  Appendix E Excitation at formation and o/p conversion on grains

For the excitation at formation, we follow an approach similar to Le Petit et
al. (2006) and consider an equipartition of the energy released at formation
(i.e. 4.5 eV) with 1/3 of the energy being used as internal excitation of H2,
1/3 being transferred to the grain and 1/3 being converted into kinetic
energy. The fraction of H2 formed in each ro-vibrational levels follows

$f_{vJ}(T)=\frac{g_{J}\exp(-E_{vJ}/T)}{\sum g_{J}\exp(-E_{vJ}/T)}$ (E1)

where $T=8734$K. As noted in Le Petit et al. (2006), this is different than
assuming a Boltzmann distribution at $T=17322$ K (i.e. $\sim 1.5$eV), which
would be equivalent to assume that approximately half (i.e. $\sim 2.5$ eV) of
the energy released at formation is used as internal excitation of H2. For the
ortho/para conversion on grains, we follow the approach of Le Bourlot (2000)
and consider that conversion occurs with an efficiency set by

$\eta_{c}=\exp\big{(}-\tau_{c}k_{\text{ev}})$ (E2)

where $\tau_{c}$ is the conversion timescale and $k_{\text{ev}}$ the
evaporation rate of H2 from the grains. The o/p conversion rate is then
calculated by

$k_{\text{o/p}}=\alpha_{\text{H}_{2}}v_{\text{H}_{2}}\langle\sigma_{d}n_{d}\rangle\eta_{c}$
(E3)

where $\alpha_{\text{H}_{2}}$ is the H2 sticking probability taken from Matar
et al. (2010) and Chaabouni et al. (2012), $v_{\text{H}_{2}}$ is the thermal
speed of H2 molecules and $\langle\sigma_{d}n_{d}\rangle$ the product of the
dust cross sectional area and the dust density averaged over the grain size
distribution.

## References

  * Abgrall et al. (2000) Abgrall, H., Roueff, E., & Drira, I. 2000, A&AS, 141, 297, doi: 10.1051/aas:2000121
  * Abgrall et al. (1993a) Abgrall, H., Roueff, E., Launay, F., Roncin, J. Y., & Subtil, J. L. 1993a, A&AS, 101, 323
  * Abgrall et al. (1993b) —. 1993b, A&AS, 101, 273
  * Agúndez et al. (2010) Agúndez, M., Goicoechea, J. R., Cernicharo, J., Faure, A., & Roueff, E. 2010, ApJ, 713, 662, doi: 10.1088/0004-637X/713/1/662
  * Agúndez et al. (2018) Agúndez, M., Roueff, E., Le Petit, F., & Le Bourlot, J. 2018, A&A, 616, A19, doi: 10.1051/0004-6361/201732518
  * Anderson et al. (2017) Anderson, D. E., Bergin, E. A., Blake, G. A., et al. 2017, ApJ, 845, 13, doi: 10.3847/1538-4357/aa7da1
  * Bergin et al. (2016) Bergin, E. A., Du, F., Cleeves, L. I., et al. 2016, ApJ, 831, 101, doi: 10.3847/0004-637X/831/1/101
  * Bergner et al. (2019) Bergner, J. B., Öberg, K. I., Bergin, E. A., et al. 2019, ApJ, 876, 25, doi: 10.3847/1538-4357/ab141e
  * Black & Dalgarno (1976) Black, J. H., & Dalgarno, A. 1976, ApJ, 203, 132, doi: 10.1086/154055
  * Bosman et al. (2021) Bosman, A. D., Alarcon, F., Zhang, K., & Bergin, E. A. 2021, arXiv e-prints, arXiv:2101.12502.  https://arxiv.org/abs/2101.12502
  * Cazzoletti et al. (2018) Cazzoletti, P., van Dishoeck, E. F., Visser, R., Facchini, S., & Bruderer, S. 2018, A&A, 609, A93, doi: 10.1051/0004-6361/201731457
  * Chaabouni et al. (2012) Chaabouni, H., Bergeron, H., Baouche, S., et al. 2012, A&A, 538, A128, doi: 10.1051/0004-6361/201117409
  * Cleeves et al. (2018) Cleeves, L. I., Öberg, K. I., Wilner, D. J., et al. 2018, ApJ, 865, 155, doi: 10.3847/1538-4357/aade96
  * Dabrowski (1984) Dabrowski, I. 1984, Canadian Journal of Physics, 62, 1639, doi: 10.1139/p84-210
  * Draine & Bertoldi (1996) Draine, B. T., & Bertoldi, F. 1996, ApJ, 468, 269, doi: 10.1086/177689
  * Federman et al. (1979) Federman, S. R., Glassgold, A. E., & Kwan, J. 1979, ApJ, 227, 466, doi: 10.1086/156753
  * Flower & Roueff (1998) Flower, D. R., & Roueff, E. 1998, Journal of Physics B Atomic Molecular Physics, 31, 2935, doi: 10.1088/0953-4075/31/13/012
  * Flower et al. (1998) Flower, D. R., Roueff, E., & Zeippen, C. J. 1998, Journal of Physics B Atomic Molecular Physics, 31, 1105, doi: 10.1088/0953-4075/31/5/017
  * France et al. (2014) France, K., Schindhelm, E., Bergin, E. A., Roueff, E., & Abgrall, H. 2014, ApJ, 784, 127, doi: 10.1088/0004-637X/784/2/127
  * Gerlich (1990) Gerlich, D. 1990, J. Chem. Phys., 92, 2377, doi: 10.1063/1.457980
  * Gerlich et al. (1987) Gerlich, D., Disch, R., & Scherbarth, S. 1987, J. Chem. Phys., 87, 350, doi: 10.1063/1.453580
  * Gredel & Dalgarno (1995) Gredel, R., & Dalgarno, A. 1995, ApJ, 446, 852, doi: 10.1086/175843
  * Heays et al. (2017) Heays, A. N., Bosman, A. D., & van Dishoeck, E. F. 2017, A&A, 602, A105, doi: 10.1051/0004-6361/201628742
  * Kamp et al. (2017) Kamp, I., Thi, W. F., Woitke, P., et al. 2017, A&A, 607, A41, doi: 10.1051/0004-6361/201730388
  * Krijt et al. (2020) Krijt, S., Bosman, A. D., Zhang, K., et al. 2020, ApJ, 899, 134, doi: 10.3847/1538-4357/aba75d
  * Krijt et al. (2018) Krijt, S., Schwarz, K. R., Bergin, E. A., & Ciesla, F. J. 2018, ApJ, 864, 78, doi: 10.3847/1538-4357/aad69b
  * Le Bourlot (2000) Le Bourlot, J. 2000, A&A, 360, 656
  * Le Bourlot et al. (1999) Le Bourlot, J., Pineau des Forêts, G., & Flower, D. R. 1999, MNRAS, 305, 802, doi: 10.1046/j.1365-8711.1999.02497.x
  * Le Petit et al. (2006) Le Petit, F., Nehmé, C., Le Bourlot, J., & Roueff, E. 2006, ApJS, 164, 506, doi: 10.1086/503252
  * Matar et al. (2010) Matar, E., Bergeron, H., Dulieu, F., et al. 2010, J. Chem. Phys., 133, 104507, doi: 10.1063/1.3484867
  * Miotello et al. (2019) Miotello, A., Facchini, S., van Dishoeck, E. F., et al. 2019, arXiv e-prints, arXiv:1909.04477.  https://arxiv.org/abs/1909.04477
  * Nagy et al. (2013) Nagy, Z., Van der Tak, F. F. S., Ossenkopf, V., et al. 2013, A&A, 550, A96, doi: 10.1051/0004-6361/201220519
  * Nomura et al. (2007) Nomura, H., Aikawa, Y., Tsujimoto, M., Nakagawa, Y., & Millar, T. J. 2007, ApJ, 661, 334, doi: 10.1086/513419
  * Pascual et al. (2002) Pascual, R. Z., Schatz, G. C., Lendvay, G., & Troya, D. 2002, Journal of Physical Chemistry A, 106, 4125, doi: 10.1021/jp0133079
  * Ruaud & Gorti (2019) Ruaud, M., & Gorti, U. 2019, ApJ, 885, 146, doi: 10.3847/1538-4357/ab4996
  * Semenov & Wiebe (2011) Semenov, D., & Wiebe, D. 2011, ApJS, 196, 25, doi: 10.1088/0067-0049/196/2/25
  * Sternberg & Dalgarno (1989) Sternberg, A., & Dalgarno, A. 1989, ApJ, 338, 197, doi: 10.1086/167193
  * Sultanov & Balakrishnan (2005) Sultanov, R. A., & Balakrishnan, N. 2005, ApJ, 629, 305, doi: 10.1086/431356
  * Tielens & Hollenbach (1985) Tielens, A. G. G. M., & Hollenbach, D. 1985, ApJ, 291, 722, doi: 10.1086/163111
  * Tiné et al. (1997) Tiné, S., Lepp, S., Gredel, R., & Dalgarno, A. 1997, ApJ, 481, 282, doi: 10.1086/304048
  * van Dishoeck & Black (1986) van Dishoeck, E. F., & Black, J. H. 1986, ApJS, 62, 109, doi: 10.1086/191135
  * Visser et al. (2018) Visser, R., Bruderer, S., Cazzoletti, P., et al. 2018, A&A, 615, A75, doi: 10.1051/0004-6361/201731898
  * Wakelam et al. (2015) Wakelam, V., Loison, J. C., Herbst, E., et al. 2015, ApJS, 217, 20, doi: 10.1088/0067-0049/217/2/20
  * Walsh et al. (2014) Walsh, C., Millar, T. J., Nomura, H., et al. 2014, A&A, 563, A33, doi: 10.1051/0004-6361/201322446
  * Wolniewicz et al. (1998) Wolniewicz, L., Simbotin, I., & Dalgarno, A. 1998, ApJS, 115, 293, doi: 10.1086/313091
  * Wrathmall & Flower (2007) Wrathmall, S. A., & Flower, D. R. 2007, Journal of Physics B Atomic Molecular Physics, 40, 3221, doi: 10.1088/0953-4075/40/16/003
  * Zanchet et al. (2013a) Zanchet, A., Agúndez, M., Herrero, V. J., Aguado, A., & Roncero, O. 2013a, AJ, 146, 125, doi: 10.1088/0004-6256/146/5/125
  * Zanchet et al. (2013b) Zanchet, A., Godard, B., Bulut, N., et al. 2013b, ApJ, 766, 80, doi: 10.1088/0004-637X/766/2/80
  * Zanchet et al. (2019) Zanchet, A., Lique, F., Roncero, O., Goicoechea, J. R., & Bulut, N. 2019, A&A, 626, A103, doi: 10.1051/0004-6361/201935471