# Modeling Nitrogen Fractionation in the Protoplanetary Disk around TW Hya:
Model Constraints on Grain Population and Carbon-to-Oxygen Elemental Abundance
Ratio

Seokho Lee  National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka,
Tokyo 181-8588, Japan<EMAIL_ADDRESS>Hideko Nomura  National
Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
Kenji Furuya  National Astronomical Observatory of Japan, 2-21-1 Osawa,
Mitaka, Tokyo 181-8588, Japan  Jeong-Eun Lee  School of Space Research, Kyung
Hee University, 1732 Deogyeong-daero, Giheung-gu, Yongin-si, Gyeonggi-do,
Korea

###### Abstract

Observations conducted using the Atacama Large Millimeter/submillimeter Array
on the protoplanetary disk around TW Hya show the nitrogen fractionation of
HCN molecules in HC14N/HC15N $\sim$120 at a radius of $\sim$20 AU. In this
study, we investigated the physical and chemical conditions that control this
nitrogen fractionation process. To this end, a new disk model was developed,
in which the isotope-selective photodissociation of N2 and isotope-exchange
chemical reactions have been incorporated. Our model can successfully
reproduce the observed HCN column density when the elemental abundances of the
gas-phase carbon and oxygen are depleted by two orders of magnitude relative
to those in the interstellar medium and carbon is more abundant than oxygen
([C/O]elem$>$ 1). The isotope-selective photodissociation of N2 is the
dominant nitrogen fractionation process in our models. The observed
HC14N/HC15N ratio, which increases outwards, can also be reproduced by the
model by assuming that the small dust grains in the atmosphere of the outer
disk are depleted more than those in the inner disk. This is consistent with
grain evolution models, according to which small dust grains are continuously
replenished in the inner disk due to fragmentation of the large dust grains
that radially drift from the outer disk.

protoplanetary disk, astrochemistry, numerical

##  1 Introduction

The isotopic ratio of molecules is a powerful tool for investigating the
origin of solar system materials and for revealing the possible chemical link
between the solar system and the interstellar medium (ISM) (see e.g.,
Ceccarelli et al., 2014, for a recent review). Nitrogen has two stable
isotopes, 14N and 15N. The nitrogen isotope ratio, 14N/15N, is nonuniform in
the primitive materials of the solar system. Cometary ices, such as HCN and
NH3, have a 14N/15N ratio of $\sim$140 (e.g., Mumma & Charnley, 2011; Shinnaka
et al., 2016), which is a factor of three lower than the elemental nitrogen
isotope ratio, [14N/15N]elem, of 440 observed in solar winds (Marty et al.,
2011); that is, cometary ices are enriched in 15N. By contrast, bulk
meteorites show an isotope ratio between the two aforenoted values ($\sim$270,
Füri & Marty, 2015). Thus, the origin of nitrogen isotope fractionation, i.e.,
the mechanism and location of fractionation initiation, remains unclear.

The 14N/15N ratios of various molecules have been quantified in the low-mass
star-forming regions at different evolutionary stages, from the cold dense
cores to the protoplanetary disks. The [14N/15N]elem ratio in the ISM is
$\sim$300 (Lucas & Liszt, 1998; Ritchey et al., 2015; Hily-Blant et al.,
2017). In prestellar cores, nitriles such as HCN, HC3N, and HC5N have a
14N/15N ratio of $\sim$250 – 400, which is close to [14N/15N]elem when the
ratios are derived using the direct method rather than the double-isotope
method (e.g., Taniguchi & Saito, 2017; Hily-Blant et al., 2018; Magalhães et
al., 2018). On the other hand, the observed 14N/15N ratios in N2H+ and N2D+ in
the prestellar cores are higher ($\sim$500 – 1000) than [14N/15N]elem
(Bizzocchi et al., 2013; Furuya et al., 2018; Redaelli et al., 2018). In more
evolved sources, Class 0/I protostellar sources and Class II protoplanetary
disks, the 14N/15N ratio in HCN has been observed (e.g., Guzmán et al., 2017;
Bergner et al., 2020). The 14N/15N ratio in HCN is lower than [14N/15N]elem
and decreases from 290 to 90 with the evolutionary stages traced by bolometric
temperature (Bergner et al., 2020).

Two chemical process have been proposed to explain the observed 15N
fractionation: isotope-exchange reactions (Terzieva & Herbst, 2000) and
isotope-selective photodissociation of N2 (Liang et al., 2007; Heays et al.,
2014). 15N fractionation by isotope-exchange reactions has been thoroughly
studied by many researchers (e.g., Terzieva & Herbst, 2000; Rodgers &
Charnley, 2008; Roueff et al., 2015). A 15N-bearing molecule has a slightly
lower zero-point energy than the corresponding 14N isotopologue. This results
in endothermicity for the exchange of 15N for 14N, leading to the enrichment
of 15N in molecules at low temperatures. However, astrochemical models have
found that under the molecular cloud/core conditions (i.e., at temperatures as
low as 10 K), the isotope-exchange reactions of nitrogen are inefficient
(Roueff et al., 2015; Wirström & Charnley, 2018; Loison et al., 2019) and the
molecular 14N/15N ratios simply reflect the [14N/15N]elem values when
chemistry reaches the steady state (Roueff et al., 2015).

Astrochemical models suggest that isotope-selective photodissociation of N2 is
the more dominant process as compared to the isotope-exchange reactions in
molecular clouds and protoplanetary disks (Furuya & Aikawa, 2018; Visser et
al., 2018). Around N2 photodissociation front, N15N is photodissociated more
than 14N2, which results in the 15N enrichment of atomic nitrogen, and thus
also the other N-bearing molecules produced from atomic nitrogen. The isotope-
selective photodissociation of N2 can explain, at least qualitatively, the
depletion of 14N in N2H+ in prestellar cores (Furuya & Aikawa, 2018) and the
enrichment of 14N in HCN in Class II disks (Visser et al., 2018). As pointed
out in the paper by Bergner et al. (2020), the decrease in the 14N/15N ratio
in HCN with the evolutionary stages could be naturally explained by the
isotope-selective photodissociation of N2. This is because the ultraviolet
(UV) exposure from the central stars increases with the dissipation of the
envelope materials.

The observations, performed using the Atacama Large Millimeter/submillimeter
Array (ALMA), of the TW Hya disk provide the spatial variation of nitrogen
fractionation in gas-phase HCN (Hily-Blant et al., 2019). The 14N/15N ratio in
the comet-forming region ($\sim$20 AU) is 121 $\pm$ 11, which is a factor of
three lower than the [14N/15N]elem ratio (Hily-Blant et al., 2017). This
feature is similar to that recorded in comets, implying that HCN could trace
the fractionated reservoir observed in comets. However, whether the observed
gas-phase isotope ratio reflects the isotope ratio in the ice phase or not
remains unclear. The observed isotope ratio increases in the outward direction
of the disk, which might indicate that the isotope-selective photodissociation
of N2 is a dominant chemical process, as reported based on the disk model by
Visser et al. (2018). However, the disk model is unable to reproduce the
spatial distribution of nitrogen fractionation because in that study, typical
disk models were used rather than a dedicated disk model for the TW Hya disk,
even though the radial distribution of nitrogen fractionation strongly depends
on the physical parameters of the model (Visser et al., 2018). Furthermore,
they used only the initial abundance, similar to the ISM value, whereas the
gas-phase carbon and oxygen are depleted in the TW Hya disk (e.g., Hogerheijde
et al., 2011; Bergin et al., 2013; Favre et al., 2013; Du et al., 2015; Kama
et al., 2016), which could also affect the nitrogen fractionation.

In this study, we have developed a new thermochemical disk model incorporating
the nitrogen isotope chemistry and have attempted to reproduce the nitrogen
fractionation of the gas-phase HCN observed in the TW Hya disk. In Section 2,
we describe the physical and chemical models for the TW Hya disk. The
abundance and isotope profiles obtained using our models have been presented
in Section 3. Our models have been compared with experimental observations,
and the interpretation of the best-fit model and the effects of a few
parameters have been discussed in Section 4. A summary of the conclusions
drawn from this study has been given in Section 5.

##  2 Model

In this study, we have investigated the nitrogen fractionation observed in the
TW Hya disk using a newly developed axisymmetric two-dimensional (2D)
thermochemical disk model, packages of unified modeling for radiative
transfer, gas energetics, and chemistry (PURE-C). This model has been updated
from the thermochemical model for the UV-heated outflow cavity walls (Lee et
al., 2014, 2015). A short summary of the model is provided here and a detailed
description is given in Appendix A.

PURE-C calculates the radiation field and the dust temperature using the Monte
Carlo method for a given density profile. The gas thermal balance and
chemistry are solved iteratively. In this stage, a simple chemical network is
used since only specific species such as H, H2, CO, H2O, electrons, and some
metal species affect the gas thermal balance. The chemistry is evolved with
the full chemical network for 10 Myrs using the obtained density, temperature,
and radiation field profiles.

###  2.1 Parameters for the TW Hya disk

The parameters of the TW Hya disk were adopted from the literature. The
important parameters along with the symbols used in this paper and their
adopted values are listed in Table 1. The density profiles of the gas and dust
grains were adopted from Cleeves et al. (2015), which fits the spectral energy
distribution of the TW Hya disk. In their model, the small dust grains having
radii $r_{g}$= 0.005 $\mu$m – 1 $\mu$m were coupled with gas and the large
dust grains ($r_{g}$= 0.005 $\mu$m – 1 mm) were concentrated near the
midplane. The gas density profile of the disk ($\rho_{g}$) was described in
the following parameterized form in cylindrical coordinates ($R$, $z$)
(Hartmann et al., 1998; Du & Bergin, 2014; Woitke et al., 2016):

$\rho_{g}(R,z)=\frac{\Sigma_{g}(R)}{\sqrt{2\pi}h(R)}\exp\left[-\frac{1}{2}\left(\frac{z}{h(R)}\right)^{2}\right],$
(1)

where the surface density, $\Sigma_{g}(R)$, and the scale height, $h(R)$, are
given by

$\displaystyle\Sigma_{g}(R)$ $\displaystyle=$
$\displaystyle\Sigma_{c}\left(\frac{R}{R_{c}}\right)^{-\gamma}\exp\left[-\left(\frac{R}{R_{\rm
c}}\right)\right],$ (2) $\displaystyle h(R)$ $\displaystyle=$ $\displaystyle
h_{c}\left(\frac{R}{R_{c}}\right)^{\beta}.$ (3)

The gas density profile of our model is shown in the left panel of Figure 1.

The densities of the small ($\rho_{s}$) and large ($\rho_{l}$) dust grain
populations were described by

$\rho_{s}(R,z)=\frac{(1-f)\Sigma_{d}(R)}{\sqrt{2\pi}h(R)}\exp{\left[-\frac{1}{2}\left(\frac{z}{h(R)}\right)^{2}\right]},$
(4)

and

$\rho_{l}(R,z)=\frac{f\Sigma_{d}(R)}{\sqrt{2\pi}\chi_{d}h(R)}\exp{\left[-\frac{1}{2}\left(\frac{z}{\chi_{d}h(R)}\right)^{2}\right]},$
(5)

respectively, where $f$ is the mass fraction of the large dust grains and
$\chi_{d}$ is the concentration factor of the large dust grains. We assumed a
radially constant gas-to-dust ratio of $\delta_{g}=100$, which was defined by
comparing the surface density between the gas and dust grains
($\delta_{g}=\Sigma_{g}/\Sigma_{d}$).

In this study, we ran three models with different mass fractions of the large
dust grain. According to Visser et al. (2018), the nitrogen fractionation of
HCN is affected predominantly by the grain population among the disk
parameters because the grain population affects the UV flux for a given gas
column density and changes the importance of the isotope-selective
photodissociation of N2, which is the dominant process for nitrogen
fractionation. For the TW Hya disk, $f=0.9$ was used by Cleeves et al. (2015)
and $f=0.99$ was used by Kama et al. (2016), where the parameters in the
density profile were also different from those used in Cleeves et al. (2015).
Both models were able to fit the spectral energy distribution of the TW Hya
disk. The model developed by Cleeves et al. (2015) is unable to account for
the observed radial concentration of large dust grains because the gas-phase
species depend mainly on the total surface area of the small dust grains in
the atmosphere, which are responsible for the observed molecular emission
lines. Therefore, we adopted the simplest model of Cleeves et al. (2015) and
ran the models with two $f$ values of 0.9 (ML) and 0.99 (MS). We also
introduced a hybrid model (MH) in order to satisfactorily fit the observations
(see Sect. 4.1);

$f=\left\\{\begin{array}[]{cl}0.9,&\;if\,\,R\leq 20{\rm~{}AU}\\\
0.9+0.09\times(R-20{\rm~{}AU})/20{\rm~{}AU},&\;if\,\,20{\rm AU}\leq R\leq
40{\rm~{}AU}\\\ 0.99,&\;if\,\,R>40{\rm~{}AU}.\end{array}\right.$ (6)

It should be noted here that the three models exhibit similar spectral energy
distributions, except for the flux around 200 $\mu$m with the maximum
difference of a factor of two.

###  2.2 Chemistry for the nitrogen fractionation of HCN in the TW Hya disk

Certain features experimentally observed in the disk were taken into account
in PURE-C. In the midplane of the disk, chemical reactions on the dust grains
are important for forming complex organic molecules (e.g., Furuya & Aikawa,
2014; Walsh et al., 2014). The line-shielding effect of molecules by H, H2, C,
and CO was also considered in the model by using the tables provided by Heays
et al. (2017), because UV photons can be attenuated by abundant molecules as
well as the dust grains present in the atmosphere of the disk. The column
densities used in self-shielding and line-shielding were calculated in 2D
spaces by averaging the radial and vertical column densities weighted by the
UV fluxes along each direction. X-ray ionization is another important process
occurring in the disk (Cleeves et al., 2015; Offner et al., 2019), and thus we
adopted an analytical model to calculate the X-ray ionization rate in the disk
(Igea & Glassgold, 1999). We assumed that X-ray ionization works similar to
cosmic rays in order to incorporate its chemistry (see Walsh et al., 2012),
and the cosmic ray ionization rate was taken to be 1.0$\times 10^{-19}$ s-1
(Cleeves et al., 2015).

The nitrogen isotope chemistry was also included in PURE-C. We reduced the
chemical network used in Furuya & Aikawa (2018), which is modeled for nitrogen
fractionation in the molecular cloud. The original network included the spin
state of H2, the isotope-selective photodissociation of N2 (Heays et al.,
2014), and the isotope-exchange reactions (Roueff et al., 2015). However, this
network is too large to solve the 2D disk model. Thus, we have focused only on
the nitrogen fractionation of HCN in this work. Visser et al. (2018) used the
2D disk model for nitrogen isotope chemistry with a small network including
species listed in Table B.1 of their paper, which were sufficient for
investigating the nitrogen fractionation of HCN. Therefore, we reduced the
chemical network of Furuya & Aikawa (2018) by choosing the species used in
Visser et al. (2018) except for the polycyclic aromatic hydrocarbon (PAH)
species and 13C-bearing species. It should be noted here that PURE-C adopts a
two-phase model (gas and grain surface) using the binding energy on water ice,
whereas the model of Furuya & Aikawa (2018) adopted the three-phase model
(gas, grain surface, and ice mantle) using the binding energy corrected by the
composition of the ice surface. In our model, the binding energies (on water
ice) of N2, HCN, CN, and NH3 were 1000 K, 3700 K, 2800 K, and 5500 K
respectively.

The stellar spectra of the TW Hya disk (e.g., Nomura & Millar, 2005; France et
al., 2014) were used for obtaining the unshielded photodissociation/ionization
rate ($\alpha_{0}$) and the line-shielding function (Heays et al., 2017). The
dust shielding function, which is generally calculated for an infinite slab,
is inadequate for our 2D disk model that considers the settling of dust
grains. We calculated the dust-attenuated UV fluxes ($\chi$) at one
representative wavelength of 9.8 eV by solving the UV radiative transfer, and
the photodissociation/ionization rate was simply calculated as
$\alpha_{0}\chi$. For the CO (Visser et al., 2009) and N2 (Heays et al., 2014)
self-shielding functions, we chose the excitation temperatures and Doppler
widths of 20 K and 0.3 km s-1 for CO and 30 K and 0.13 km s-1 for N2
respectively.

In this study, we have investigated the effect of the initial elemental
abundances of carbon and oxygen. Elemental gas-phase carbon and oxygen
deficiencies of a factor of 10–100 and their abundance ratios, [C/O]elem$>$ 1,
were already inferred for the TW Hya disk (Hogerheijde et al., 2011; Bergin et
al., 2013; Favre et al., 2013; Du et al., 2015; Kama et al., 2016). These
deficiencies could be understood as a chemical and grain evolution (Kama et
al., 2016; Bergin et al., 2016): CO and H2O freeze-out onto the dust grains
and the grains covered with ices settle down towards the midplane and drift
inwards. Chemical reactions also lead to the depletion of the gas-phase CO
(e.g., Walsh et al., 2010; Furuya & Aikawa, 2014; Schwarz et al., 2018, 2019).
We adopted [14N/15N]elem = 330 for the TW Hya disk (Hily-Blant et al., 2017).
The [C/O]elem ratio affects the HCN column density (Cleeves et al., 2018) and
could also possibly affect the nitrogen fractionation of HCN. There is no
study available on the elemental gas-phase nitrogen abundance or the [N/O]elem
ratio in the TW Hya disk, and the gas-phase nitrogen element does not deplete
in the IM Lup protoplanetary disk (Cleeves et al., 2018). Therefore, we
modified the initial abundance values reported by Cleeves et al. (2015) and
ran four models with different initial abundances listed in Table 2. The I1
model represents the typical values in the ISM, and the I2 model reduces the
abundances of water ice and CO gas by two orders of magnitude. In the I3 and
I4 models, the initial abundances of water ice and atomic carbon were changed
so that the elemental abundances in the gas phase satisfy [C/O]elem$>$ 1 after
the desorption of water ice. Therefore, in total, we ran 12 models listed in
Table 3 using different grain populations (ML, MS, MH) and different initial
abundances (I1–I4).

##  3 Results

###  3.1 Physical characteristics

Figure 1 shows the 2D distributions of the gas density ($n_{\rm gas}$) and the
ionization rate ($\zeta$), and Figure 2 presents the profiles of the UV flux
($\chi$, in the units of Draine field), dust temperature, gas temperature, and
N2 column density from the central star to a given position (hereafter $N_{\rm
N_{2}}$)111 In this work, two types of column densities have been used. One
($N_{\rm N_{2}}$) is measured from the central star to a given position and
used for calculating the self- and line-shielding effects (see Figures 2 and
3, and Appendix). The other is integrated vertically to compare with
observations, and is the default column density used in this study (see
Figures 8 and 9). for the MLI3, MSI3, and MHI3 models. Figure 3 shows the
vertical cuts of Figures 1 and 2 for the MLI3 model at three radii of 20, 40,
and 60 AU. The column density ratio of HC14N/HC15N observed in the TW Hya disk
(Hily-Blant et al., 2019) is the lowest at 20 AU and peaks around 40 AU. The
60 AU point is just outside the radius where the ratio is measured (55 AU).

Grain populations affect the UV fluxes and the dust and gas temperatures.
Small dust grains in the atmosphere in the MS model are observed to be
depleted by an order of magnitude compared to those in the ML model.
Therefore, UV photons are attenuated less in the MS model than in the ML
model. Photons in the other wavelengths from the central star are also
observed to be attenuated less in the MS model, and thus a higher radiation
flux in the MS model produces a higher dust temperature than those in the ML
model. At a height of $z/R<$ 0.2, collision between the gas and dust grains is
a dominant heating and cooling source. Thus, the gas temperature also exhibits
a trend similar to the dust temperature.

The grain populations also influence $N_{\rm N_{2}}$ as shown in the bottom
panels in Figure 2. $N_{\rm N_{2}}$ depends on the path of the UV photons as
well as the N2 number density (abundance) distribution. The white dotted lines
in the bottom panels in Figure 2 indicate the height at which the UV fluxes
directly irradiated from the central star are the same as those scattered by
the dust grains in the atmosphere of the disk and descending vertically. Most
UV photons pass via the radial direction from the central star above the white
dotted lines. Although the N2 number density distribution is also affected by
$N_{\rm N_{2}}$, the latter tends to be larger for UV photons passing through
the radial direction than for those passing via the vertical direction in our
models. The height of the white dotted line in the MS model is lower than that
in the ML model because the small dust grains are less abundant in the former
as compared to the latter. This causes less scattering of the UV photons in
the atmosphere in the MS model than in the ML model. Therefore, for a given
radius and height, the MS model has a higher $N_{\rm N_{2}}$ than the ML
model, and near the midplane in the outer disk, the MS model has a lower N2
photodissociation rate than the ML model even though the UV flux is higher in
the former than in the latter. For example, at a radius of 60 AU and a height
of $z/R$= 0.2, the MSI3 model has a UV flux and $N_{\rm N_{2}}$ that are a
factor of five higher than the MLI3 model. As a result, the N2
photodissociation rate in the MSI3 model is observed to be a factor of five
lower than that in the MLI3 model, indicating that the N2 self-shielding
effect is dominant than the dust attenuation of UV photons (see the dotted and
solid lines in the top-right panel of Figure 3).

###  3.2 Cyanide abundances and isotope ratios

Figures 4 and 5 show the distributions of the abundance and nitrogen isotope
ratio respectively for atomic N, N2, HCN, and CN in the MLI3 model and the
corresponding vertical cuts at 20, 40, and 60 AU are plotted in Figure 6.
Figure 7 plots the same quantities as Figure 6 except for the MSI3 model. The
gas-phase N2 is observed to be abundant in the warm molecular layer ($z/R\geq
0.1$ and UV flux $\chi\geq 1$), and both CN- and HCN-abundant layers are
observed to be located within the N2 abundant layer. An upper boundary of the
N2 abundant layer is confined by the N2 photodissociation front. As for the
lower boundary, all N-bearing species freeze-out on the dust grains. The
evaporation temperature of N2 is around 20 K in our models, and the freeze-out
onto the dust grains depletes the gas-phase N2 near the midplane in the outer
disk ($>$20 AU). The gas-phase N2 is also depleted above the N2 snowline with
the dust temperature of $\sim$20 K, as shown in the top-right plot in Figure
4. In this region, the ionization rate is higher than $\sim$10-17 s-1 (see the
right panel of Figure 1 and the bottom-left panel of Figure 3), and the gas-
phase CO is converted to CO2 ice within a few Myrs (Furuya & Aikawa, 2014;
Bergin et al., 2014). The gas-phase N2 becomes NH3 ice after the gas-phase CO
depletes (see Furuya & Aikawa, 2014). The lower boundary of the N2-abundant
layer moves from the dense region closer to the midplane to the less dense
region in the upper layer with the passage of time. Notably, in our models
with [C/O]elem$>1$, even in the warm molecular layer, HCN ice is more abundant
than the gas-phase HCN and it is the most dominant ice among the N-bearing
species.

Heavy atomic nitrogen (15N) is enriched near the atomic N and N2 transition
layer (hereafter, NT layer), where the atomic N and N2 abundances with respect
to the total hydrogen nuclei ($n_{\rm gas}$ = $n$(H) + 2 $n$(H2)) are higher
than $\sim$10-7 due to the isotope-selective photodissociation of N2.
Formation of CN and HCN generally begins with atomic N, and thus the nitrogen
isotope ratios of the two molecules follow the atomic 14N/15N ratio near the
NT layer. These results are similar to those reported by Visser et al. (2018).

The HCN-abundant layer could be composed of two layers: (i) an NT layer where
the HC14N/HC15N ratio follows the atomic 14N/15N ratio and (ii) a “lower
molecular layer” (hereafter LM layer) just below the NT layer with an
abundance of the atomic N lower than $\sim$10-7. In this layer, the isotope-
exchange reactions affect the HC14N/HC15N ratio. In our model, we included the
reaction of 15N + CN $\leftrightarrow$ 14N + C15N, which was ignored in the
work of Visser et al. (2018). For the rate coefficient of this reaction, we
used the upper limit value proposed by Roueff et al. (2015). As can be seen in
the bottom panels of Figure 6, the atomic 14N/15N ratio is higher than
[14N/15N]elem, whereas the HC14N/HC15N ratio is lower than [14N/15N]elem in
the LM layer. This is due to the abovementioned isotope-exchange reactions.
However, the HC14N/HC15N ratio in the LM layer is closer to the [14N/15N]elem
ratio compared to that in the NT layer and the isotope-selective
photodissociation of N2 is the dominant process for nitrogen fractionation of
HCN. The CN-abundant layer also exhibits a trend similar to the HCN-abundant
layer. Furthermore, the CN-abundant layer has an additional “upper molecular
layer” just above the NT layer, where the C14N/C15N ratio is close to the
[14N/15N]elem ratio because the $N_{\rm N_{2}}$ is too low, and thus the
isotope-selective photodissociation of N2 does not work.

###  3.3 Effect of the initial abundance and the grain population

The column densities of HCN, CN, and C2H are affected more by the [C/O]elem
ratio than by the elemental abundances of carbon and oxygen and the dust grain
population. Figure 8 shows the column densities of HCN, CN, and C2H and the
column density ratios of HC14N/HC15N and C14N/C15N obtained from our models
with different grain populations and initial abundances (gas-phase elemental
abundances and the [C/O]elem ratio). The column density of C2H is more
sensitive to the [C/O]elem ratio than those of HCN and CN (Cleeves et al.,
2018). Further, the models with [C/O]elem$>$1 can reproduce the
observationally detectable column density of C2H, as reported by Bergin et al.
(2016). The bright C2H emission in the protoplanetary disk is used as an
indicator of [C/O]elem$>$1, and thus we also checked the effect of the initial
abundance and grain population on the C2H column density and compared them
with observations in this study (see Section 4.1). The HCN and CN column
densities for the models with [C/O]elem$>$ 1 are a factor of 10 – 100 higher
than those for the models with [C/O]elem$<$ 1 (Cleeves et al., 2018). For the
models with [C/O]elem$>$ 1, the MS models show a factor of 2 – 10 higher
column densities of HCN, CN, and C2H than the ML models. The UV photons
penetrate closer to the midplane of the disk in the MS model than in the ML
model and this affects the gas-phase CO and N2 depletion processes. When the
models have the same [C/O]elem ratio, the gas-phase elemental abundances of
carbon and oxygen do not affect the HCN and CN column densities, whereas the
C2H column density obtained from the MLI4 model is an order of magnitude
higher than that obtained from the MLI3 model (see the green and red lines in
Figure 8).

On the other hand, the column density ratios of HC14N/HC15N and C14N/C15N are
affected by the grain population as well as the [C/O]elem ratio. The column
density ratio of the isotopologues is determined by the contribution of the NT
layer relative to that of the LM layer. The NT layer reduces the column
density ratio of the isotopologues, whereas the LM layer leads to the ratio to
be close to the [14N/15N]elem ratio. For HCN, the models with [C/O]elem$>$ 1
have higher ratios than the models with [C/O]elem$<$ 1\. In the models with
[C/O]elem$<$ 1, more atomic O exists in the gas phase and atomic N becomes NO
rather than HCN and CN. This effect is found to be greater in the LM layer
than in the NT layer, and thus the contribution of the NT layer to the total
column density in the model with [C/O]elem$>$ 1 is smaller than in the model
with [C/O]elem$<$ 1\. The contribution of the “upper molecular” layer to the
total CN column density is small in our models except for the models having an
initial abundance of I1. Therefore, the I1 models have a column density ratio
of C14N/C15N that is higher or comparable to that in the models with
[C/O]elem$>$ 1.

The column density ratios of HC14N/HC15N and C14N/C15N obtained from the MS
model are closer to [14N/15N]elem than those obtained from the ML model. As
mentioned above, the UV photons penetrate closer to the midplane of the disk
in the MS model as compared to that in the ML model, and the LM layer is
extended to a lower height in the MS model than in the ML model (see Figures 6
and 7). Furthermore, $N_{\rm N_{2}}$ in the NT layer increases faster from the
disk atmosphere to the deeper regions in the MS model than in the ML model, as
shown in the bottom panels in Figure 2, because the UV photons in the NT layer
tend to propagate radially in the MS model, whereas they propagate vertically
in the ML model (see the dotted vertical lines in Figures 6 and 7). Therefore,
the isotope-selective photodissociation of N2 works in the narrower height
ranges in the MS model than in the ML model. This indicates that the width of
the NT layer in the MS model is narrower than that in the ML model (see
Figures 6 and 7). As a result, the contribution of the NT layer in the MS
model is smaller than that in the ML model, and thus the isotope ratios are
closer to [14N/15N]elem in the MS model than in the ML model (see Figure 8).
It should be noted that this trend is contrary to the result obtained by
Visser et al. (2018) where the model having a larger fraction of mass in the
large dust grains has a lower column density ratio of HC14N/HC15N than the
model having a smaller fraction of mass in the large dust grains (see Fig. 14
in Visser et al. (2018)). The main difference between our model and theirs is
the method of calculating the column density under the self-shielding effect.
In their model, the column density is measured as the minimum of the
radial/vertical (inward/upward) column density (see also Miotello et al.,
2014). On the other hand, the column density weighted by the UV flux in the
radial and vertical directions has been used in our model (see Appendix for
more details). Therefore, the value of $N_{\rm N_{2}}$ obtained using our
model tends to be higher than that obtained in their work and depends on the
path of the UV photons, which could result in opposite trends in the two
models.

##  4 Discussion

###  4.1 Comparison with observations

A comparison of the results obtained from our models has been done with the
observations on the TW Hya disk. The observed column densities derived from
the single dish observations (gray bars, Kastner et al., 2014) and the ALMA
observations (see below) are plotted in Figure 8. The observed isotope column
density ratios of HCN and CN were adopted from (Hily-Blant et al., 2017,
2019), plotted as filled circles in the 2nd row and gray hatched bars in the
4th row, respectively.

Kastner et al. (2014) carried out a line survey around 300 GHz on the TW Hya
disk using the 12 m Atacama pathfinder experiment (APEX) telescope. The column
densities of CN and C2H derived from their hyperfine analyses were (9.6 $\pm$
1.0) $\times 10^{13}$ cm-2 and (5.1 $\pm$ 3.0) $\times 10^{15}$ cm-2
respectively, when averaged within 5 arcsec. The HCN column density of (0.3 –
9.3) $\times 10^{13}$ cm-2 was derived under the assumption of local thermal
equilibrium with an excitation temperature similar to those of the CN and C2H
lines (5 K – 13 K).

ALMA observations and modeling with the HCN, CN, and C2H lines was carried out
for the TW Hya disk. The CN and C2H lines were observed to exhibit ring-like
emission distributions with emission peaks around 40 and 60 AU, respectively
(Bergin et al., 2016; Cazzoletti et al., 2018; Hily-Blant et al., 2017;
Kastner et al., 2015; Nomura et al., 2016). In particular, [C/O]elem$>$ 1 was
required to reproduce bright C2H emission rings (Bergin et al., 2016). The
HC15N ALMA observation indicated that the line intensity declines towards
outer disk (Hily-Blant et al., 2019). When gas temperatures of 30 K – 50 K and
gas density of 106 – 108 cm-3 (Bergin et al., 2016; Hily-Blant et al., 2019)
were adopted for the RADEX model (van der Tak et al., 2007), the observed
fluxes of the optically thin lines matched with the column densities of HCN,
CN, and C2H of 1014 – 1015 cm-2 (see the black bars in the 1st row, the gray
hatched bars in the 3rd row, and the black bars in the 5th row, respectively,
in Figure 8). Here, a brightness temperature of 1 K was used for HC15N 4-3 at
40 AU (Hily-Blant et al., 2019), the flux of 0.166 Jy km s-1 was integrated
over two C15N N=3-2, J=7/2$-$5/2 lines (Hily-Blant et al., 2017) and the peak
flux of 0.089 Jy beam-1 for the C2H N = 4 $-$ 3, J = 7/2 $-$ 5/2, F = 3 $-$ 3
line (Bergin et al., 2016), respectively.

We note that the molecular column densities become higher under the assumption
of low gas density. The HCN column densities at 40 AU are $\sim$1014 cm-2 and
up to $\sim$1015 cm-2 when the gas density is higher than 107 cm-3 and 106
cm-3, respectively. The SMA observation of the HCN 3–2 line of the TW Hya disk
showed that the HCN column density at 40 AU was $\sim$6$\times 10^{13}$ cm-2
(Qi et al., 2008). The HCN-abundant layer had a gas density higher than 107
cm-3 and a gas temperature around 30 K in the model used in Qi et al. (2008).

The column densities constrain the initial gas-phase elemental abundances of
carbon and oxygen. As shown in the bottom panels of Figure 8, the C2H column
density is most sensitive to the [C/O]elem ratio and only those models having
[C/O]elem$>$ 1 can reproduce the observed column density reported by Bergin et
al. (2016). Furthermore, the models with [C/O]elem$>$ 1 also reproduce the
column density derived from the HC15N ALMA observations (see the black bars in
the top panels of Figure 8, Hily-Blant et al., 2019). However, the observed CN
column densities (the dark gray horizontal bars in the middle rows of Figure
8, Hily-Blant et al., 2017) are a factor of 10 lower than those from the
models with [C/O]elem$>$ 1.

The ML and MS models were unable to reproduce the observed column density
ratio of HC14N/HC15N, and thus the MH model with Equation 6 was introduced to
satisfactorily fit the observations. In the observations (Hily-Blant et al.,
2018), the column density ratio of HC14N/HC15N increases from $\sim$120 up to
$\sim$330 as the distance from the central star increases from 20 AU to 50 AU.
In our models with the ML and MS grain populations and [C/O]elem$>$ 1, the
column density ratios of HC14N/HC15N are roughly constant or decrease with the
radius in the outer disk ($>$20 AU). The observed HC14N/HC15N ratio is similar
to the ratio around 25 AU in the MLI3 model and to that around 40 AU in the
MSI3 model. Therefore, we can infer that the hybrid model and the column
density ratio of HC14N/HC15N in the MHI3 model are consistent with the
observations, as seen from Figure 8.

The MHI3 model could reproduce the observed ring-like emission distribution of
C2H. The C2H emission peaks at 60 AU and the central emission hole has an
integrated intensity that is a factor of six weaker than that at the ring peak
(Bergin et al., 2016). Their models with the spatially varied [C/O]elem ratio
could reproduce this ring-like emission (Bergin et al., 2016). The C2H column
densities in the MLI3 and MSI3 models were spatially constant in the outer
disk ($>$ 10 AU). On the other hand, the results of the MHI3 model showed that
the column density of C2H peaks around 40 AU and column density of C2H at 40
AU is a factor of 10 higher than that in the inner disk ($<$ 20 AU).
Therefore, the model with the spatially varied grain population could be
another solution for the ring-like C2H emission. However, it needs additional
fine-tuning.

###  4.2 The best-fit model

The grain population of the best fit model can be understood using grain
evolution models. As mentioned above, the small dust grains in the atmosphere
control the UV flux and $N_{\rm N_{2}}$ in the HCN-abundant layer, and change
the HC14N/HC15N ratio. Therefore, the best fit model indicates that the small
dust grains in the atmosphere of the outer disk ($>$ 20 AU) are depleted more
than those of the inner disk. However, it is difficult to obtain the amount of
small dust grains from the infrared observations of the (optically thick)
light scattered off the dust grains toward TW Hya (for example, van Boekel et
al., 2017). The small dust grains can be replenished when grain fragmentation
is efficient and dominates the grain growth (e.g., Dullemond & Dominik, 2005;
Birnstiel et al., 2012). Generally, the inner disks are dominated by
fragmentation, whereas the outer disks are growth-dominated and the small dust
grains can survive for a few Myrs in the case of strong turbulence (e.g.,
Birnstiel et al., 2012). The grains become larger due to collisional sticking
in the outer disk. They migrate to the inner disk and fragment into small dust
grains if the turbulent motion is sufficiently strong. Furthermore, this model
can also explain the multiwavelength observations toward the TW Hya disk (Menu
et al., 2014).

The [C/O]elem ratio also affects the HC14N/HC15N ratio. In this study, we
assumed a spatially constant [C/O]elem ratio. However, the infrared molecular
line observations were interpreted using a model with a [C/O]elem of 0.5
within the disk radius of 2.4 AU (Bosman & Banzatti, 2019), whereas the
molecular line observations from ALMA were explained using the models with
[C/O]elem$>$ 1 in the outer disk (e.g., Bergin et al., 2016, this work).
Furthermore, models considering both grain and chemical evolution also showed
that the [C/O]elem ratio varies in both radial and vertical directions (Krijt
et al., 2018, 2020). When the [C/O]elem ratio varies only along the distance
from the central star, the column densities of HCN and C2H change
significantly (Cleeves et al., 2018). In the best fit model, the HC14N/HC15N
ratio at 20 AU is slightly higher than the ratio obtained from experimental
observations. When the [C/O]elem ratio increases with the height, the
contribution of the NT layer to the total HCN column density can be increased.
Then, the HC14N/HC15N ratio can be fitted satisfactorily. Furthermore, our
model is unable to reproduce the column density and the nitrogen isotope ratio
in CN, which might be solved by vertically varying the [C/O]elem ratio and
small [C/O]elem ratio in the surface layer of the disk.

###  4.3 Effect of other parameters

We have investigated the effects of a few parameters on the results obtained
from our model. The line-shielding effect has been included in this work and
this effect could work in the disk atmosphere where small dust grains are
depleted due to, for example, growth of the dust grains and their settling in
the midplane of the disk. The isotope-exchange reaction and the vibrationally
excited H2 also affect the HC14N/HC15N ratio. We ran the MHI3 model without
the line-shielding effect, the isotope-exchange reaction, and the
vibrationally excited H2 in order to investigate their effect on the resulting
molecular abundances and nitrogen fractionation. For the X-ray ionization
rate, we used the analytical model of Igea & Glassgold (1999). However,
depletion of small dust grains due to dust settling could reduce the X-ray
opacity (Bethell & Bergin, 2011), resulting in an increased X-ray ionization
rate near the midplane of the disk. We also ran the model using the X-ray
ionization rate a factor of ten higher than that in the reference model. For
the sake of simplicity, we assumed that the gas temperature is the same as
that in the reference model, MHI3. The results thus obtained are plotted in
Figure 9.

The line-shielding effect increases the column densities of HCN, CN, and C2H,
but it does not change the column density ratio of the isotopologues
significantly. When the line-shielding effect is considered, a larger number
of molecules can survive the photodissociation process in the entire warm
molecular layer. This effect is not significantly different for the NT layer
and the LM layer, and thus the contributions of both layers to the total
column density remain unchanged. Furthermore, the isotope ratios of HCN and CN
in the NT layer exhibit a trend similar to atomic N, which is affected by the
N2 self-shielding effect and not by the line-shielding effect.

The isotope-exchange reaction leads to only a minor change in the column
density ratio of the isotopologues. As mentioned above, we included the
reaction 15N + CN $\leftrightarrow$ 14N + C15N, which was ignored in the study
by Visser et al. (2018). However, the isotope-selective photodissociation of
N2 predominantly affects the column density ratio of the isotopologues rather
than the isotope-exchange reactions, as reported by Visser et al. (2018), even
if the abovementioned exchange reaction is included. However, in some cases,
the isotope-exchange reaction contributes up to 60 %: around 40 AU, the
difference between the column density ratio of HCN from [14N/15N]elem in the
reference model is 50, whereas that in the model without the abovementioned
exchange reaction is 30 (it corresponds to the difference between the colored
lines and the gray horizontal lines in the 2nd row in Figure 9).

The vibrationally excited H2 increases the CN and HCN column densities and
reduces the nitrogen isotope ratio of both molecules. The endothermic reaction
between the vibrationally excited H2 and atomic N powers the formation of CN
and HCN (e.g., Visser et al., 2018). This reaction is important in the NT
layer, but not in the LM layer. Thus, when the vibrationally excited H2 is
considered, the contribution of the NT layer to the total column density
increases and the column density ratio of the isotopologues is reduced
further.

X-ray ionization is important in the LM layer and increases the column density
ratios of the isotopologues as well as the column densities of CN and HCN.
Higher ionization increases the abundances of assorted ions, resulting in
increasing abundances of molecules other than CO and N2 (Aikawa & Herbst,
1999). UV photons are the dominant ionization sources rather than X-rays in
the NT layer (above $z/R>$ 0.2), where the isotope-selective photodissociation
of N2 reduces the nitrogen isotope ratio of atomic N, CN, and HCN. Thus, in
this layer, there is no difference between the models incorporating different
X-ray ionization rates. On the other hand, in the LM layer, X-ray photons form
a larger number of ions, and thus larger quantities of CN and HCN. Thus, when
the X-ray ionization rate increases, most CN and HCN molecules exist in the LM
layer and both the column density and the isotope ratio increase.

##  5 Summary

We have developed a disk model for reproducing the nitrogen fractionation of
the gas-phase HCN observed in the TW Hya disk. In the model, the isotope-
exchange reactions and isotope-selective photodissociation of N2 were
included, which are important reactions that cause nitrogen fractionation. The
effect of the grain population and initial abundances on the nitrogen
fractionation process and the column density of HCN was investigated in this
study.

The observed column density and the nitrogen isotope ratio of HCN could be
reproduced when the [C/O]elem ratio was larger than unity and the small dust
grains in the atmosphere were depleted considerably in the outer disk than in
the inner disk. The HC14N/HC15N ratio was lower than 100 around the NT layer
due to the isotope-selective photodissociation of N2. The column density ratio
of HC14N/HC15N depended on the relative contribution of the NT layer to the
total HCN column density. The higher [C/O]elem ratio increased the HCN column
density as well as the column density ratio of HC14N/HC15N because it led to a
broader LM layer for HCN. The depletion of small dust grains in the atmosphere
induced a narrower NT layer and a broader LM layer, which resulted in an
increase the HCN column density and the column density ratio of HC14N/HC15N.
The radially constant grain population models showed that the column density
ratio of HC14N/HC15N is spatially constant or decreases outwards. Thus, the
observed HC14N/HC15N ratio can be explained by a model in which the small dust
grains in the outer disk are depleted more than those in the inner disk. This
population of small dust grains is consistent with grain evolution models
(e.g., Birnstiel et al., 2012). In this work, the line-shielding effects for
photodissociation by stellar UV photons were included, which enhanced the
column density of HCN, but did not change the column density ratio of
HC14N/HC15N.

We would like to thank the referee for comments that improved our paper. This
work is supported by the NAOJ leadership program, MEXT Grants-in-Aid for
Scientific Research 18H05441 and 19K03910, the Basic Science Research Program
through the National Research Foundation of Korea (grant no.
NRF-2018R1A2B6003423), and NAOJ ALMA Scientific Research Grant Number of
2018-10B.

## References

  * Abbas et al. (2006) Abbas, M. M., Tankosic, D., Craven, P. D., et al. 2006, ApJ, 645, 324, doi: 10.1086/504281
  * Aikawa & Herbst (1999) Aikawa, Y., & Herbst, E. 1999, A&A, 351, 233
  * Bai & Goodman (2009) Bai, X.-N., & Goodman, J. 2009, ApJ, 701, 737, doi: 10.1088/0004-637X/701/1/737
  * Bergin et al. (2014) Bergin, E. A., Cleeves, L. I., Crockett, N., & Blake, G. A. 2014, Faraday Discussions, 168, 61, doi: 10.1039/C4FD00003J
  * 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
  * Bergin et al. (2013) Bergin, E. A., Cleeves, L. I., Gorti, U., et al. 2013, Nature, 493, 644, doi: 10.1038/nature11805
  * Bergner et al. (2020) Bergner, J. B., Öberg, K. I., Bergin, E. A., et al. 2020, ApJ, 898, 97, doi: 10.3847/1538-4357/ab9e71
  * Bertin et al. (2012) Bertin, M., Fayolle, E. C., Romanzin, C., et al. 2012, Physical Chemistry Chemical Physics (Incorporating Faraday Transactions), 14, 9929, doi: 10.1039/C2CP41177F
  * Bertin et al. (2013) —. 2013, ApJ, 779, 120, doi: 10.1088/0004-637X/779/2/120
  * Bethell & Bergin (2011) Bethell, T. J., & Bergin, E. A. 2011, ApJ, 740, 7, doi: 10.1088/0004-637X/740/1/7
  * Birnstiel et al. (2012) Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148, doi: 10.1051/0004-6361/201118136
  * Bizzocchi et al. (2013) Bizzocchi, L., Caselli, P., Leonardo, E., & Dore, L. 2013, A&A, 555, A109, doi: 10.1051/0004-6361/201321276
  * Bjorkman & Wood (2001) Bjorkman, J. E., & Wood, K. 2001, ApJ, 554, 615, doi: 10.1086/321336
  * Bosman & Banzatti (2019) Bosman, A. D., & Banzatti, A. 2019, A&A, 632, L10, doi: 10.1051/0004-6361/201936638
  * Bruderer et al. (2012) Bruderer, S., van Dishoeck, E. F., Doty, S. D., & Herczeg, G. J. 2012, A&A, 541, A91, doi: 10.1051/0004-6361/201118218
  * Burke & Hollenbach (1983) Burke, J. R., & Hollenbach, D. J. 1983, ApJ, 265, 223, doi: 10.1086/160667
  * 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
  * Ceccarelli et al. (2014) Ceccarelli, C., Caselli, P., Bockelée-Morvan, D., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 859
  * 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. (2015) Cleeves, L. I., Bergin, E. A., Qi, C., Adams, F. C., & Öberg, K. I. 2015, ApJ, 799, 204, doi: 10.1088/0004-637X/799/2/204
  * 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
  * Dorschner et al. (1995) Dorschner, J., Begemann, B., Henning, T., Jaeger, C., & Mutschke, H. 1995, A&A, 300, 503
  * Du & Bergin (2014) Du, F., & Bergin, E. A. 2014, ApJ, 792, 2, doi: 10.1088/0004-637X/792/1/2
  * Du et al. (2015) Du, F., Bergin, E. A., & Hogerheijde, M. R. 2015, ApJ, 807, L32, doi: 10.1088/2041-8205/807/2/L32
  * Dullemond & Dominik (2005) Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971, doi: 10.1051/0004-6361:20042080
  * Favre et al. (2013) Favre, C., Cleeves, L. I., Bergin, E. A., Qi, C., & Blake, G. A. 2013, ApJ, 776, L38, doi: 10.1088/2041-8205/776/2/L38
  * Fayolle et al. (2011) Fayolle, E. C., Bertin, M., Romanzin, C., et al. 2011, ApJ, 739, L36, doi: 10.1088/2041-8205/739/2/L36
  * Fayolle et al. (2013) —. 2013, A&A, 556, A122, doi: 10.1051/0004-6361/201321533
  * 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
  * Füri & Marty (2015) Füri, E., & Marty, B. 2015, Nature Geoscience, 8, 515, doi: 10.1038/ngeo2451
  * Furuya & Aikawa (2014) Furuya, K., & Aikawa, Y. 2014, ApJ, 790, 97, doi: 10.1088/0004-637X/790/2/97
  * Furuya & Aikawa (2018) —. 2018, ApJ, 857, 105, doi: 10.3847/1538-4357/aab768
  * Furuya et al. (2018) Furuya, K., Watanabe, Y., Sakai, T., Aikawa, Y., & Yamamoto, S. 2018, A&A, 615, L16, doi: 10.1051/0004-6361/201833607
  * Garrod (2013) Garrod, R. T. 2013, ApJ, 765, 60, doi: 10.1088/0004-637X/765/1/60
  * Garrod & Pauly (2011) Garrod, R. T., & Pauly, T. 2011, ApJ, 735, 15, doi: 10.1088/0004-637X/735/1/15
  * Geers et al. (2006) Geers, V. C., Augereau, J.-C., Pontoppidan, K. M., et al. 2006, A&A, 459, 545, doi: 10.1051/0004-6361:20064830
  * Groenewegen (1994) Groenewegen, M. A. T. 1994, A&A, 290, 531
  * Guzmán et al. (2017) Guzmán, V. V., Öberg, K. I., Huang, J., Loomis, R., & Qi, C. 2017, ApJ, 836, 30, doi: 10.3847/1538-4357/836/1/30
  * Hartmann et al. (1998) Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385, doi: 10.1086/305277
  * Hasegawa et al. (1992) Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167, doi: 10.1086/191713
  * 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
  * Heays et al. (2014) Heays, A. N., Visser, R., Gredel, R., et al. 2014, A&A, 562, A61, doi: 10.1051/0004-6361/201322832
  * Hily-Blant et al. (2018) Hily-Blant, P., Faure, A., Vastel, C., et al. 2018, MNRAS, 480, 1174, doi: 10.1093/mnras/sty1562
  * Hily-Blant et al. (2017) Hily-Blant, P., Magalhaes, V., Kastner, J., et al. 2017, A&A, 603, L6, doi: 10.1051/0004-6361/201730524
  * Hily-Blant et al. (2019) Hily-Blant, P., Magalhaes de Souza, V., Kastner, J., & Forveille, T. 2019, A&A, 632, L12, doi: 10.1051/0004-6361/201936750
  * Hogerheijde et al. (2011) Hogerheijde, M. R., Bergin, E. A., Brinch, C., et al. 2011, Science, 334, 338, doi: 10.1126/science.1208931
  * Igea & Glassgold (1999) Igea, J., & Glassgold, A. E. 1999, ApJ, 518, 848, doi: 10.1086/307302
  * Ishimaru (1978) Ishimaru, A. 1978, Wave propagation and scattering in random media. Volume I \- Single scattering and transport theory
  * Kama et al. (2016) Kama, M., Bruderer, S., van Dishoeck, E. F., et al. 2016, A&A, 592, A83, doi: 10.1051/0004-6361/201526991
  * Kamp & van Zadelhoff (2001) Kamp, I., & van Zadelhoff, G.-J. 2001, A&A, 373, 641, doi: 10.1051/0004-6361:20010629
  * Kastner et al. (2014) Kastner, J. H., Hily-Blant, P., Rodriguez, D. R., Punzi, K., & Forveille, T. 2014, ApJ, 793, 55, doi: 10.1088/0004-637X/793/1/55
  * Kastner et al. (2015) Kastner, J. H., Qi, C., Gorti, U., et al. 2015, ApJ, 806, 75, doi: 10.1088/0004-637X/806/1/75
  * Krijt et al. (2020) Krijt, S., Bosman, A. D., Zhang, K., et al. 2020, arXiv e-prints, arXiv:2007.09517.  https://arxiv.org/abs/2007.09517
  * 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
  * Lee et al. (2015) Lee, S., Lee, J.-E., & Bergin, E. A. 2015, ApJS, 217, 30, doi: 10.1088/0067-0049/217/2/30
  * Lee et al. (2014) Lee, S., Lee, J.-E., Bergin, E. A., & Park, Y.-S. 2014, ApJS, 213, 33, doi: 10.1088/0067-0049/213/2/33
  * Liang et al. (2007) Liang, M.-C., Heays, A. N., Lewis, B. R., Gibson, S. T., & Yung, Y. L. 2007, ApJ, 664, L115, doi: 10.1086/520881
  * Loison et al. (2019) Loison, J.-C., Wakelam, V., Gratier, P., & Hickson, K. M. 2019, MNRAS, 484, 2747, doi: 10.1093/mnras/sty3293
  * London (1978) London, R. 1978, ApJ, 225, 405, doi: 10.1086/156502
  * Lucas & Liszt (1998) Lucas, R., & Liszt, H. 1998, A&A, 337, 246
  * Lucy (1999) Lucy, L. B. 1999, A&A, 344, 282
  * Magalhães et al. (2018) Magalhães, V. S., Hily-Blant, P., Faure, A., Hernand ez-Vera, M., & Lique, F. 2018, A&A, 615, A52, doi: 10.1051/0004-6361/201832622
  * Marty et al. (2011) Marty, B., Chaussidon, M., Wiens, R. C., Jurewicz, A. J. G., & Burnett, D. S. 2011, Science, 332, 1533, doi: 10.1126/science.1204656
  * Matar et al. (2010) Matar, E., Bergeron, H., Dulieu, F., et al. 2010, J. Chem. Phys., 133, 104507, doi: 10.1063/1.3484867
  * Menu et al. (2014) Menu, J., van Boekel, R., Henning, T., et al. 2014, A&A, 564, A93, doi: 10.1051/0004-6361/201322961
  * Min et al. (2009) Min, M., Dullemond, C. P., Dominik, C., de Koter, A., & Hovenier, J. W. 2009, A&A, 497, 155, doi: 10.1051/0004-6361/200811470
  * Min et al. (2016) Min, M., Rab, C., Woitke, P., Dominik, C., & Ménard, F. 2016, A&A, 585, A13, doi: 10.1051/0004-6361/201526048
  * Miotello et al. (2014) Miotello, A., Bruderer, S., & van Dishoeck, E. F. 2014, A&A, 572, A96, doi: 10.1051/0004-6361/201424712
  * Molano & Kamp (2012a) Molano, G. C., & Kamp, I. 2012a, A&A, 537, A138, doi: 10.1051/0004-6361/201015868
  * Molano & Kamp (2012b) —. 2012b, A&A, 547, A7, doi: 10.1051/0004-6361/201219943
  * Mumma & Charnley (2011) Mumma, M. J., & Charnley, S. B. 2011, ARA&A, 49, 471, doi: 10.1146/annurev-astro-081309-130811
  * Nomura & Millar (2005) Nomura, H., & Millar, T. J. 2005, A&A, 438, 923, doi: 10.1051/0004-6361:20052809
  * Nomura et al. (2016) Nomura, H., Tsukagoshi, T., Kawabe, R., et al. 2016, ApJ, 819, L7, doi: 10.3847/2041-8205/819/1/L7
  * Öberg et al. (2009) Öberg, K. I., Linnartz, H., Visser, R., & van Dishoeck, E. F. 2009, ApJ, 693, 1209, doi: 10.1088/0004-637X/693/2/1209
  * Offner et al. (2019) Offner, S. S. R., Gaches, B. A. L., & Holdship, J. R. 2019, ApJ, 883, 121, doi: 10.3847/1538-4357/ab3e02
  * Pinte et al. (2009) Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967, doi: 10.1051/0004-6361/200811555
  * Pinte et al. (2006) Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797, doi: 10.1051/0004-6361:20053275
  * Qi et al. (2008) Qi, C., Wilner, D. J., Aikawa, Y., Blake, G. A., & Hogerheijde, M. R. 2008, ApJ, 681, 1396, doi: 10.1086/588516
  * Redaelli et al. (2018) Redaelli, E., Bizzocchi, L., Caselli, P., et al. 2018, A&A, 617, A7, doi: 10.1051/0004-6361/201833065
  * Ritchey et al. (2015) Ritchey, A. M., Federman, S. R., & Lambert, D. L. 2015, ApJ, 804, L3, doi: 10.1088/2041-8205/804/1/L3
  * Robitaille (2010) Robitaille, T. P. 2010, A&A, 520, A70, doi: 10.1051/0004-6361/201015025
  * Rodgers & Charnley (2008) Rodgers, S. D., & Charnley, S. B. 2008, ApJ, 689, 1448, doi: 10.1086/592195
  * Roueff et al. (2015) Roueff, E., Loison, J. C., & Hickson, K. M. 2015, A&A, 576, A99, doi: 10.1051/0004-6361/201425113
  * Schwarz et al. (2018) Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2018, ApJ, 856, 85, doi: 10.3847/1538-4357/aaae08
  * Schwarz et al. (2019) —. 2019, ApJ, 877, 131, doi: 10.3847/1538-4357/ab1c5e
  * Semenov et al. (2010) Semenov, D., Hersant, F., Wakelam, V., et al. 2010, A&A, 522, A42, doi: 10.1051/0004-6361/201015149
  * Shinnaka et al. (2016) Shinnaka, Y., Kawakita, H., Jehin, E., et al. 2016, MNRAS, 462, S195, doi: 10.1093/mnras/stw2410
  * Sternberg et al. (2014) Sternberg, A., Le Petit, F., Roueff, E., & Le Bourlot, J. 2014, ApJ, 790, 10, doi: 10.1088/0004-637X/790/1/10
  * Taniguchi & Saito (2017) Taniguchi, K., & Saito, M. 2017, PASJ, 69, L7, doi: 10.1093/pasj/psx065
  * Terzieva & Herbst (2000) Terzieva, R., & Herbst, E. 2000, MNRAS, 317, 563, doi: 10.1046/j.1365-8711.2000.03618.x
  * Tielens (2008) Tielens, A. G. G. M. 2008, ARA&A, 46, 289, doi: 10.1146/annurev.astro.46.060407.145211
  * Tielens & Hollenbach (1985) Tielens, A. G. G. M., & Hollenbach, D. 1985, ApJ, 291, 722, doi: 10.1086/163111
  * van Boekel et al. (2017) van Boekel, R., Henning, T., Menu, J., et al. 2017, ApJ, 837, 132, doi: 10.3847/1538-4357/aa5d68
  * van der Tak et al. (2007) van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627, doi: 10.1051/0004-6361:20066820
  * van Zadelhoff et al. (2003) van Zadelhoff, G.-J., Aikawa, Y., Hogerheijde, M. R., & van Dishoeck, E. F. 2003, A&A, 397, 789, doi: 10.1051/0004-6361:20021592
  * Vasyunin et al. (2011) Vasyunin, A. I., Wiebe, D. S., Birnstiel, T., et al. 2011, ApJ, 727, 76, doi: 10.1088/0004-637X/727/2/76
  * Visser et al. (2018) Visser, R., Bruderer, S., Cazzoletti, P., et al. 2018, A&A, 615, A75, doi: 10.1051/0004-6361/201731898
  * Visser et al. (2009) Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323, doi: 10.1051/0004-6361/200912129
  * Walsh et al. (2010) Walsh, C., Millar, T. J., & Nomura, H. 2010, ApJ, 722, 1607, doi: 10.1088/0004-637X/722/2/1607
  * Walsh et al. (2014) Walsh, C., Millar, T. J., Nomura, H., et al. 2014, A&A, 563, A33, doi: 10.1051/0004-6361/201322446
  * Walsh et al. (2012) Walsh, C., Nomura, H., Millar, T. J., & Aikawa, Y. 2012, ApJ, 747, 114, doi: 10.1088/0004-637X/747/2/114
  * Wirström & Charnley (2018) Wirström, E. S., & Charnley, S. B. 2018, MNRAS, 474, 3720, doi: 10.1093/mnras/stx3030
  * Woitke (2015) Woitke, P. 2015, in European Physical Journal Web of Conferences, Vol. 102, European Physical Journal Web of Conferences, 00011
  * Woitke et al. (2016) Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103, doi: 10.1051/0004-6361/201526538
  * Zubko et al. (1996) Zubko, V. G., Mennella, V., Colangeli, L., & Bussoletti, E. 1996, MNRAS, 282, 1321, doi: 10.1093/mnras/282.4.1321

Table 1: Model parameters, and values for the reference model. Quantity | Symbol | Value
---|---|---
stellar mass | $M_{\star}$ | $0.74\,M_{\odot}$
effective temperature | $T_{\star}$ |  $4110\,$K
stellar luminosity | $L_{\star}$ | $0.28\,L_{\odot}$
UV luminosity | $L_{\rm UV}$ | $0.017\,L_{\odot}$
X-ray luminosity | $L_{X}$ | $1.6\times 10^{30}\,{\rm erg\,s^{-1}}$
strength of interstellar UV | $\chi^{\rm ISM}$ | 1 Draine field
cosmic ray H2 ionization rate | $\zeta_{\rm H_{2}}$ |  $\\!\\!1.0\times 10^{-19}$ s${}^{-1}\\!\\!\\!$
disk mass | $M_{\rm disk}$ | $0.01\,M_{\odot}$
gas/dust mass ratio | $\delta_{g}$ | 100
concentration factor of the large dust grain | $\chi_{d}$ | 0.2
inner disk radius | $R_{\rm in}$ | 0.07 AU
outer disk radius | $R_{\rm out}$ | 200 AU
reference and tapering-off radius | $R_{c}$ | 150 AU
column density power index | $\gamma$ | 1
reference scale height | $H_{\rm c}$ | 15 AU
flaring power index | $\beta$ | 1.3
Table 2: Initial abundances relative to the total hydrogen nuclei. Species | Abundanceaa z(y) means z$\times$10y. | Species | Abundanceaa z(y) means z$\times$10y.
---|---|---|---
H2 | $5.000(-01)$ | He | $1.400(-01)$
HCNbb The initial abundance of 15N bearing species is by a factor of
330 lower than that of 14N bearing species except for N2 (165).  | $1.000(-08)$ | NH3 icebb The initial abundance of 15N bearing species is by a factor of
330 lower than that of 14N bearing species except for N2 (165).  | $9.900(-06)$
Nbb The initial abundance of 15N bearing species is by a factor of
330 lower than that of 14N bearing species except for N2 (165).  | $5.100(-06)$ | N2 bb The initial abundance of 15N bearing species is by a factor of
330 lower than that of 14N bearing species except for N2 (165).  | $1.000(-06)$
C | $7.000(-07)$ | CH4 | $1.000(-07)$
CNbb The initial abundance of 15N bearing species is by a factor of
330 lower than that of 14N bearing species except for N2 (165).  | $6.600(-08)$ | H${}_{3}^{+}$ | $1.000(-08)$
HCO+ | $9.000(-09)$ | C2H | $8.000(-09)$
C+ | $1.000(-09)$ | Mg+ | $1.000(-11)$
Si+ | $1.000(-11)$ | S+ | $1.000(-11)$
Fe+ | $1.000(-11)$ |  |
Model | Species / Abundancesaa z(y) means z$\times$10y.
name | H2O ice | CO | C
I1 | $2.500(-04)$ | $1.000(-04)$ | $7.000(-07)$
I2 | $2.500(-05)$ | $1.000(-06)$ | $7.000(-07)$
I3 | $2.500(-07)$ | $1.000(-06)$ | $7.000(-07)$
I4 | $2.500(-06)$ | $1.000(-06)$ | $7.000(-06)$
Table 3: Free parameters in our models  Name |  $f$aaMass fraction of large dust grains. | Initial elemental abundancesbbz(y) means z$\times$10y.
---|---|---
|  | IMccModel names in Table 2. | C | O | C/O
MLI1 | 0.9 | I1 | 1.0(-4) | 3.5(-4) | 0.3
MLI2 | 0.9 | I2 | 1.9(-6) | 3.5(-6) | 0.5
MLI3 | 0.9 | I3 | 1.9(-6) | 1.3(-6) | 1.5
MLI4 | 0.9 | I4 | 1.8(-5) | 1.3(-5) | 1.5
MSI1 | 0.99 | I1 | 1.0(-4) | 3.5(-4) | 0.3
MSI2 | 0.99 | I2 | 1.9(-6) | 3.5(-6) | 0.5
MSI3 | 0.99 | I3 | 1.9(-6) | 1.3(-6) | 1.5
MSI4 | 0.99 | I4 | 1.8(-5) | 1.3(-5) | 1.5
MHI1 | Equation 6 | I1 | 1.0(-4) | 3.5(-4) | 0.3
MHI2 | Equation 6 | I2 | 1.9(-6) | 3.5(-6) | 0.5
MHI3 | Equation 6 | I3 | 1.9(-6) | 1.3(-6) | 1.5
MHI4 | Equation 6 | I4 | 1.8(-5) | 1.3(-5) | 1.5

Figure 1: Gas density and ionization rate (sum of X-ray and cosmic ray
contribution) in all models (see Section 2.1).

Figure 2: UV flux (in the unit of Draine field), dust temperature, gas
temperature, and N2 column density ($N_{\rm N_{2}}$) from the central star to
a given position (from top to bottom) for the three grain population models:
MLI3 ($f$=0.9; left), MSI3 ($f$=0.99; middle), and MHI3 (see Equation 6;
right). The white dotted lines in the bottom panels indicate the height where
the UV flux from the central star equals that from the upper atmosphere. Most
UV photons directly originated from the central star above the white dotted
line while the photons scattered by the dust grains are dominant in the disk
atmosphere below the white dotted line.  Figure 3:  Vertical cuts of Figures 1
and 2 for the MLI3 model at 20 AU (black), 40 AU (red), and 60 AU (blue). Top
right panel: N2 photodissociation rates with (solid) and without (dotted) the
self-shielding effect.

Figure 4:  Gas-phase abundance of N-bearing species for the MLI3 model. The
abundances relative to the total hydrogen nuclei in a log scale are plotted in
color image and contours.

Figure 5:  Nitrogen isotope ratio for the species presented in Figure 4 for
the MLI3 model. The contours indicate the abundances shown in Figure 4. The
regions where the abundance is lower than 10-10 are masked and plotted in
white. Note that the [14N/15N]elem ratio is 330 in our model.  Figure 6:
Vertical cuts of Figures 4 (top) and 5 (bottom) for the MLI3 model at 20 AU
(left), 40 AU (middle), and 60 AU (right). The middle rows plot the number
density normalized to the maximum number density of each species along the
vertical direction in order to show the contribution to the column density.
The dotted vertical lines indicate the height where the UV flux from the
central star equals that from the upper atmosphere. [14N/15N]elem (=330) is
presented in the dotted horizontal lines in the bottom panels.  Figure 7: Same
as Figure 6 except for the MSI3 model. Figure 8: HCN column density, column
density ratio of HC14N/HC15N, CN column density, column density ratio of
C14N/C15N, and C2H column density (from top to bottom) for the ML (left),
MS(center), and MH (right) models. The red, green, cyan, and purple lines
indicate the result of models with different initial abundances of I4, I3, I2,
and I1, respectively. The gray bars indicate the column densities derived from
the APEX 12 m telescope observations (Kastner et al., 2014). The column
densities and column density ratios of isotopologue derived from the ALMA
observations are presented in the black bars, filled circles, and gray hatched
bars (see text).  Figure 9: HCN column density, column density ratio of
HC14N/HC15N, CN column density, column density ratio of C14N/C15N, and C2H
column density (from top to bottom) for the test models. Reference model is
the MHI3 model (thick black line). The cyan, red, and green lines indicate the
results for the models without the vibrationally excited H2, the isotope-
exchange reactions, and the line-shielding effect, respectively. The blue
lines represent the results for the model in which the X-ray ionization rate
was a factor of 10 higher than that in the reference model. The observed data
are the same as those in Figure 8.

##  Appendix A Details of the PURE-C model

###  A.1 Dust grains

The dust size distribution, $f_{0}(a)$ $\rm[cm^{-4}]$, was assumed to follow a
power law of particle radius $a\,\rm[cm]$ as

$f_{0}(a)\;\propto\;a^{-3.5}\quad\mbox{with}\quad a\in[a_{\rm min},a_{\rm
max}]\ .$ (A1)

In this work, we used two grain populations: small dust grains with radii
$r_{g}$= 0.005 $\mu$m – 1 $\mu$m and large dust grains with radii $r_{g}$=
0.005 $\mu$m – 1 mm. The representative radii ($a_{\rm dust}$) of the small
and large dust grains are 0.1 $\mu$m and 10 $\mu$m following the method given
by Vasyunin et al. (2011).

The dust opacity was calculated using the OPACITY-TOOL222
https://dianaproject.wp.st-andrews.ac.uk/data-results-downloads/fortran-
package/ with the discrete dipole approximation (Min et al., 2016) and the
dust grains were considered to be a mixture of amorphous laboratory silicates
(Dorschner et al., 1995, $\rm Mg_{0.7}Fe_{0.3}SiO_{3}$) with amorphous carbon
(Zubko et al., 1996, BE-sample). We adopted the reference parameters for the
dust grains as given in Woitke et al. (2016): dust material density of 2.09 g
cm-3, volume fractions of silicate, amorphous carbon, and porosity of 60 %, 15
%, and 25 % respectively, and the maximum hollow volume ratio of 0.8.

###  A.2 Dust continuum radiative transfer

The dust continuum radiative transfer problem was solved using a Monte Carlo
method described in the article by Lucy (1999). The dust temperature was
calculated by assuming a radiative equilibrium:

$4\pi\int_{0}^{\infty}\kappa_{\nu}B_{\nu}(T)d\nu=4\pi\int_{0}^{\infty}\kappa_{\nu}J_{\nu}d\nu,$
(A2)

where $\kappa_{\nu}$ is the absorption coefficient of the dust grain at the
frequency $\nu$, $B_{\nu}(T)$ is the Planck function with the dust temperature
T, and $J_{\nu}$ is the mean intensity in the cell. For a given dust grain,
the left-hand side of Equation A2 is the total emission from the dust grain in
the cell with the dust temperature. The right-hand side is the total energy
absorbed by the dust grain and is derived numerically by

$\displaystyle 4\pi\int_{0}^{\infty}{\kappa_{\nu}J_{\nu}d\nu}=\frac{L_{\rm
tot}/N_{\gamma}}{V}\sum_{i}\Delta\tau_{i}$ (A3)

where $L_{\rm tot}$ is the total luminosity in the system, $N_{\gamma}$ is the
number of model photons, $V$ is the volume of the cell, and $\Delta\tau_{i}$
is the optical depth of the trajectory of a model photon $i$ within the grid
cell. The dust temperature is converged within a few iterations.

In the very optically thick midplane of the disk, where the Monte Carlo method
was not efficient, a modified random walk method (Min et al., 2009;
Robitaille, 2010) was applied. When a photon escapes from the optically thick
and homogeneous sphere of radius R and dust grain density $\rho_{\rm dust}$,
the dust grain within the sphere absorbs an energy per unit dust grain density
of

$E=(L_{\rm tot}/N_{\gamma})\kappa_{\rm
P}\left(\frac{R^{2}}{2}\overline{\chi^{-1}}\rho_{\rm dust}\right)$ (A4)

where

$\kappa_{\rm
P}=\int_{0}^{\infty}\kappa_{\nu}B_{\nu}(T)d\nu/\int_{0}^{\infty}B_{\nu}(T)d\nu$
(A5)
$\overline{\chi^{-1}}=\int_{0}^{\infty}\chi_{\nu}^{-1}B_{\nu}(T)d\nu/\int_{0}^{\infty}B_{\nu}(T)d\nu.$
(A6)

Here, the dust extinction coefficient $\chi_{\nu}$ is defined as
$\kappa_{\nu}+(1-g_{\nu})\sigma_{\nu}$ with an asymmetry parameter $g_{\nu}$
and the scattering coefficient $\sigma_{\nu}$ (Ishimaru, 1978; Min et al.,
2009). It should be noted that $B_{\nu}(T)$ in both Equations A5 and A6 is
replaced by $\partial B_{\nu}(T)/\partial T$ when the Bjorkman-Wood method
(Bjorkman & Wood, 2001) is used.

We calculated the dust temperature for the protoplanetary disks in the
benchmark test (Pinte et al., 2009)333https://ipag.osug.fr/~pintec/benchmark/.
Two cases were tested for isotropic scattering with disk masses of $3\times
10^{-8}$ $M_{\odot}$ ($\tau=10^{3}$) and $3\times 10^{-5}$ $M_{\odot}$
($\tau=10^{6}$). Figure 10 shows the dust temperature distribution at the
midplane of the disk corresponding to our code (red) and for the MCFOST (blue,
Pinte et al., 2006). In both cases, the agreement is very good and our results
are within the deviations of other codes participating in the benchmark test.

###  A.3 X-ray ionization rate

X-ray photons are an important source of ionization at the intermediate column
densities (Igea & Glassgold, 1999; Offner et al., 2019). In this work, we
adopted the analytic formula given in Igea & Glassgold (1999). The ionization
rate can be expressed as a function of the column density of hydrogen nuclei
vertically above ($N_{\rm a}$) and below ($N_{\rm b}$) at a given location in
the disk (Igea & Glassgold, 1999; Bai & Goodman, 2009):

$\displaystyle\zeta_{X}$ $\displaystyle=$ $\displaystyle
L_{X,29}\left(\frac{r}{1\,{\rm AU}}\right)^{-2.2}\times$ (A7)
$\displaystyle(\zeta_{1}[e^{-(N_{\rm a}/N_{\rm 1})^{\alpha_{1}}}+e^{-(N_{\rm
b}/N_{\rm 1})^{\alpha_{1}}}]+$ $\displaystyle\zeta_{2}[e^{-(N_{\rm a}/N_{\rm
2})^{\alpha_{2}}}+e^{-(N_{\rm b}/N_{\rm 2})^{\alpha_{2}}}]),$

where $L_{X,29}\equiv L_{X}/10^{29}$ erg s-1, $L_{\rm X}$ is the X-ray
luminosity, $\zeta_{1}=6\times 10^{-12}$ s-1, $\zeta_{2}=10^{-15}$ s-1,
$\alpha_{1}=0.4$, $\alpha_{2}=0.65$, $N_{1}=1.5\times 10^{21}$ cm-2, and
$N_{2}=7\times 10^{23}$ cm-2.

###  A.4 Chemistry

Abundances of species were derived from the equations of chemical kinetics
describing the formation and destruction of the species:

$\frac{dn_{i}}{dt}=\sum_{l,m}k_{lm}n_{l}n_{m}-n_{i}\sum_{l}k_{il}n_{l}+\sum_{l}k_{l}n_{l}-n_{i}\sum_{l}k_{l}+k_{i}^{\rm
des}n_{i}^{s}-k_{i}^{\rm ads}n_{i}$ (A8)
$\frac{dn_{i}^{s}}{dt}=\sum_{l,m}k^{s}_{lm}n_{l}^{s}n_{m}^{s}-n_{i}^{s}\sum_{l}k_{il}^{s}n_{l}^{s}-k_{i}^{\rm
des}n_{i}^{s}+k_{i}^{\rm ads}n_{i}$ (A9)

where $n_{i}$ and $n_{i}^{s}$ are the $i$-th species (cm-3) in the gas and ice
phases, respectively, $k_{lm}$ ($k_{il}$) and $k_{i}$ ($k_{l}$) are the gas-
phase reaction rates (in units of cm3 s-1 and s-1, respectively), $k_{i}^{\rm
ad}$ and $k_{i}^{\rm des}$ denote the adsorption and desorption rates (s-1),
respectively, and $k_{lm}^{s}$ ($k_{il}^{s}$) is the surface reaction rate
(cm3 s-1). All reactions except the surface reaction are described in Lee et
al. (2014).

The chemical code of Lee et al. (2014) was developed from the Heidelberg
“ALCHEMIC” code (Semenov et al., 2010), which is based on the publicly
available DVODPK (Differential Variable-coefficient Ordinary Differential
equation solver with the Preconditioned Krylov method GMRES for the solution
of linear systems) ODE package444http://www.netlib.org/ode/vodpk.f and a high-
performance sparse asymmetric MA28 solver from the Harwell Mathematical
Software Library 555http://www.hsl.rl.ac.uk/. We replaced the MA28 solver with
an updated solver version, MA48, which is faster than MA28 by a factor of two.

We updated some reactions by considering the recent studies. When neutral
species are adsorbed on the grain surface, we assume a sticking coefficient of
1.0, except for H and H2. For H and H2, we follow the gas-temperature-
dependent formulae derived from laboratory experiments (Matar et al., 2010;
Chaabouni et al., 2012). The wavelength-dependent photodesorption experiments
have shown that photodesorption is induced by photoabsorption in the dominant
ice component and energy transfer to the surface molecules (Fayolle et al.,
2011, 2013; Bertin et al., 2012, 2013). Therefore, we used a single
photodesorption rate of H2O ($Y_{\rm ph}=10^{-3}$ molecule photon-1, Öberg et
al., 2009) for all species.

Changes in the dust opacity in the UV range due to the coagulation of dust
grains was considered in the UV radiative transfer process. However, this also
affects the UV photons induced by cosmic rays (Molano & Kamp, 2012a, b, CRP).
Although the opacities of some molecules can be comparable to the dust opacity
(Molano & Kamp, 2012b), we simply enhanced the photodissociation rates by CRP
in the network using a scale factor of $\sigma_{\rm UV}/2\times 10^{-21}$
cm-2, where $\sigma_{\rm UV}$ is the UV extinction cross section for the dust
grain in the model (Molano & Kamp, 2012a). In addition, we also calculated the
number of CRPs ($N_{\rm CU}$) following the expression given by Molano & Kamp
(2012a):

$\begin{split}N_{\rm
CU}=12500\times\left(\frac{1}{1-\omega}\right)\left(\frac{\zeta_{\rm
H_{2}}}{5\times 10^{-17}\,{\rm s}^{-1}}\right)\\\ \times\left(\frac{2\times
10^{-21}\,{\rm cm}^{-2}}{\sigma_{\rm UV}}\right)\left(\frac{n_{\rm
H_{2}}/n_{\rm gas}}{0.5}\right),\end{split}$ (A10)

where $\omega$ is the grain albedo and $n_{\rm H_{2}}$ and $n_{\rm gas}$ are
the number densities of molecular hydrogen and hydrogen nuclei, respectively.
It should be noted here that this approximation is available when the gas
opacity is much lower than the dust opacity (see Molano & Kamp, 2012b).

We also included the reactions with the vibrationally excited H2 (v-H2) in our
model. Following Tielens & Hollenbach (1985) and Bruderer et al. (2012), the
v-H2 is a vibrationally excited pseudo-level with an energy of 30163 K
(London, 1978). When two-body reactions with H2 had an activation barrier, the
exponential factor $\gamma$ was replaced with max(0,$\gamma$\- 30163 K) in the
reactions with v-H2. The UV pumping rate was 8 times the H2 photodissociation
rate (Sternberg et al., 2014), and the spontaneous decay rate was 2$\times
10^{-7}$ s-1 (London, 1978). The collision rates of H and H2 were adopted from
(Tielens & Hollenbach, 1985).

For the grain surface chemistry of dust grains having a size of
$a_{\text{dust}}$ and number density of $n_{\text{dust}}$, we considered only
the Langmuir–Hinshelwood mechanism:

$\begin{split}k_{lm}^{s}=&\kappa_{lm}\Bigg{(}k_{\text{hop}}^{s}(l)+k_{\text{hop}}^{s}(m)\Bigg{)}\\\
&\times\frac{1}{N_{\text{site}}4\pi
a_{\text{dust}}^{2}n_{\text{dust}}}~{}[\text{cm}^{3}\text{s}^{-1}]\end{split}$

where $\kappa_{lm}$ is the probability that the reaction occurs. The thermal
hopping rate are given by

$k_{\text{hop}}^{s}(l)=\sqrt{\frac{2N_{\text{site}}k_{\text{B}}E_{\text{b}}(l)}{\pi^{2}m_{l}}}\,\exp\Bigg{(}-\frac{E_{\text{diff}}(l)}{T_{\text{dust}}}\Bigg{)},$
(A11)

where $N_{\text{site}}$ is the number density of the surface site ($\simeq
1.5\times 10^{15}\,\text{cm}^{-2}$), $k_{\text{B}}$ is the Boltzmann constant,
$E_{\text{b}}(l)$ and $E_{\text{diff}}(l)$ are the binding and diffusion
energies (in K) of the species $l$, respectively, and $m_{l}$ is its mass. We
assumed that $E_{\text{diff}}(l)=0.5E_{\text{b}}(l)$ (Garrod & Pauly, 2011).

When the reaction was exothermic and barrierless, we set $\kappa_{lm}=1$
(Hasegawa et al., 1992). For exothermic reactions with an activation barrier,
denoted by $E_{\text{A},lm}$, we calculated $\kappa_{lm}$ as the result of the
competition among reactions, hopping, and evaporation, as suggested by Garrod
& Pauly (2011):

$\kappa_{lm}=\frac{\nu_{lm}\kappa_{lm}^{*}}{\nu_{lm}\kappa_{lm}^{*}+k_{\text{hop}}^{s}(l)+k_{\text{hop}}^{s}(m)}.$

where $\kappa_{lm}^{*}$ can be expressed as
$\exp(-E_{\text{A},lm}/T_{\text{dust}})$ or the quantum mechanical probability
for tunneling through a rectangular barrier of thickness $a$:
$\kappa_{lm}^{*}=\exp[-2(a/\hbar)(2\mu E_{\text{A},lm})^{1/2}]$, where $\mu$
is the reduced mass (see Hasegawa et al., 1992, for details). Here, $\nu_{lm}$
was taken to be equal to the larger value of the characteristic frequencies of
the two reactants $l$ and $m$ (Garrod, 2013).

###  A.5 Self-Shielding and Line-Shielding

Far ultraviolet (FUV) radiation is emitted from the central star and the
interstellar radiation field (ISRF). The FUV radiative transfer was calculated
using a Monte Carlo method considering anisotropic scattering (see Lee et al.,
2014), and the FUV fluxes from the two sources were treated separately. We
used the opacities of the dust grains as given in Section A.1. The
contributions of the UV flux along the radial, $G_{\rm dust}^{in}(ir,iz)$, and
vertical direction, $G_{\rm dust}^{up}(ir,iz)$, were calculated for a given
grid cell, $(ir,iz)$, in the disk. Furthermore, a local column density,
$N_{\rm local}(ir,iz)$, within a grid cell, $(ir,iz)$, was calculated as
follows:

$N_{\rm local}(ir,iz)=n_{\rm gas}(ir,iz)\frac{\sum_{i}I_{i}(s)\Delta
s}{\sum_{i}I_{i}(s)},$ (A12)

where $I_{i}(s)$ is the dust-attenuated intensity of a model photon $i$
penetrating the grid cell $(ir,iz)$, and $\Delta s$ is the path length
traveled within the grid cell.

When taking the self-shielding effects of H2, CO, C, and N2 into account,
their column densities along the path of UV photons should also be calculated.
Thus, prior knowledge of their abundances in all grid cells is required. This
problem could be solved using an iterative method (van Zadelhoff et al., 2003;
Bruderer et al., 2012). In the 2D disk model, the UV photons from the central
star mainly enter the grid cells through their inner as well as upper
boundaries. Therefore, we solved the problem in the consecutive order from the
upper-innermost grid to the lower-outermost grid.

Then, the column densities of species X along the photon path from the central
star to a given grid cell $(ir,iz)$ were derived as follows:

$N_{\rm X}(ir,iz)=X_{\rm X}(ir,iz)\cdot N_{\rm local}(ir,iz)+N_{\rm X}^{\rm
ext}(ir,iz),$ (A13)

where the abundance of the species X is $X_{\rm X}(ir,iz)$, the local column
density is $N_{\rm local}(ir,iz)$ (see Equation A12), and the external column
density of species X is

$N_{\rm X}^{\rm ext}(ir,iz)=\left\\{\begin{array}[]{cl}N_{\rm
X}(ir-1,iz),&\;if\,\,G_{\rm dust}^{in}(ir,iz)/G_{\rm dust}^{up}(ir,iz)>10,\\\
\frac{N_{\rm X}(ir-1,iz)\cdot G_{\rm dust}^{in}(ir,iz)+N_{\rm X}(ir,iz-1)\cdot
G_{\rm dust}^{up}(ir,iz)}{G_{\rm dust}^{in}(ir,iz)+G_{\rm
dust}^{up}(ir,iz)},&\;if\,\,0.1\leq G_{\rm dust}^{in}(ir,iz)/G_{\rm
dust}^{up}(ir,iz)\leq 10,\\\ N_{\rm X}(ir,iz-1),&\;if\,\,G_{\rm
dust}^{in}(ir,iz)/G_{\rm dust}^{up}(ir,iz)<0.1,\end{array}\right.$ (A14)

Here, we used the average value when the contributions of the vertical and
radial directions to the dust-attenuated UV flux were comparable because the
extreme cases could lead to incorrect values of the column density. For
example, near the midplane in the outer disk, $G_{\rm dust}^{in}(ir,iz)$ is
lower than $G_{\rm dust}^{up}(ir,iz)$, whereas $N_{\rm X}(ir-1,iz)$ is much
higher than $N_{\rm X}(ir,iz-1)$, which leads to an unreasonably high value of
the column density. Furthermore, $N_{\rm X}(ir,iz+1)$ and $N_{\rm X}(ir-1,iz)$
are comparable, as shown in the bottom panels of Figure 2, when $G_{\rm
dust}^{up}(ir,iz)$ and $G_{\rm dust}^{in}(ir,iz)$ are similar. Therefore, the
photodissociation rate of N2 with the self-shielding effect in our calculation
is consistent with the sum of those in the vertical and radial directions even
though the self-shielding function of N2 is nonlinear.

The left panel in Figure 11 shows the total hydrogen column densities from the
central star to a given position along the vertical direction obtained using
different methods at the radius of 40 AU. The black line ($N_{\rm MC}$) is the
column density along the path of the photons from the central star to the grid
$(ir,iz)$ obtained using the Monte Carlo simulation:

${\rm ln}\,N_{\rm MC}(ir,iz)=\frac{\sum_{i}I_{i}(s){\rm
ln}\,N_{i}(s)}{\sum_{i}I_{i}(s)},$ (A15)

where $N_{i}(s)$ is the column density along the path of the model photon $i$
from the central star to the grid $(ir,iz)$. In this method, ${\rm
ln}\,N_{i}(s)$ is averaged instead of $N_{i}(s)$ because the dynamical range
of $N_{i}(s)$ is significantly large depending on the photon paths, owing to
which it is difficult to obtain the averaged $N_{i}(s)$ numerically. Our
method using Equation A13 (the red line in the left panel in Figure 11)
derives the column density between the column densities integrated radially
from the central star (green) and vertically from the upper atmosphere (blue)
and is consistent with $N_{\rm MC}$ (black) within a factor of two. Therefore,
our method was able to trace a realistic column density.

The N2 photodissociation rates obtained using the different methods are shown
in the right panel of Figure 11. The red and blue lines indicate the N2
photodissociation rates obtained using the column densities shown with the
same color in the left panel. The black line indicates the exact solution
because the self-shielding effect is calculated along the path of the UV
photons using the Monte Carlo method when the abundance distribution of atomic
H, H2, and N2 for the MLI3 model is adopted. Our method reproduces the same
photodissociation rate in the atmosphere, whereas it slightly underestimates
the photodissociation rate at a lower height where the N2 self-shielding
effect is dominant as compared to the dust attenuation of UV photons.

The line-shielding effect for photodissociation by the stellar UV photons was
also considered in our model. Heays et al. (2017)
666https://home.strw.leidenuniv.nl/~ewine/photo/ provides tables for the line-
shielding effect due to H, H2, C, CO, and N2. They assumed an excitation
temperature of 100 K and a Gaussian Doppler broadening width of b = 3 km s-1.
The combined line-shielding functions were assumed to be a product of the
line-shielding functions corresponding to the abovementioned five species. The
column density was calculated using Equation A13.

###  A.6 Gas energetics

In the disk, the properties of PAH and dust grains are different from those in
the ISM. For the photoelectric heating of large dust grains, we adopted the
approximate formula for silicate grains in Kamp & van Zadelhoff (2001),

$\Gamma_{\rm pe}=2.5\times 10^{-4}\,\sigma_{\rm UV}^{\rm
tot}\,\epsilon\,G_{\rm dust}$ (A16)

with

$\epsilon=\frac{6\times 10^{-2}}{1+1.8\times
10^{-3}x^{0.91}}+\frac{y\left(10^{-4}\,T\right)^{1.2}}{1+10^{-2}x},$ (A17)

and

$y=\left\\{\begin{array}[]{cl}0.7,&\;if\,\,x\leq 10^{-4}\\\
0.36,&\;if\,\,10^{-4}\leq x\leq 1\\\ 0.15,&\;if\,\,x>1.\end{array}\right.$
(A18)

Here, $G_{\rm dust}$ is the dust-attenuated UV strength in the Habing field,
$\sigma_{\rm UV}^{\rm tot}$ is the total FUV absorption coefficient per unit
volume, and $x\equiv G_{\rm dust}\sqrt{T_{\rm gas}}/n_{\rm e}$ is the grain
charge parameter, where $n_{\rm e}$ is the electron number density. It should
be noted here that a new laboratory experiment has shown that the
photoelectric yields of large dust grains are much larger than the
theoretically expected values (Abbas et al., 2006; Woitke, 2015).

The abundances of PAHs in the disk are much smaller (order of 1–2) than those
in the ISM (Geers et al., 2006). Thus, we scaled down the photoelectric
heating rate by the PAH abundance ratio with respect to the standard PAH
abundance in the ISM (Tielens, 2008), $f_{\rm PAH}$ and adopted its value to
be 0.01.

The dust thermal accommodation was calculated, taking into account the dust
grain $i$ of size $a^{i}_{\rm dust}$, number density, $n^{i}_{\rm dust}$, and
dust temperature, $T^{i}_{\rm dust}$ (Burke & Hollenbach, 1983; Groenewegen,
1994; Woitke, 2015), as follows:

$\displaystyle\Gamma_{\rm dust}-\Lambda_{\rm dust}$ $\displaystyle=$
$\displaystyle\sum_{i}\pi(a^{i}_{\rm dust})^{2}n^{i}_{\rm
dust}\sqrt{\frac{8kT_{\rm gas}}{\pi m_{\rm H}}}n_{\rm gas}\alpha_{\rm acc}$
(A19) $\displaystyle\times\left(2kT^{i}_{\rm dust}-2kT_{\rm gas}\right)$
$\displaystyle\approx$ $\displaystyle 4\times 10^{-12}n_{\rm
gas}\sum_{i}\pi(a^{i}_{\rm dust})^{2}n^{i}_{\rm dust}\alpha_{\rm acc}$
$\displaystyle\times\sqrt{T_{\rm gas}}(T^{i}_{\rm dust}-T_{\rm gas}),$

with the efficiency for inelastic collision given by

$\alpha_{\rm acc}\approx 0.1+0.35\,\exp\left(-\sqrt{\frac{T_{\rm
gas}+T^{i}_{\rm dust}}{500K}}\right).$ (A20)

Figure 10: Dust temperature distribution at the mid-plane in the optically
thin case (left; $\tau=10^{3}$) and in the optically thick case (right;
$\tau=10^{6}$). The red and blue lines indicate the results from our model and
MCFOST (Pinte et al., 2006), respectively.

Figure 11:  Total hydrogen column density from the central star ($N({\rm
H_{tot}})$) along the vertical direction (left) and the photodissociation rate
of N2 (right) at the radius of 40 AU in the ML model. Left: the black and red
lines indicate the $N({\rm H_{tot}})$ calculated by the Monte Carlo method
(Equation A15) and with Equation A13, respectively. The column densities
integrated radially from the central star and vertically from the upper
atmosphere are presented in the green and blue lines, respectively. The black
dotted line represents the contribution of UV photons along the radial
direction to the UV flux for a given height ($G_{\rm dust}^{in}$). Right: the
red and blue lines indicate the photodissociation rate of N2 with the self-
shielding effect using the column densities in the left panel with the same
colors. The black line represents the photodissociation of N2 with the self-
shielding effect, which is calculated along the path of model photons adopting
the abundance distributions of atomic H, H2, and N2 by the Monte Carlo method.