# Far-ultraviolet and X-ray irradiated protoplanetary disks: a grid of models. I. The disk structure

R. Meijerink, G. Aresu, I. Kamp, M. Spaans, W. -F. Thi, P. Woitke

► To cite this version:

R. Meijerink, G. Aresu, I. Kamp, M. Spaans, W. -F. Thi, et al. Far-ultraviolet and X-ray irradiated protoplanetary disks: a grid of models. I. The disk structure. Astronomy and Astrophysics - A&A, 2012, 547, 10.1051/0004-6361/201219863. insu-03612398

HAL Id: insu-03612398

https://insu.hal.science/insu-03612398

Submitted on 17 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Far-ultraviolet and X-ray irradiated protoplanetary disks: a grid of models

## I. The disk structure*

R. Meijerink<fnref n="1" />,<fnref n="2" />, G. Aresu<fnref n="1" />, I. Kamp<fnref n="1" />, M. Spaans<fnref n="1" />, W.-F. Thi<fnref n="3" />, and P. Woitke<fnref n="4" />

Received 22 June 2012 / Accepted 1 August 2012

### ABSTRACT

**Context.** Planets are thought to eventually form from the mostly gaseous (~99% of the mass) disks around young stars. The density structure and chemical composition of protoplanetary disks are affected by the incident radiation field at optical, far-ultraviolet (FUV), and X-ray wavelengths, as well as by the dust properties.

**Aims.** The effect of FUV and X-rays on the disk structure and the gas chemical composition are investigated. This work forms the basis of a second paper, which discusses the impact on diagnostic lines of, e.g., C+, O, H₂O, and Ne+ observed with facilities such as *Spitzer* and *Herschel*.

**Methods.** A grid of 240 models is computed in which the X-ray and FUV luminosity, minimum grain size, dust size distribution, and surface density distribution are varied in a systematic way. The hydrostatic structure and the thermo-chemical structure are calculated using Protoplanetary Disk Model (ProDiMo), with the recent addition of X-rays.

**Results.** The abundance structure of neutral oxygen is very stable to changes in the X-ray and FUV luminosity, and the emission lines will thus be useful tracers of the disk mass and temperature. The C+ abundance distribution is sensitive to both X-rays and FUV. The radial column density profile shows two peaks, one at the inner rim and a second one at a radius $r = 5–10$ AU. The minimum is caused by shadowing from the inner rim. The fluctuations in value of the column density as a function of radius are smoothed out when FUV and X-ray luminosities increase. Ne+ and other heavy elements with an ionization potential higher than *IP* > 13.6 eV have a very strong response to X-rays, and the column density in the inner disk increases by two orders of magnitude from the lowest ($L_X = 10^{29}$ erg s⁻<fnref n="1" />) to the highest considered X-ray flux ($L_X = 10^{32}$ erg s⁻<fnref n="1" />). FUV confines the Ne+ ionized region to areas closer to the star at low X-ray luminosities ($L_X = 10^{29}$ erg s⁻<fnref n="1" />). This is indirectly caused by changes in the disk structure. The radial column densities of Ne+ are higher than $10^{12}$ cm⁻<fnref n="1" /> out to radii $r > 100$ AU (at $L_{FUV} \ge 10^{31}$ erg s⁻<fnref n="1" />), whereas the column density already drops below this value at radii $r > 20$ AU at $L_{FUV} = 10^{32}$ erg s⁻<fnref n="1" />. H₂O abundances are enhanced by X-rays due to higher temperatures in the inner disk than in the FUV only case, thus leading to a more efficient neutral-neutral formation channel. Also, the higher ionization fraction provides an ion-molecule route in the outer disk. The line fluxes and profiles are affected by the effects on these species, thus providing diagnostic value in the study of FUV and X-ray irradiated disks around T Tauri stars.

**Key words.** protoplanetary disks – X-rays: ISM – ultraviolet: ISM

## 1. Introduction

New observing facilities in the past decade pushed our understanding of protoplanetary disks from a rough picture of a vertically layered structure to a wealth of details on the composition and two-dimensional structure of the gas and dust of these disks. In the infrared, the *Spitzer* Space Telescope performed systematic studies of nearby star-forming regions. The spectral energy distributions (SEDs) revealed the physical structure of disks, e.g., the presence of gaps, source-to-source variations and important statistics on SED types, which allows evolutionary scenarios to be built (cf. Merín et al. 2010). Extensive gas-phase emission line studies with *Spitzer* provide a first indication of chemical diversity across stellar spectral types (cf. Najita et al. 2009;

Pontoppidan et al. 2010; Lahuis et al. 2007). In the near-infrared, there have been large ground-based studies with, e.g., the VLT and the Keck telescope resolving the lines to study the kinematics (spatial origin of the lines in hot disk surfaces and winds) and their excitation mechanisms (thermal and fluorescence; Brittain et al. 2010; Pontoppidan et al. 2011; Fedele et al. 2011). In the past two years, the *Herschel* Space Observatory extended the spectral window to the far-IR and submm, with spectral line scans within the Dust, Ice, and Gas In Time (DIGIT) key program (PI: N. Evans) for many Herbig disks, e.g., HD 100546 (Sturm et al. 2010), water studies of a few selected targets, such as TW Hya (Hogerheijde et al. 2011) from the Water in Star-forming Regions with *Herschel* (WISH) key program (PI: E. van Dishoeck), and large statistical gas surveys targeting the dominant cooling lines (e.g., Mathews et al. 2010; Meeus et al. 2012; Dent et al., in prep.; Riviere-Marichalar et al. 2012) in

<footnote n="1">¹ Kapteyn Astronomical Institute, PO Box 800, 9700 AV Groningen, The Netherlands
e-mail: meijerink@astro.rug.nl</footnote>

<footnote n="2">² Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, Netherlands</footnote>

<footnote n="3">³ UJF-Grenoble 1/CNRS-INSU, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG) UMR5274, 38041 Grenoble, France</footnote>

<footnote n="4">⁴ SUPA, School of Physics & Astronomy, University of St. Andrews, North Haugh, St. Andrews KY16 9SS, UK</footnote>

<footnote>* Appendix B is available in electronic form at
http://www.aanda.org</footnote>
the GAS in Protoplanetary Systems (GASPS) key program (PI: B. Dent). The HIFI and some PACS line detections deal with the cold water, while the 63.3 µm H₂O line possibly probes the inner water reservoir.

Observational studies report trends in emission line strength with the irradiation of the central star. Güdel et al. (2010) find that the [Ne II] emission at 12.81 µm correlates with the X-ray luminosity of the central star; the slope of the correlation for nonjet sources is 0.44 ± 0.18. Pinte et al. (2010) and Meeus et al. (2012) show that the [O I] 63 µm line flux increases with stellar luminosity and that this is largely driven by the strength of the far-ultraviolet (FUV) luminosity. Riviere-Marichalar et al. (2012) report a tentative trend of the 63.3 µm water line flux with X-ray luminosity.

It is obvious from previous works that the most relevant parameter for the SED is the total stellar luminosity; in most cases, $L_X$ and $L_{FUV}$ present only a small fraction of this and hence by themselves do not cause significant SED changes. However, recent near-IR high-contrast imaging with HiCIAO on the Subaru 8.2 m telescope allowed probing of the inner disk structure on scales below 0.1" (Thalmann et al. 2010; Hashimoto et al. 2011). This offers the possibility of a direct probe, e.g., the height of the inner rim as predicted by protoplanetary disk models. Thi et al. (2011) predicted that the height of the inner rim could be a factor two higher when the vertical hydrostatic structure takes into account the gas temperature at the inner rim. Aresu et al. (2011) show that the height of the inner rim increases in disk models with the impinging stellar $L_X$.

In the last decade, much theoretical progress has taken place to support the interpretation of the wealth of new observations, specifically the gas observations. The foundation was laid by Chiang & Goldreich (1997), D'Alessio et al. (1998), and Dullemond et al. (2001), whose layered dust disk models have helped to interpret SEDs for irradiated disks (see Dullemond et al. 2007, for a review of disk structure models). Based on this, several groups focussed on the effects of FUV and X-rays on the thermal and chemical structure of the gas in a pre-prescribed protoplanetary disk model, e.g., van Zadelhoff et al. (2003), Kamp & Dullemond (2004), Jonkheid et al. (2004), and Glassgold et al. (2004). In a subsequent step, Nomura & Millar (2005), Gorti & Hollenbach (2008), and Woitke et al. (2009) studied the chemical structure of disks, while self-consistently solving for the vertical hydrostatic equilibrium using the gas temperature. Most recently, Aresu et al. (2011) performed an exploratory study to assess the relative importance of FUV and X-rays by expanding the existing Protoplanetary Disk Model (ProDiMo) code to include X-ray processes. The modeling efforts that solve for the vertical disk structure are computationally expensive because the chemical networks, heating/cooling balance, 2D continuum radiative transfer, and hydrostatic equilibrium have to be solved iteratively. Hence, these studies have largely focussed on a single representative disk model, or at most a handful of models, varying one specific parameter. However, large-model grids are in principle required to understand the potential diagnostic power of certain observables. So far examples for this are large grids (of the order of 200 000–300 000 models) of parametrized dust disk models to study SED diagnostics (Robitaille et al. 2006) and of parametrized gas/dust disks to study the gas emission line diagnostics (Woitke et al. 2011; Kamp et al. 2011).

In this paper, we perform for the first time an extensive analysis of a grid of 240 self-consistent disk models (including the vertical disk structure) to study the effects of X-rays, FUV, and the relative importance of grain size and gas surface density distribution on the thermal, chemical, and physical structure of

disks around T Tauri stars. This paper builds on the implementation of X-rays into the ProDiMo code as described in Aresu et al. (2011). While we discuss here foremost the physical and chemical structure of the models and how they change with irradiation, a companion paper (Aresu et al. 2012, Paper II), discusses the power of line diagnostic in disentangling these effects. The paper is structured in the following way: Sect. 2 describes the updates on the code and the range of parameters used in the grid. The effects on disk temperature and density structure will be discussed in Sect. 3. The resulting distribution of various key species and its key reactions are extensively described in Sect. 4. Section 5 shows the radial column density distributions. Those are key in understanding the line profiles, which is the topic of the accompanying Paper II. Section 6 summarizes the conclusions and implications of the paper.

## 2. Updates on ProDiMo and the calculated grid

The original ProDiMo code (Woitke et al. 2009) is based on the code as discussed by Kamp & Dullemond (2004). The original model includes (1) frequency-dependent two-dimensional dust continuum radiative transfer, (2) kinetic gas-phase and FUV photo-chemistry, (3) ice formation, and (4) detailed non-local thermodynamic equilibrium (non-LTE) heating and cooling with (5) a consistent calculation of the hydrostatic disk structure. The models are characterized by a high degree of consistency among the various physical, chemical, and radiative processes, since the mutual feedbacks are solved iteratively. Aresu et al. (2011) included X-ray heating and chemical processes, guided by work from Maloney et al. (1996), Glassgold et al. (2004), and Meijerink & Spaans (2005), while most recently the X-ray chemistry processes are updated following Ádámkovics et al. (2011). This includes an extension of our chemical network with species such as Ne, Ar, and their singly and doubly ionized states, as well as other heavy elements. We added $Ne^+$, $Ne^{2+}$, $Ar^+$, and $Ar^{2+}$ in the heating/cooling balance and also implemented an extended sulfur chemistry (following Leen & Graff 1988; Meijerink et al. 2008, e.g., SO, $SO_2$, HS, and their related reactions were added) to achieve a proper calculation of the sulfur-based species abundances. This allowed us to make more realistic predictions of, e.g., the sulfur fine-structure lines. A detailed description on the treatment of X-ray chemistry in the code is given in Appendix A.

Our input spectrum is composed of the stellar spectrum and a FUV excess as described in Woitke et al. (2009), along with an X-ray component as outlined in Aresu et al. (2011). The FUV excess (from 92.5 to 205 nm) spectrum is a power law with $I_v(\lambda) \propto \lambda^{1.2}$. The X-ray spectrum (0.1–10 keV) has been chosen following Glassgold et al. (2007) and is the same we used in Aresu et al. (2011). It is a bremsstrahlung spectrum with $I_v(E) \propto 1/E \times \exp(-E/kT)$ (see Fig. 1). We varied those quantities that could potentially affect the penetration of X-ray and FUV radiation and change the energy deposition through the disk. For example, we use two different values for the dust size distribution power law. The first value, generally used in the literature for the ISM, is a power law index of $a_{pow} = 3.5$. It is well known, though, that dust coagulation is the seed process that in timescales of a few million years leads to the formation of planetesimals. To model a disk that could have undergone dust coagulation, we also adopt a shallower power law index of $a_{pow} = 2.5$. For the same reason, we vary the minimum dust size: 0.1, 0.3, and 1 µm. These two parameters ($a_{min}$ and $a_{pow}$) determine the FUV opacity and, as a result, impact the amount of energy reprocessed by the photo-electric heating.
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Input spectrum for model with $L_X = 10^{30}$ erg s$^{-1}$ and $L_{FUV} = 10^{31}$ erg s$^{-1}$.</figcaption></figure>

Ultimately, the parameters affect the gas thermal balance in the disk. The changing opacity and optical depth also impact the dust temperature $T_{dust}$.

Aresu et al. (2011) performed an exploratory study on the influence of X-rays and FUV irradiation on the structure of protoplanetary disks, for X-ray luminosities ranging from $L_X = 10^{29}-10^{32}$ erg s$^{-1}$. The range of values for this particular parameter is based on the outcome of the Taurus survey (Güdel et al. 2007), which showed that X-ray-active T Tauri stars emit in the range $L_X = 10^{29}-10^{31}$ erg s$^{-1}$. The model with $L_X = 10^{32}$ erg s$^{-1}$, on the other hand, is an extreme case with interesting implications in terms of modeling. FUV excess emission from classic T Tauri stars between 6 and 13.6 eV is believed to arise from accretion spots on the surface of the star and stellar activity. The disk is truncated to several stellar radii from the star magnetic field, which channels matter toward the stellar photosphere. Shocks in the outer layers of the star cause the temperature to rise and emit FUV photons, which results in a FUV excess with respect to the stellar luminosity in the same band (Gullbring et al. 2000). Gorti & Hollenbach (2008) infer the FUV excess luminosity (91.2 nm < $h\nu$ < 205 nm) to be between $L_{FUV}/L_* \sim 10^{-2}-10^{-3}$. In our grid, we decided to scale the FUV flux so that we obtain the same range of luminosities as we use for the X-rays. The adopted values of $L_{FUV}/L_*$ are then $2.6 \times 10^{-2}$, $2.6 \times 10^{-3}$, $2.6 \times 10^{-4}$, and $2.6 \times 10^{-5}$. This way we can directly compare the input energy radiation in the FUV and X-ray band as well as its effects on the disk. We also explore two different values for the surface density power law distribution $\epsilon$, which is defined as $\Sigma(r) = \Sigma(r_0) \times (r/r_0)^{-\epsilon}$: Hayashi (1981) derived the value $\epsilon = 1.5$ in his model for the minimum mass solar nebular (MMSN) and their diagnostic value for our own solar system, while Hartmann (1998) suggested $\epsilon = 1$ for objects older than 1 Myr.

In summary, the current paper discusses the effects of variations in the following parameters: X-ray luminosity (five values), FUV luminosity (four), surface density profile (two), and dust size distribution, varying both the minimum grain size $a_{min}$ (three) as well as power law indices (two). This yields a total of 240 models; the summary of the model parameters is given in Table 1.

<table><caption>**Table 1.** Parameters used in the models.</caption><thead><tr><th>Quantity</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Stellar mass</td><td>M<sub>*</sub></td><td>1 M<sub>☉</sub></td></tr><tr><td>Effective temperature</td><td>T<sub>eff</sub></td><td>5770 K</td></tr><tr><td>Stellar luminosity</td><td>L<sub>*</sub></td><td>1 L<sub>☉</sub></td></tr><tr><td>Disk mass</td><td>M<sub>disk</sub></td><td>0.01 M<sub>☉</sub></td></tr><tr><td>X-ray luminosity (0.1–50 keV)</td><td>L<sub>X</sub></td><td>0, 10<sup>29</sup>, 10<sup>30</sup><br>10<sup>31</sup>, 10<sup>32</sup> erg s<sup>-1</sup></td></tr><tr><td>FUV luminosity</td><td>L<sub>FUV</sub></td><td>10<sup>29</sup>, 10<sup>30</sup>,<br>10<sup>31</sup>, 10<sup>32</sup> erg s<sup>-1</sup></td></tr><tr><td>Inner disk radius</td><td>r<sub>in</sub></td><td>0.5 AU</td></tr><tr><td>Outer disk radius</td><td>r<sub>out</sub></td><td>500 AU</td></tr><tr><td>Surface density power law index</td><td>ε</td><td>1.0, 1.5</td></tr><tr><td>Dust-to-gas mass ratio</td><td>ρ<sub>d</sub>/ρ</td><td>0.01</td></tr><tr><td>Min. dust particle size</td><td>a<sub>min</sub></td><td>0.1, 0.3, 1.0 μm</td></tr><tr><td>Max. dust particle size</td><td>a<sub>max</sub></td><td>10 μm</td></tr><tr><td>Dust size distribution power index</td><td>a<sub>pow</sub></td><td>2.5, 3.5</td></tr><tr><td>Dust material mass density</td><td>ρ<sub>gr</sub></td><td>2.5 g cm<sup>-3</sup></td></tr><tr><td>Strength of incident ISM FUV</td><td>χ<sup>ISM</sup></td><td>1</td></tr><tr><td>Cosmic ray ionization rate of H<sub>2</sub></td><td>ζ<sub>CR</sub></td><td>5 × 10<sup>-17</sup> s<sup>-1</sup></td></tr><tr><td>Abundance of PAHs relative to ISM</td><td>f<sub>PAH</sub></td><td>0.12</td></tr><tr><td>Viscosity parameter</td><td>α</td><td>0</td></tr></tbody></table>

### 3. Disk thermal and density structure

The coupling of gas and dust to the radiation field determines the heating and cooling rates, which in return determine the pressure balance and structure of a disk. FUV and X-rays couple differently to the gas: FUV is absorbed by dust grains and ejects an electron with a small surplus of energy in the form of kinetic energy, which then consequently heats the gas. X-rays, on the other hand, are absorbed by the K-, L-, or M-shell of atoms, where no distinction is made as to whether this atom is part of a molecule, dust grain, or PAH. This assumption might overestimate the X-ray absorption cross section by a factor of approximately ~2 at energies $E < 1$ keV, since we do not include self-shielding effects by large grains (see, e.g., Fig. 1 of Fireman 1974). On the other hand, not much is known about gas and dust phase elemental abundances in disks. The observed abundance of, e.g., neon (dominating the X-ray absorption cross section between ~0.9 to 1.2 keV, and, as such, the total absorption of our 1 keV thermal source), already varies over a larger range (see Glassgold et al. 2004, for a discussion on this topic). The electron that is ejected has a large energy capable of exciting, ionizing and heating the gas. The efficiency for X-ray heating is much larger, about 10 to 40 percent, depending on the ionization fraction of the gas, while that for FUV heating does not exceed one to three percent. The coupling of X-rays to gas is weaker than for FUV, due to the smaller cross section for the absorption, but the larger efficiency makes sure that it is an important heating source, even when the FUV excess is large. Their relative contribution to the energy budget will result in a different structure and emitted line fluxes. X-rays do not much affect the continuum fluxes, certainly less than $a_{min}$, $a_{pow}$, and $\epsilon$ (see Fig. 2). It shows that the fluxes change at most by a factor of three between $\lambda = 5-100$ μm due to X-rays. The models with $L_X = 10^{32}$ erg s$^{-1}$ are affected most, and these high X-ray fluxes are rarely observed for T Tauri stars.

Woitke et al. (2009) have already pointed out that the cooling and heating balance is very important in determining the physical structure of the disk. They present in their Fig. 9 the density structures for (1) a model where the gas temperature is decoupled from the dust and (2) where the gas and dust temperature are coupled. Model 1 shows a very complex structure, where the
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Average spectral energy distributions of models with $L_X = 0$ (red), $10^{29}$ (blue), $10^{30}$ (orange), $10^{31}$ (green), and $10^{32}$ erg s$^{-1}$ (yellow).</figcaption></figure>

gas is puffed up at the inner rim and a second bump in the density is seen at a radius $r \sim 10$ AU. Aresu et al. (2011) extended this model by including X-rays and showed for the particular model of Woitke et al. (2009) that the vertical extension inward of 10 AU becomes smoothed out for increasing X-ray fluxes and merges with the inner rim for the largest X-ray fluxes. A minimum in the vertical extension remains visible, though.

In Figs. 3 and 4, the density and temperature of the disk are shown for FUV luminosities ranging from $L_{FUV} = 10^{29}$, $10^{30}$, $10^{31}$, and $10^{32}$ erg/s (left to right) and for X-ray luminosities $L_X = 0$, $10^{29}$, $10^{30}$, $10^{31}$, and $10^{32}$ erg/s. There are 12 different possible variations with these X-ray and FUV luminosities (due to the variations in $a_{min}$, $a_{max}$, and $\epsilon$); as a baseline model, we always show the models with $a_{min} = 0.1$ µm, $\epsilon = 1.5$, and $a_{pow} = 3.5$.

## 3.1. Density distribution

Because density and temperature structure are directly related and we assume hydrostatic equilibrium, many of the effects that are seen in the density structure will show up in the temperature structure as well. The FUV-only models show that the inner rim becomes higher and broader for increasing luminosity. The disk right behind the inner wall is shielded, and a second, vertically more extended, density structure thus appears at radii between $r = 2$ and 6 to 30 AU (depending on the FUV luminosity).

The X-rays have a different effect on the density structure. While the peak of the second bump in the density structure appears at larger radii for higher FUV luminosities, the X-rays mostly affect the region within 5 AU. The second bump behind the inner rim expands to smaller and smaller radii for increasing X-ray luminosities (while keeping FUV the same), and the absolute densities become higher (from $n \sim 10^6$ to $10^8$ cm$^{-3}$). Eventually, the puffed-up inner rim and the second bump merge. At the highest FUV luminosity, the merging of the two density structures occurs at the smallest X-ray luminosity.

## 3.2. Temperature structure

When only FUV is included in the models, we find that the temperature structure of the disk shows an inversion in the vertical direction: the temperature reaches a maximum at a certain relative height $z/r$ and shows a fast drop above it. For example, the models with $L_{FUV} = 10^{29}$ erg/s show a region with temperatures $T > 1000$ K out to a radius $r \sim 3$ AU and between relative heights $z/r \sim 0.1$ to 0.4. Noticeable is also that the hot region extends to a slightly higher relative height at the inner rim (up to $z/r = 0.45$) and at a radius $r = 3$ AU (out to $z/r = 0.55$). This is the result of shielding; the disk surface is less heated right behind the inner rim. For larger FUV luminosities, shielding effects are even more prominent and the drop in temperature right behind the inner rim is increasingly pronounced.

Adding X-ray heating changes this picture. These models show a much more extended region with temperatures higher than $T > 1000$ K, even when the X-ray luminosity is much smaller than the FUV luminosity. At radii $r < 5-10$ AU, the inversion layer disappears and the temperature just smoothly increases toward larger $z/r$. At larger radii ($r > 10-20$ AU), the inversion layer only disappears when the X-ray luminosity is at least ten percent of the FUV luminosity. In models where the X-ray luminosity is much smaller, temperature distribution is also affected (or even dominated) by the FUV. An example is $L_{FUV} = 10^{32}$ erg/s in combination with $L_X = 10^{29}$ erg/s. Although the inversion layer disappears at small radii ($r < 5$ AU), it is still present at larger radii.

## 3.3. Heating and cooling processes

Temperature and density structure are in the end determined by the balance between the heating and cooling processes. Line cooling processes (except Lyman $\alpha$ and [OI] 6300 Å line cooling) are treated in an escape probability approximation. The code provides an option to calculate [OI] 6300 Å line cooling in an escape probability approximation as well, but this was not used in the present calculation. Before solving the equations for statistical equilibrium, a global continuum radiation transfer calculation is done to estimate the background mean intensities for the radiative excitation and de-excitation rates. Other heating and cooling processes are approximated by an analytical expression. Over 50 different heating and 50 different cooling processes are included in the code. Except for the treatment of X-ray related heating processes, which are described in Aresu et al. (2011), we refer to Woitke et al. (2009, 2011) for a more detailed description. The reason that so many processes are included is that they each play a significant role, depending on the ambient density and radiation field (FUV, X-rays, cosmic rays), which vary over many orders of magnitude. It is beyond the scope of this paper to describe all these processes in detail for the models. However, we do present the dominant heating and cooling processes in Figs. B.1 and B.2 for the models for which we are showing the temperature and density structure.

First, we consider the FUV-only models. At the lowest FUV luminosity $L_{FUV} = 10^{29}$ erg s$^{-1}$, there are three main heating process in the unattenuated part of the disk, namely, background heating by $C^+$ at the lowest densities, $n_H = 10^5$ cm$^{-3}$ (dark blue), PAH heating within the inner rim and the second extension of the disk (orange), and heating by carbon ionization in the outer disk (blue-green). Increasing the FUV field makes the picture more complicated. Heating by collisional de-excitation of $H_2$ becomes important in the second bump (light blue), while
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Total density: FUV luminosity increasing from $L_{\text{FUV}} = 10^{29}$ (left) to $10^{32}$ erg/s (right). X-ray luminosity increasing from $L_{\text{X}} = 0$ (top) to $10^{32}$ erg/s (bottom).</figcaption></figure>

background heating by FeII is dominant at the inner rim. Locally, other processes, such as photo-electric heating, are important. Going to the more shielded regions of the disk, there is roughly a three-layered structure of dominating heating processes: (1) infrared background heating CO ro-vibrational transitions (black), (2) slightly deeper in the disk heating by thermal accomodation on grains (white), and (3) in the mid-plane cosmic ray heating (red). This three-layered structure becomes more confined toward the mid-plane for higher FUV luminosities.

X-ray Coulomb heating dominates in the upper part of the disk, when the FUV to X-ray luminosity $L_{\text{FUV}}/L_{\text{X}} \le 1$. When it is higher than one, X-ray Coulomb heating is restricted to the higher regions of the disk. Adding X-rays increases the temperature significantly in the disk, making [FeII] also dominant in large parts of the second-density extension at the lowest X-ray luminosity, $L_{\text{X}} = 10^{29}$ erg s⁻¹. For large X-ray luminosities, the structure of heating sources becomes less complicated. The entire unshielded part of the disk is dominated by X-ray
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Gas temperature distribution: FUV luminosity increasing from $L_{\text{FUV}} = 10^{29}$ (left) to $10^{32}$ erg/s (right). X-ray luminosity increasing from $L_x = 0$ (top) to $10^{32}$ erg/s (bottom). White contours represent the $A_v = 1$ and 10. The black contours represent the gas temperatures.</figcaption></figure>

Coulomb heating, followed by the three-layered structure described above (CO ro-vib, thermal accommodation on grains, cosmic ray). For the highest two X-ray fluxes, an additional layer with X-ray H$_2$ dissociation heating as the main heating source is located on top of this three-layered structure.

In the FUV-only case, there are three main coolants in the unattenuated part of the disk: C$^+$ line cooling (yellow), [FeII] line cooling (blue-green), and [OI] line cooling. The size of the region, where [FeII] line cooling dominates expands for higher FUV luminosities, while reducing the size of region where [OI] line cooling dominates. When X-rays are added, the

region where C$^+$ line cooling dominates is pushed to the outer part of the disk. Temperatures are much higher in the upper part of the disk, and as a result Lyman $\alpha$ cooling (black) dominates in increasingly larger regions of disk, when X-ray luminosities become larger. At the highest X-ray luminosities, we find a three-layered structure of Lyman $\alpha$ cooling, [FeII] line cooling, and [OI] line cooling. The shielded region of the disk shows a layered structure. CO ro-vibrational (red) and H$_2$O rotational cooling (green-blue) are on top. Closer to mid-plane, several smaller regions have their own major coolant, such as HCN (purple), and HNC (orange).
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Scale height of the disk at $z/r = 0.5$. The black lines show the 12 different models ($a_{min}$, $a_{max}$, and $\epsilon$), at a particular $L_{FUV}$ and $L_X$. The vertical green lines indicate a maximum or minimum in the scale height. The red line is the average of the 12 models.</figcaption></figure>

### 3.4. Vertical scale height

Another way to show the combined effects of X-rays and FUV on the disk structure is by comparing the scale vertical height. Our definition is the same as the one used by Woitke et al. (2009), and is approximately (i.e., assuming $z \ll r$) given by

$$ \left(\frac{H}{r}\right)^2 \approx \frac{2rc_T}{GM_*} \qquad (1) $$

where $H$ is defined as $\rho(z) \approx \rho(0) \exp(z^2/H^2)$ and $c_T$ is the isothermal sound speed. As mostly the disk atmosphere is affected by the FUV and X-ray irradiation, we plot $H/r$ at the relative height $z/r = 0.5$ in Fig. 5. In each panel, we overplot the scale height for all the different parameters in black (twelve models in total), and the red line is the average. Four blue lines are overplotted, indicating flaring index $H/r \propto r^p$ with $p = 0.15-0.45$, where the lines are normalized to the vertical scale height at $r = 100$ AU.

A number of things stand out in these plots: (i) the vertical scale height of the inner rim is mostly affected by the X-rays, and almost no response is seen for different FUV luminosities. The scale height is $H/r = 0.15$ without X-rays at the inner rim, and it increases to ~0.4 at $L_X = 10^{32}$ erg/s. (ii) The second bump also shows up when plotting $H/r$ and slowly moves outward for increasing X-ray luminosities. (iii) The flaring index $p$ is 0.25 in the outer disk ($r > 100$ AU), similar to the 2/7 exponent found by Chiang & Goldreich (1997). The vertical scale height

in the inner disk is elevated, compared to the outer disk. (iv) The $H_0/r \approx 0.1$ and $0.18$ at $r_0 = 100$ AU and $z/r = 0.1$ and $0.5$, respectively. For a direct comparison to Chiang & Goldreich (1997), we measure $H_0/r$ at $r_0 = 10$ AU, because the approximation of the disk dust temperature being vertically isothermal does not hold in our models at the larger distances. We find a value $H_0/r \approx 0.04$, which is very similar to Chiang & Goldreich (1997), after correcting for the stellar mass and their definition of the scale height. (v) The merging of the inner rim with the second bump is clearly visible, as changes in the vertical scale height in the inner disk ($r = 1-10$ AU) become more gradual for larger X-ray luminosities.

## 4. The chemical balance in the disk

The combined effects of X-rays and FUV irradiation on the disk are discussed in this section. The density structure is in pressure equilibrium as discussed in the previous section and is altered due to irradiation effects. The chemical rates are density dependent, and as such the pressure balance also influences the abundances of the species.

### 4.1. Electron abundances

The main electron donor in a FUV-only chemistry is atomic carbon, because the incident FUV radiation field is not capable of ionizing atomic hydrogen, leaving cosmic rays as the only ion-
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Electron abundance: FUV luminosity is fixed at $10^{31}$ erg/s. X-ray luminosity increases from $L_x = 0$ (left) to $10^{32}$ erg/s (right). The black dashed line indicates $A_V = 1$. The blue lines mark the ionization fractions $x_e = 10^{-2}$ and $10^{-1}$.</figcaption></figure>

ization source of atomic hydrogen. As a result, the maximum electron abundance will not be higher than $x_e \sim 10^{-4}-10^{-3}$, as can be seen in Figs. 6 and B.3. The maximum electron abundances occur where the disk is unshielded to the radiation source. The inner rim has significant ionization fractions ($x_e \sim 10^{-4}-10^{-3}$) all the way down to the midplane of the disk, while further out the electron abundances are only this high at densities $n < 10^7$ cm$^{-3}$. The drop of the electron abundances below $x_e < 10^{-6}-10^{-5}$ for higher densities nicely coincides with the $A_V = 1$. Closer to the midplane, when the radiation becomes more shielded and densities are higher, the electron abundance exhibits a fast drop. This is the result of a combination of a decreasing ionization rate, $\zeta \propto n$, and an increasing recombination rate, $k_{\text{rec}} \propto n^2$. The abundance drop becomes more gradual at larger radii, which is a result of lower ambient densities and thus lower recombination rates. The region where the transition occurs is moving slightly toward smaller relative heights, when FUV fluxes are higher. The change in the transition location is not much due to the sharp increase of the density toward the midplane ($k_{\text{rec}} \propto n^2$): a relatively small increase in the density already compensates for larger FUV fluxes.

This picture changes quite significantly when the radiation fields also contain X-rays. The absorption of an X-ray photon by an arbitrary species creates a fast electron that can produce many ionizations, e.g., a 1 keV electron is able to produce ~27 hydrogen ionizations. As a result, the electron abundances are much higher when X-rays are included. In the unattenuated part of the disk, ionization fractions are of the order $x_e \sim 10^{-2}-10^{-1}$ (contours for $x_e = 10^{-2}$ and $10^{-1}$ are indicated). They can even be higher than $x_e > 10^{-1}$ at very small radii ($r < 0.5$ AU) and high relative heights ($z/r \gtrsim 0.5$). Here, low densities ($n_H < 10^7$ cm$^{-3}$) and high temperatures ($T \gg 5000$ K) reduce the recombination rates. Recombination rates decrease with temperature up to $T \sim 10\,000$ K because they are dominated by radiative recombination, which scales as $k_{\text{rec}} \propto n_e H_i (T/10^4)^{-X_{\text{rad}}}$, with $X_{\text{rad}} \sim 0.6-0.9$. We find slightly higher ionization fractions than, e.g., Glassgold et al. (2004). Their calculation at $r = 1$ AU (see their Fig. 4) is truncated at vertical column density $N_H \approx 2 \times 10^{18}$ cm$^{-2}$ and density $n_H \approx 10^7$ cm$^{-3}$. Ercolano et al. (2008) (their Figs. 1 and 2) show calculations for lower column densities and densities, and they find ionization fractions slightly higher than $x_e > 1$ at 0.07 AU at vertical column densities $N_H < 10^{16}$ cm$^{-2}$. Both calculations are consistent with those presented here. The electron abundance structure looks a little counterintuitive, especially for the lowest two X-ray fluxes $L_X = 10^{29}$ and $10^{30}$ erg/s. The region with very high ($x_e \approx 10^{-2}-10^{-1}$) becomes smaller for increasing FUV luminosities (Fig. B.3). This is because the inner rim is puffed up more and even merged with the second bump for the highest FUV fluxes. The densities are elevated to $n_H = 10^8$ cm$^{-3}$ in these regions, and the recombina-

tion rates are orders of magnitude higher as a result. The outer disk is shielded from radiation. The regions with very high ionization fractions ($x_e > 10^{-2}$) are of similar size for all FUV luminosities only for the highest X-ray luminosity (see Fig. B.3).

When the ionization fraction is increased by orders of magnitude, the gas shows a radically different chemistry, as it will be dominated by ion-molecule reactions. This will be discussed for some of the key species, and water in particular.

## 4.2. H and $H_2$ abundances

The abundance structure of atomic hydrogen is shown in the top panel of Fig. 7 (and also in Figs. B.4 and B.5). The highest abundances ($x_H \gtrsim 10^{-2}$) are obtained in the unshielded region of the disk, where FUV and X-rays have high fluxes. At radii $r < 0.6$ AU, the disk is atomic all the way down to the midplane of the disk, as it is directly exposed to the central source. At larger radii, the atomic fraction drops below $x_H \sim 10^{-2}$ and makes a transition to molecular hydrogen at a relative height of $z/r \sim 0.1$ between $r \sim 0.6-5$ AU. The transition occurs at increasingly larger relative heights at larger distances to the central source, and it is at approximately $z/r \sim 0.4-0.5$ (depending on the X-ray flux) at a distance $r \sim 100$ AU. The transition occurs at slightly smaller relative heights at larger X-ray fluxes.

The formation of $H_2$ on grains is extremely efficient (~1) when the dust temperatures are not high, $T_{\text{dust}} \sim 10-50$ K (Cazaux & Tielens 2004). At dust temperatures higher than $T_{\text{dust}} > 100$ K, the $H_2$ efficiency drops rapidly and $H_2$ formation on dust does not occur at $T_{\text{dust}} > 1000$ K. $H_2$ is also formed through the $H^-$ route, $H^- + H \rightarrow H_2 + e^-$, especially when X-rays are present, and the fractional abundance of electrons is relatively high (and thus also $H^-$). It is possible to maintain an efficient route to form molecular hydrogen in the gas phase at high temperatures, $T_{\text{gas}} > 300$ K and $T_{\text{dust}} > 100$ K (see contours in Fig. 7). Also FUV is much more efficient in destroying $H_2$, than X-rays. The FUV photo-dissociates $H_2$ ($H_2+$ FUV photon $\rightarrow H_2^+ \rightarrow H+H$). X-rays predominantly ionize molecular hydrogen indirectly by collisions with fast electrons produced after an X-ray absorption ($H_2 + e^* \rightarrow H_2^+ + e^- + e^*$), while only a small fraction of $H_2$ is dissociated in this process. After ionization, $H_2$ is able to reform $H_2$ through $H_2^+ + H \rightarrow H_2 + H^+$, but it will also be able to form $H_3^+$ through the reaction $H_2^+ + H_2 \rightarrow H_3^+ + H$. This particular species is key in forming molecules through ion-molecule reactions. Because X-rays are destroying $H_2$ less efficiently and can provide an environment to form $H_2$ in the gas-phase, a fast ion-neutral chemistry results, which makes it possible to form molecules (and maintain significant abundances) at much higher temperatures than when only FUV is present.
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Abundance structure of atomic (top) and molecular (bottom) hydrogen. FUV and X-ray fluxes are the same as in Fig. 6. The blue contours indicate the gas temperature $T_{\text{gas}} = 300$ and 1000 K. The black line indicates $T_{\text{dust}} = 100$ K.</figcaption></figure>

The bottom panel of Fig. 7 presents the abundances of H$_2$. It shows that the H$_2$ abundance in the warm atmosphere of the disk is largest for the highest X-ray fluxes (as high as $x_{H2} \sim 10^{-6}$ in the region directly exposed to the central source at $r = 0.6$ AU). This is the result not only of the processes described above, but also of the scale height of the inner rim, which increases with X-ray luminosities, an effect described in the previous section.

### 4.3. $H^+, H^-, H_2^+, \text{ and } H_3^+$ abundances

The main drivers of the chemistry in the atmospheres of disks exposed to irradiation are ionic species, starting with $H^+, H_2^+$, and $H_3^+$ (see also the discussion of the water chemistry). The $H^-$ abundance structure is shown for completeness, since it can be important in the production of H$_2$ (see Figs. 8 and B.6 to B.9).

The abundance structure of ionized hydrogen shows abundances as high as $x_{H^+} \sim 10^{-3}-10^{-2}$, even if there is only FUV irradiation. FUV is not able to photo-ionize atomic hydrogen, and hydrogen ionization is solely done by cosmic rays. This only happens at very low densities ($n < 10^4 \text{ cm}^{-3}$), though. Once densities of the order of $n \sim 10^5 \text{ cm}^{-3}$ are reached, the ionized hydrogen fractions are closer to $x_e \sim 10^{-5}-10^{-4}$ or smaller. The ionization balance is regulated by the charge exchange reaction $H + O^+ \leftrightarrow H^+ + O$. The backward reaction has an energy barrier of $T = 227$ K, so the ratio of $O^+$ with respect to $H^+$ decreases with temperature. The charge exchange rates are much faster than the photo-ionization reaction or the recombination rates of the species. The amount of $H^+$ is thus directly coupled to the production rate of $O^+$, which is also produced by cosmic ray ionization. Even when X-ray ionization determines the ionization balance of the gas, the charge exchange rates are much faster than the primary and secondary X-ray ionization rates of the species, and the ion abundances are directly coupled. X-rays produce much higher fractions of $H^+$ (as high as $x_{H^+} > 0.1$ at the lowest densities). The abundance structure strongly resembles that of the electrons, although the $H^+$ abundance drops faster closer to the mid-plane, and other species such as $Na^+$ take over as electron donor (see also Ádámkovics et al. 2011).

Once the gas gets more shielded to the radiation, it has higher abundances of molecular species (see, e.g., Fig. B.5). In those regions, $H^+$ can also be produced by other reaction paths, such as $CO^+ + H \rightarrow CO + H^+$ when only FUV is present and secondary X-ray ionization of H$_2$ when X-rays are present as well. Once cosmic ray ionization is the dominant source of ionization, $H^+$ is produced by cosmic ray ionization of H$_2$ or charge exchange with $He^+$.

The $H^-$ ion is produced by cosmic ray ionization ($H_2+CR \rightarrow H^+ + H^-$) or by radiative recombination ($H + e^- \rightarrow H^- + h\nu$). This last reaction is fairly slow at low temperatures due to temperature barriers, but it becomes important in the warm part of the atmosphere of the disk. When only FUV is irradiating the disk, the $H^-$ abundance is high ($x_{H^-} \sim 10^{-11}-10^{-10}$) at the inner rim and in the second bump extending out of the disk, and to a lesser extent in the transition zone from atomic to molecular gas. There are two reasons that there is not as much $H^-$ at high altitudes in the disk. The first one is that the electron densities are two to three orders of magnitude lower. The second one is that the ambient temperatures are between 50 and 100 K, which is not very favorable as the reaction rate contains an energy barrier. At the inner rim and in the second bump, both temperatures ($T > 1000$ K) and densities ($n > 10^7 \text{ cm}^{-3}$, and thus also electron densities) are much higher. When X-ray irradiation plays a role, the abundances structure changes a lot. The temperatures are in excess of $T > 1000$ K to large relative heights (see Fig. 4), and electron abundances are much higher (due to higher ionization rates and slower recombination rates). Especially when X-ray luminosities exceed $L_X > 10^{30} \text{ erg/s}$, the relative abundances of $H^-$ are one to three orders of magnitude higher (close to $x_{H^-} \sim 10^{-8}$). In the FUV case, the $H^-$ route to form H$_2$ is at least three orders of magnitude lower than the formation route on dust grains. It becomes a significant contributor (~30%) to the total H$_2$ formation rate at the highest X-ray luminosities when these extreme ($x_{H^-} \sim 10^{-8}$) abundances are reached, and at the same time, H$_2$ formation on dust is quenched by high dust temperatures.

$H_2^+$ is formed by cosmic ray ionization ($H_2 + CR \rightarrow H_2^+ + e^-$) in the regions where the gas is molecular and shielded from FUV
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Abundance structure of H+, H−, H+₂, and H+₃. FUV and X-ray fluxes are the same as in Fig. 6. Black contour indicates $n_H = 10^5$ cm⁻³.</figcaption></figure>

and X-ray radiation. When the H₂ abundance is high ($x_{H2} \sim 0.5$), however, the H+₂ is quickly transformed to H+₃ (H+₂ + H₂ → H+₃ + H). It is above this layer, that the H+₂ abundance is increasing. Here, the gas is more exposed to radiation (FUV or X-rays). If X-rays are present, H₂ is ionized through secondary ionizations. In the FUV-only case, where we do not have this last reaction, other routes contribute to the production of H+₂. One example is S⁺ + H₂ → H+₂ + S. In the unshielded regions of the disk, where H₂ abundances are low ($x_{H2} < 10^{-6}$), the only efficient way to form H+₂ is through radiative association H⁺ + H → H+₂ + *hv*. However, this is only efficient at relatively high temperatures (*T* > 5000 K) and significant abundance levels of H⁺. In Fig. 8, the cation H⁺ is shown to be more abundant by orders of magnitude when X-rays are present. It is obvious that the production of H+₂ is much more efficient as well, and hence abundances $x_{H+2} \sim 10^{-8}-10^{-7}$ occur in the high warm atmosphere of the disk. In the FUV-only case, there is only a significant abundance ($x_{H+2} > 10^{-12}$) in the inner rim and the second bump.

The dominant way to form H+₃ is through the ion-molecule reaction H+₂ + H₂ → H+₃ + H. As a result, the H/H₂ transition

layer and the region below is very suitable to produce large abundances: H₂ is present in reasonable amounts ($x_{H2} > 10^{-5}-10^{-3}$, see Fig. 7), and radiation is available to produce H+₂ (see previous paragraph, and Fig. 8). Once the H+₂ drops below an abundance $x_{H+2} < 10^{-12}$, the H+₃ abundance also drops to values below $x_{H+3} < 10^{-11}$. Consequently, there is a distinct layer where H+₃ is formed in the most optimal way. We also saw that the H+₂ is more efficiently formed when X-rays are abundantly present in the disk. The H+₃ abundances thus also significantly increase for larger X-ray luminosities and can be as high as $x_{H+3} \gtrsim 10^{-8}$. In the FUV-only case, the size of the region where H+₃ is present in large abundances is smaller, and the maximum abundance is at least an order of magnitude lower compared to the models that include X-rays. It is a key species in driving ion-molecule chemistry, and so the chemistry with and without X-rays will obviously yield significantly different abundance structures.

### 4.4. C²⁺, C⁺, C and CO abundances

The ionization potentials (IP) of the species C and C⁺ are IP = 11.26 and 24.38 eV, respectively. While FUV photons
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** C$^{2+}$, C$^{+}$, C, CO abundances. The black dashed line indicates $log(\chi/n_H) = -5$. Fluxes are the same as Fig. 6.</figcaption></figure>

can ionize neutral carbon, X-rays (or collisions with the resulting fast electrons) and/or cosmic rays are needed to reach a higher degree of ionization. As a result, the abundance patterns of the two species C$^+$ and C$^{2+}$ respond radically differently to the incident radiation field on the disk. The C$^{2+}$ abundance pattern is, of course, strongly correlated with the X-ray luminosity, and C$^{2+}$ is absent in models without X-rays. The region with significant C$^{2+}$ abundances ($x_{C^{2+}} > 10^{-9}-10^{-8}$) extends to larger radii and smaller relative heights for higher X-ray luminosities, simply because there is more ionizing radiation available (and not inhibited by dust, since the dominant opacity is caused by the gas). However, the FUV plays an important role in shaping the disk and therefore indirectly affects the abundance distribution as well. The bump in the total H number density (see Fig. 3) extends to larger relative heights for larger FUV luminosities, thereby increasing the density at larger relative heights. This causes a larger attenuating column between the central source and the outer part of the disk, thus reducing the amount of ionizing photons. Furthermore, it allows H$_2$ to form at smaller radii. H$_2$ is very efficient in reducing C$^{2+}$ by charge exchange reactions, such as C$^{2+}$ + H$_2$ → C$^+$ + H. As a result, the

FUV luminosity confines the C$^{2+}$ to smaller radii (see Fig. B.10). C$^+$ is present throughout the unattenuated part of the disk. It has an abundance $x_{C^+} \sim 10^{-4}$ down to the mid-plane in the inner rim ($r < 0.6$ AU). Beyond this radius, carbon is only significantly ionized above a relative height $z/r \sim 0.1$ (at $r \sim 0.6$ AU), increasing to $z/r \sim 0.4$ at $r \sim 200$ AU. The width $\Delta z/r$, over which the transition from C$^+$ to C and CO occurs, is larger toward the outer regions of the disk. Absolute densities of C$^+$ and e$^-$ in the inner disk are larger, and recombination rates scale with $n^2$. As a result, the transition from ionized to neutral carbon occurs abruptly, similar to the abundance drop of electrons (at approximately $z/r \sim 0.01$ at $r = 1$ AU). However, in the outer disk ($r = 200$ AU), the transition stretches out from $z/r \sim 0.04$ to $0.01$. Although there are variations of a factor a few in the absolute abundance of C$^+$, the overall appearance of the abundance structure is not significantly affected for different values of the X-ray and FUV luminosity (see Fig. B.11).

The neutral carbon abundance pattern is very much affected by variations in the X-ray and FUV luminosities. In the case where only FUV is irradiating the disk, one finds a clear transition from C$^+$ to C to CO (see Fig. 9), which is expected
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** O<sup>2+</sup>, O<sup>+</sup>, and O abundances. Fluxes are the same as in Fig. 6. The black contour in the lower panel indicates the critical density of the [OI] 63 µm line.</figcaption></figure>

in a FUV-dominated PDR (cf. Hollenbach & Tielens 1999). To guide the eye, we added a contour with value log(χ/n<sub>H</sub>) = -5 (with χ the Draine field). This contour follows the C<sup>+</sup>/C/CO transition very well, indicating PDR physics. Neutral carbon is confined in a layer between the C<sup>+</sup> and CO. This picture is more complicated when X-rays are added. Perhaps not entirely intuitive, the overall trend with increasing X-rays is that neutral carbon is more and more abundant in the unattenuated part of the atmosphere (as high as x<sub>C</sub> ~ 10<sup>-5</sup>-10<sup>-4</sup>). When X-rays are present, neutral carbon is higher by at least one to two orders of magnitude in abundance compared to models that exclude X-rays. The C and C<sup>+</sup> ionization rates are of course higher, but the recombination rates are also higher owing to the increased electron abundances (although the higher temperatures reduce the recombination rates, as already pointed out in Sect. 4.1). This is because electrons are now also produced by, e.g., ionization of hydrogen and helium. The fact that C<sup>+</sup> and C coexist in more or less equal amounts, when X-rays are dominating the ionization fraction of gas clouds, was already pointed out in papers by Maloney et al. (1996) (see their Figs. 3a and 4a) and Meijerink & Spaans (2005) (see their Figs. 3 and 4).

The bulk of the CO is situated in the shielded, optically thick part of the disk, and the total CO gas mass does not change when varying the FUV and X-ray radiation field. This molecule could thus serve as a tracer of the total mass of a protoplanetary disk, keeping in mind that the conversion factor depends on the amount of ice formation. A few details should be noted. There is a significant correlation between CO and OH in the unattenuated part of the disk (compare Figs. B.13 with B.17). A route to form CO is through the reaction C + OH → CO + H. While it is

possible to obtain an abundance level of x<sub>CO</sub> ~ 10<sup>-7</sup> through this channel, it will never go up to x<sub>CO</sub> ~ 10<sup>-4</sup> because the OH abundance is x<sub>OH</sub> ~ 10<sup>-7</sup>. As OH becomes less abundant for higher FUV luminosities (when X-rays are fixed), the CO also becomes less abundant and slightly more confined to the midplane. In addition, the average abundance of CO drops slowly with increasing X-ray luminosity at the inner rim by destruction through ion-molecule reactions, such as CO+He<sup>+</sup> → C<sup>+</sup>+O+He. This could have an effect on the ro-vibrational lines that are predominantly produced in these regions.

### 4.5. O<sup>2+</sup>, O<sup>+</sup>, and O abundances

The ionization potentials of O and O<sup>+</sup> are IP = 13.68 and 35.12 eV. Even the IP of neutral oxygen is above the threshold of IP = 13.6 eV, where neutral hydrogen blocks the radiation efficiently. The implications of this are seen in Fig. 10, where the abundances of neutral, singly and doubly ionized oxygen are shown for a FUV luminosity L<sub>FUV</sub> = 10<sup>31</sup> erg s<sup>-1</sup> and all considered X-ray luminosities except L<sub>X</sub> = 10<sup>32</sup> erg s<sup>-1</sup>. The second row of Fig. 10 (and also Fig. B.15 for variations with FUV luminosity) shows the abundance structure of singly ionized oxygen. Unlike the C<sup>+</sup>, the abundances of O<sup>+</sup> do not exceed x<sub>O+</sub> ~ 10<sup>-8</sup>-10<sup>-7</sup> in the FUV-only case. The abundance of ionized oxygen is mostly set by the charge exchange balance with atomic hydrogen, while the ionization of both hydrogen and oxygen is entirely due to cosmic ray ionization. When X-rays are added, the abundances increase three to four orders of magnitude compared to the FUV-only case, even at low X-ray luminosities (L<sub>X</sub> = 10<sup>29</sup> erg s<sup>-1</sup>). Important to note is that the O<sup>+</sup> is
confined to smaller radii for larger FUV luminosities, which is
also the case for C²⁺. This was discussed earlier in Sect. 4.4 and
explained by the larger recombination rates due to higher densi-
ties in the upper atmosphere of the disk and the lower ionization
rates due to shielding. The confinement to smaller radii is even
stronger for the case of O²⁺, because this species is not affected
by cosmic rays.

The bottom panel of Fig. 10 shows the neutral oxygen abun-
dance structure. The critical density is n crit = 5 × 10 5 cm⁻³
for the [OI] 63 µm line, which is indicated by a black dashed
contour. This species does not change with X-ray luminosity.
Neutral oxygen is the dominant oxygen carrier throughout the
disk, except for those regions, where large fractions of the gas
are frozen onto dust grains (i.e., the midplane). This species is
very insensitive to changes in FUV and the ambient chemistry
(see Fig. B.16). The bulk of atomic oxygen is in LTE for the
commonly observed [OI] 63 and 145 µm fine-structure lines and
thus sensitive to the average temperature of the region that it is
probing. As a result, the [OI] fine-structure lines have the poten-
tial to probe the total energy budget of the gas (the combined
FUV and X-ray luminosity irradiating the disk). There will be
an elaborate analysis of the oxygen fine-structure line emission
in Paper II.

4.6. *OH, OH<sup>+</sup>, H<sub>2</sub>O, H<sub>2</sub>O<sup>+</sup>, and H<sub>3</sub>O<sup>+</sup> abundances*

The water chemistry and its related species are significantly af-
fected by X-ray irradiation. The physical circumstances in the
disk change with radius *r* and height *z* due to radiation shield-
ing and large differences in chemistry, which means that the
X-ray and FUV energy deposition per particle, *H*<sub>X</sub>/n, and χ/n
(with χ = 1 the interstellar radiation field as defined by Draine
1978; Draine & Bertoldi 1996, and integrated between 91.2
and 205 nm), and their relative ratio change over orders of mag-
nitude in the disk. As a result, the main chemical pathways
forming the molecules OH and H<sub>2</sub>O change throughout the disk
and the grid.

X-rays heat the gas and therefore the reaction rates with ac-
tivation barriers increase. They also drive the ionization, which
enhances formation routes through ion-molecule reactions. The
formation of the neutral species OH and H₂O is either by the
neutral-neutral reactions, H₂ + O → OH + H, and H₂ + OH →
H₂O + H, or by recombination of ionized species, such as H₃O⁺ +
e⁻ → H₂O + H and H₂O⁺ + e⁻ → OH + H. Which of these path-
ways are dominating depends thus strongly on the temperature
and ionization fraction. The routes through neutral-neutral reac-
tions have activation barriers, making them efficient only when
gas temperatures are sufficiently high, i.e., T > 200–300 K. One
way to increase the neutral-neutral route without the need of
higher temperatures is by lowering the activation barrier, which
can be done by exciting H₂ to a higher vibrational state after ab-
sorption of a FUV photon. At high relative heights (z/r > 0.5),
this reaction will dominate when only FUV photons are present.
The ion-molecule chain is started with cosmic ray or X-ray ion-
ization of H₂ and thus requires significant ionization rates.

It turns out that the neutral-neutral reaction pathways dom-
inate in the inner part of the disk (r < 5 AU). The X-rays heat
regions deeper in the disk, resulting in the production of a thick
warm water layer. Although there is a larger fraction of X-rays
ionizing the gas than when only FUV photons are present,
the X-ray heating causes the neutral-neutral reactions to domi-
nate. It is remarkable to note that in the outer region (r ≥ 20 AU)
of the disk the ion-neutral reaction network dominates. This is
because temperatures are too low for the neutral-neutral reaction

pathways to occur. These combined effects become apparent in
the distribution of the water throughout the disk, as shown in the
second panel from the top in Fig. 11 (see also Fig. B.18 for a
more elaborate view on the effects of both FUV and X-rays). In
the models without any X-rays (left-hand side), the water shows
a strong abundance peak at the inner rim (all the way down to the
midplane of the disk) and a warm water layer at a relative height
z/r ≈ 0.1, with an abundance of xH₂O ~ 10⁻⁷, and a second wa-
ter reservoir higher up in the disk. This reservoir is located at a
relative height z/r ~ 0.015 at small radii (r ~ 1 AU), and at a
relative height of z/r ~ 0.5 at a radius of r ~ 100 AU, thus fol-
lowing the flaring of the disk. These different water reservoirs
were already noted by Woitke et al. (2009) for Herbig AeBe
stars. An elaborate discussion on the formation of warm water
reservoirs through (the dominant) neutral-neutral reactions and
formation on dust grains in the inner regions of an X-ray irradi-
ated disk is also outlined by Glassgold et al. (2009).

In the models where only FUV is included, the second layer
is confined to the regions where temperatures are between T ~
300–1000 K. The regions at higher temperatures are exposed
to FUV fluxes that essentially destroy the water faster than it
is formed. When X-rays are included, the second water reser-
voir extends to increasingly larger radii. Although there is wa-
ter present in the FUV-only models out to radii r ~ 200 AU,
the abundances are ultimately one to two orders of magnitude
higher when X-rays are added. The additional water is increased
since the ion-molecule chemistry is very effective in forming the
water. This is illustrated in the third to fifth panel of Fig. 11 (and
also Figs. B.19–B.21), where the changing abundance structures
of the ionic species OH⁺, H₂O⁺, and H₃O⁺ are shown for vari-
ous combinations of X-ray and FUV luminosities (Appendix A
only). The ionic species extend to larger radii for larger X-ray
fluxes, but larger FUV fluxes decrease the extent. The fact that
FUV confines the ionic species to smaller radii is an indirect
effect, because FUV tends to puff up the upper layers of the
disk and the density becomes higher at larger relative heights
(see Fig. 3). These larger densities increase the recombination
rates and reduce the abundances of these ionic species.

The first layer described above is located in those regions of
the disk, where the conventional routes do not work anymore.
Temperatures are too low and the ionization fraction is small.
This is the region where more exotic reactions take over, such as
NH₂ + NO → N₂ + H₂O. This layer has a smaller vertical extent
in relative height for larger FUV fluxes. On the other hand, the
second separate layer becomes thicker when X-rays are added.
When the Lₓ/L<sub>FUV</sub> ≥ 1 ratio, these two layers merge. The rea-
son for the merging is the enhancement of the transient species
OH⁺, H₂O⁺, and H₃O⁺ (see also Figs. B.19–B.21). As men-
tioned earlier, these species react sensitively to the presence of an
ionization source and become increasingly abundant with higher
X-ray fluxes.

The situation with OH is very similar. OH becomes more
abundant in the outer disk, when X-rays are added. OH is located
at higher altitudes in the disk (although there is also a large over-
lap with the regions where the water is located). For this reason,
the effect of the X-rays on the abundance in the outer disk is not
as pronounced as for water. The H₂ used in the neutral-neutral
formation route of OH is on average at higher temperatures,
so the relative contribution of the ion-molecule route is smaller.

4.7. *Ar⁺, Ar²⁺, Ne⁺, and Ne²⁺ abundances*

These species have ionization potentials that are larger than
that for atomic hydrogen, with IP = 21.56, 40.96, 15.76,
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** OH, H₂O, OH⁺, H₂O⁺, and H₃O⁺ abundances. Fluxes are the same as in Fig. 6. The white contours indicate gas temperatures of 250 and 2000 K. The red contour indicates the number density n_H = 10⁶ cm⁻³.</figcaption></figure>

and 27.63 eV for Ne, Ne⁺, Ar, and Ar⁺, respectively. The only way to produce these species is through direct (photon absorption and Auger effect) or indirect (fast electron collisions) ionization by X-rays. The presence of these ionized species is thus a direct result of the X-ray irradiation of the disk.

The models include a thermal source of temperature, Tₓ = 1 keV. As a result, the cross sections for neon are more favorable for direct ionizations than for argon. The cross section for the absorption is located at approximately ~1 keV for neon, while

for argon it is ~4 keV. On the other hand, the rates for secondary ionizations are a little higher for argon, σ(A)/σ(H) = 1.1, 0.48, 3.7, and 1.8, with A equal to Ne, Ne⁺, Ar, Ar⁺, respectively. Consequently, the total ionization rates are very comparable. However, there is one big difference in the chemical network of these species. The charge transfer rate of Ne⁺ with H₂ is very low (k < 10⁻¹⁴ cm³ s⁻¹), while the other ions have significant rates for charge exchange with molecular hydrogen. The Ne⁺ abundance structure thus extends to much smaller relative
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Abundance structure of selected heavy elements: Ne⁺ (top), Ne²⁺, Ar⁺, and Ar²⁺ (bottom). The FUV luminosity is fixed at $L_{\text{FUV}} = 10^{31} \text{ erg s}^{-1}$. X-ray luminosities vary between $L_X = 10^{29}$ and $10^{32} \text{ erg s}^{-1}$.</figcaption></figure>

heights, down to $z/r \sim 0.01$, than Ar⁺, which is abundant above $z/r \sim 0.01–0.06$ at radii $r \ge 10–200$ AU. The abundance structure of Ar⁺ is very similar to those for Ne²⁺ and Ar²⁺.

## 5. Radial column density profiles

Although the abundance distributions of species well illustrate the change in the chemistry due to different combinations of X-ray and FUV luminosities, it is necessary to look at the integrated properties of the disk that are observed by our telescopes. The focus in this is on the species C⁺, O, H₂O, and Ne⁺, as their line intensities and line profiles are extensively discussed in the Paper II. The figures show the average column density in yellow, while the black lines represent the 12 models at fixed FUV and X-ray luminosities. Changing parameters other than those for FUV and X-rays does not affect the results significantly.

C⁺ column density profile (Fig. B.26): the column density of C⁺ is strongly related to the inincident FUV flux at the inner rim, where it increases from $N_{\text{C}^+} \sim 10^{17}$ to $10^{19} \text{ cm}^{-2}$, when the FUV luminosity increases from $L_{\text{FUV}} = 10^{29}$ to $10^{32} \text{ erg s}^{-1}$.

Right after the peak at the inner rim, there is a sharp drop in the column density, which is caused by the inner rim casting a shadow and therefore reducing ionizing radiation. A second peak in the column density profile is seen between $r \sim 5–10$ AU. The peaks moves to larger radii for higher FUV luminosities. The X-rays affect the column density profile in a different way. They reduce the minimum column density right after the inner rim and smooth the radial distribution. The combination of the highest FUV and X-ray fluxes gives the flattest distribution.

O column density profile (Fig. 13): the column density profiles of neutral oxygen do not change as much as those for C⁺. The maximum column density, $N_O \sim 3–10 \times 10^{19} \text{ cm}^{-2}$, is located right behind the inner rim and only varies with a factor of at most ten, which is again caused by the FUV irradiation. The width of the oxygen column density peak broadens a bit for higher FUV and X-ray fluxes. This makes the minimum in the column density less prominent in the radial column density profile. Because the column density is not very affected by radiation, it is a clean probe of the properties of a disk temperature ([OI] fine-structure lines), since there are no strong dependencies on uncertainties in the chemical network.
<figure><img src="image_14.png" /><figcaption>Fig. 13. Radial column density distribution of O. FUV and X-ray fluxes are the same as Fig. 3.</figcaption></figure>

H₂O column density profile (Fig. 14): the highest water col-
umn densities, N<sub>H2O</sub> ~ 10<sup>20</sup>–10<sup>22</sup> cm<sup>-2</sup>, are found at the inner
rim, r ~ 0.5 AU. As discussed, the water abundances are en-
hanced by higher temperatures and higher ionization fractions
throughout the disk. The FUV counteracts this to some extent,
and the most favorable situation is a high X-ray to FUV luminos-
ity ratio. FUV irradiation puffs up the inner rim, and shields the
outer disk, thus allowing less ionizing radiation to penetrate into
the outer regions of the disk. As a result, it confines the water to
smaller regions of the disk. The bottom panel (L<sub>X</sub> = 10<sup>32</sup> erg s<sup>-1</sup>)
of Fig. 14 shows that the column density distribution is flat for
the lowest FUV flux, dropping steadily as a function of radius
for the highest FUV flux. The only region in the disk where the
column density becomes smaller is at the inner rim due to re-
actions with ions such as H<sub>2</sub>O + He<sup>+</sup> → H<sup>+</sup> + He + OH. The
water column density is almost two orders of magnitude smaller
for the models with the highest X-ray luminosities compared to
those with only FUV.

Ne⁺ column density profile (Fig. B.27): Ne⁺ is only pro-
duced by X-rays and therefore only the models that include
X-rays show significant column densities. The profiles show
a peak in the column density at the inner rim and a second

bump at a few AU. The second bump smoothes out for larger
X-ray fluxes, while the maximum column density increases from
N<sub>Ne+</sub> ~ 10<sup>15</sup> to 10<sup>17</sup> cm<sup>-2</sup>. The FUV tends to confine the Ne<sup>+</sup> to
smaller radii for the lower X-ray luminosities, which becomes
apparent in the line profile (see Paper II).

**6. Conclusions**

In this paper, we discussed the combined effects of FUV and X-rays on both the density structure and the thermal and chemical balance of disks around T Tauri stars for an expected range of parameters (dust size distribution, density profile, etc.), yielding a total of 240 models. Here we highlight the main results and a few implications:

6.1. *The disk thermal and chemical structure*

*Temperature structure:* the extent of the disk where the temper-
ature is higher than *T* > 1000 K is much larger when X-rays
are included. X-rays have a much higher heating efficiency than
FUV, 30–50% compared to <3%, respectively.
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Radial column density distribution of H~2~O. FUV and X-ray fluxes are the same as Fig. 3.</figcaption></figure>

**Density structure:** increasing FUV luminosities does not change the scale height of the inner rim; it only alters the width and height of the second bump in the disk that is created at intermediate radii ($r \sim 3–10$ AU), behind the region shielded by the puffed-up inner rim. Gas temperatures at the inner rim are much higher when X-rays are included and, as a result, the inner rim is puffed up to higher and higher altitudes for increasing X-ray luminosities. As a direct extension of this theoretical work, the existence of the second bump could potentially be tested by continuum imaging face-on protoplanetary disks in the near-infrared with, e.g., VLT or Keck.

**Scale height:** considering only FUV, we see that the scale height shows a maximum in the unattenuated parts ($z/r > 0.5$) of the disk. When X-rays are added we find that this maximum is smoothed over a larger region (out to $r \sim 10$ AU). The scale height in these regions is larger than one would expect from the flaring index in the outer regions of the disk. When moving to smaller relative height, e.g., $z/r = 0.1$, the break in the flaring index disappears. Another observational possibility would be to do continuum interferometry in the near-IR (VLTI) to directly measure the physical height of the inner rim.

## 6.2. Chemical balance

**Ionization fraction:** the ionization fraction reaches values as high as $x_{e^-} \sim 10^{-2}$ in exceptional cases in our FUV-only models, whereas X-rays easily maintain these ionization fractions throughout large portions of the disk. Even when the gas becomes partially shielded, it can still maintain a significant ionization fraction, leading to an ion-molecule chemistry that can form molecules at low temperatures, which is not possible with neutral-neutral reactions as they usually have temperature barriers.

**Formation of H~2~ through the H~−~ route:** formation of H~2~ dust is usually much more efficient in environments that have solar metallicities. Because of the high ionization fraction due to X-ray irradiation, the formation route H~−~ + H → H~2~ + e~−~ is able to provide a significant addition of order 50% percent to the H~2~ on dust formation route. Overall, the H~2~ to H abundance ratio is increased by at least two orders of magnitude when X-rays are present.

**Formation of water and OH:** the OH and H~2~O abundances are more concentrated toward the inner regions of the disk,
when only FUV is irradiating the disk. This is because the temperature is only there sufficient to drive the neutral-neutral formation route. The outer disk shows significant enhancements (up to two orders of magnitude) when X-rays are added. The higher ionization fractions make it possible to form the species through ion-molecule reactions and sustain abundance levels of $x_{H_2O} \sim 10^{-6}-10^{-7}$. Such abundance levels cannot be reached in outer disk models without X-rays. Only the warm inner disks allow even higher levels of water abundance through warm neutral-neutral chemistry.

Resulting abundance structures and radial column density profiles: whereas neutral oxygen and CO are very stable to both FUV and X-rays, this is not the case for, e.g., C+, Ne+, and H$_2$O. Ne+ is strongly enhanced by X-rays and confined to the inner regions by larger FUV luminosities. This latter aspect will certainly affect the line widths. A high $L_X/L_{FUV}$ ratio is favorable for water formation, especially in the outer disk.

## 6.3. Outlook

In Paper II, we will perform a radiation transfer analysis of the aforementioned species and correlate line fluxes and line widths directly to the FUV and X-ray luminosities. This will allow a discussion of the diagnostic value of these species and provide a theoretical framework for the interpretation of observational data. We will also discuss our results in the context of data obtained within several observational efforts (Spitzer, Herschel, and ground-based observing programs).

Acknowledgements. The research leading to these results received funding from the EU Seventh Framework Programme (FP7) in 2011 under grant agreement No. 284405.

## Appendix A: X-ray chemistry

Our chemical network is composed of $N_{sp} = 110$ species (see Table A.1), linked by ~1500 reactions.

For a given species $i$, the net formation rate reads (following Woitke et al. 2009):

$$ \frac{dn_i}{dt} = \sum_{jkl} R_{jk \to il}(T_g)n_j n_k + \sum_{jl} (R_{j \to il}^{ph} + R_{j \to il}^{cr} + R_{j \to il}^{Xpr} + R_{j \to il}^{Xsec}) n_j - n_i \left( \sum_{jkl} R_{il \to jk}(T_g) + \sum_{jk} (R_{i \to jk}^{ph} + R_{i \to jk}^{cr} + R_{i \to jk}^{Xpr} + R_{i \to jk}^{Xsec}) \right). \quad (A.1) $$

The terms involved are

- $R_{jk \to il}$, the temperature-dependent rate for a two-body reaction where species $i$ and $l$ are formed out of species $j$ and $k$;

- $R_{i \to jk}^{ph}$, a photo-reaction rate that depends on the local strength of the FUV radiation field;

- $R_{i \to jk}^{cr}$, a reaction that depends on the cosmic ray ionization rate;

- $R_{j \to il}^{Xpr}$ and $R_{j \to i}^{Xsec}$, the X-ray primary and secondary ionization reaction rates.

Assuming statistical equilibrium ($\frac{dn}{dt} = 0$), we obtain $N_{sp}$ nonlinear equations for the $N_{sp}$ unknown particle densities $n_k$,

$$ F_i(n_k) = 0. \quad (A.2) $$

For the species densities $n_k$, this system of nonlinear equations is solved through a Newton-Raphson iterative method, which expands $F_i$ into a Taylor series in the neighborhood of $n_k$,

$$ F_i(n_k + \delta n_k) = F_i(n_k) + \sum_j \frac{\partial F_i}{\partial n_k} \delta n_k + O(\delta n_k^2). \quad (A.3) $$

<table><caption>**Table A.1.** Species included in the chemical network.</caption><thead><tr><th colspan="2">110 species</th></tr></thead><tbody><tr><td>H, H<sup>+</sup>, H<sup>-</sup>, H<sub>2</sub>, H<sub>2</sub><sup>+</sup>, H<sub>2</sub><sup>-</sup>, H<sub>3</sub><sup>+</sup>, He, He<sup>+</sup>, O, O<sup>+</sup>, O<sup>2+</sup>, O<sub>2</sub>, O<sub>2</sub><sup>+</sup>, OH, OH<sup>+</sup>, H<sub>2</sub>O, H<sub>2</sub>O<sup>+</sup>, H<sub>3</sub>O<sup>+</sup>, CO, CO<sup>+</sup>, CO<sub>2</sub>, CO<sub>2</sub><sup>+</sup>, HCO, HCO<sup>+</sup>, H<sub>2</sub>CO, N, N<sup>+</sup>, N<sup>2+</sup>, NO, NO<sup>+</sup>, C, C<sup>+</sup>, C<sup>2+</sup>, CH, CH<sup>+</sup>, CH<sub>2</sub>, CH<sub>2</sub><sup>+</sup>, CH<sub>3</sub>, CH<sub>3</sub><sup>+</sup>, CH<sub>4</sub>, CH<sub>4</sub><sup>+</sup>, CH<sub>5</sub><sup>+</sup>, Si, Si<sup>+</sup>, Si<sup>2+</sup>, SiO, SiO<sup>+</sup>, SiH, SiH<sup>+</sup>, SiH<sub>2</sub><sup>+</sup>, SiOH<sup>+</sup>, S, S<sup>+</sup>, S<sup>2+</sup>, SO, SO<sup>+</sup>, SO<sub>2</sub>, SO<sub>2</sub><sup>+</sup>, OCS, CS, CS<sup>+</sup>, HS, HS<sup>+</sup>, HCS<sup>+</sup>, H<sub>2</sub>S<sup>+</sup>, H<sub>3</sub>S<sup>+</sup>, Mg, Mg<sup>+</sup>, Mg<sup>2+</sup>, Fe, Fe<sup>+</sup>, Fe<sup>2+</sup>, Ne, Ne<sup>+</sup>, Ne<sup>2+</sup>, Ar, Ar<sup>+</sup>, Ar<sup>2+</sup>, Na, Na<sup>+</sup>, Na<sup>2+</sup>, NH, NH<sup>+</sup>, NH<sub>2</sub>, NH<sub>2</sub><sup>+</sup>, N<sub>2</sub>H<sup>+</sup>, NH<sub>3</sub>, NH<sub>3</sub><sup>+</sup>, NH<sub>3</sub><sup>+</sup>, N<sub>2</sub>, NH<sub>2</sub><sup>+</sup>, CN, CN<sup>+</sup>, HCN, HCN<sup>+</sup>, HNC, C<sub>2</sub>H<sub>2</sub>, C<sub>2</sub>H<sub>3</sub><sup>+</sup>, HCNH<sup>+</sup>, CO#, H<sub>2</sub>O#, CO<sub>2</sub>#, CH<sub>4</sub>#, NH<sub>3</sub>#, PAH, PAH<sup>-</sup>, PAH<sup>+</sup>, PAH<sup>2+</sup>, PAH<sup>3+</sup></td></tr></tbody></table>

**References.** References for the rate coefficients: Millar et al. (1986); Lennon et al. (1988); Landini & Fossi (1991); Anicich (1993); Badnell (2006).

The code then needs to calculate the Jacobi matrix:

$$ J_{ik} = \frac{\partial F_i}{\partial n_k}. \quad (A.4) $$

Neglecting terms of the order of $\delta n_k^2$, the math SLATEC routines are used to find a set of $n_k$ that satisfies $F_i(n_k + \delta n_k) = 0$:

$$ J_{ik} \cdot \delta n_k = -F_i. \quad (A.5) $$

When the rate of a given reaction does not depend on the particle density $n_k$, the derivative of $F_i$ is straightforward, e.g., a photo-reaction that destroys species $i$:

$$ J_{ik} = \frac{\partial F_i}{\partial n_k} = -R_{i \to j}^{ph} \cdot \delta_{ik}. \quad (A.6) $$

X-ray reaction rates, however, depend on the electron density $n_e$, atomic hydrogen density $n_H$, and/or molecular hydrogen density $n_{H_2}$. These densities determine how much of the absorbed X-ray photons will go into the different channels for heating, ionization, and excitation. If the rate $R$ shows a dependency on the particle densities, for example, the electron density ($n_{el}$), atomic hydrogen density ($n_H$), and/or molecular hydrogen density ($n_{H_2}$), the Jacobian can be expressed as:

$$ J_{ik} = \frac{dF_i}{dn_k} = \frac{\partial F_i}{\partial n_k} + \frac{\partial F_i}{\partial n_{el}} \frac{\partial n_{el}}{\partial n_k} + \frac{\partial F_i}{\partial n_H} \frac{\partial n_H}{\partial n_k} + \frac{\partial F_i}{\partial n_{H_2}} \frac{\partial n_{H_2}}{\partial n_k}. \quad (A.7) $$

ProDiMo calculates this term analytically to ensure an accurate chemical solution. Here we describe how the Jacobian terms for the X-ray reactions rates are calculated. We define the following quantities, which will be used below:

$$ n_{\langle H \rangle} = \sum_j n_j Q(j, H) \quad (A.8) $$

$$ n_{el} = \sum_j n_j q_j, \quad (A.9) $$

where $n_{\langle H \rangle}$ is the total hydrogen nuclei density, $Q(j, H)$ is the stoichiometric coefficient for the hydrogen nuclei (which is 0 for all those species that do not contain hydrogen, while for species containing hydrogen we have: $Q(H, H) = 1$, $Q(H_2, H) = 2$, $Q(H_3^+, H) = 3$, etc.). The electron density is denoted as $n_{el}$ and $q_j$ is the charge of the particle $j$.
### A.1. Primary ionization

For the generic atomic species A, A⁺ and molecular species AB, these are reactions of the kind:

$$A + X_{ph} \rightarrow A^{2+} + 2e^{-} \qquad (A.10)$$

$$A^{+} + X_{ph} \rightarrow A^{2+} + e^{-} \qquad (A.11)$$

$$ \begin{aligned} & AB + X_{ph} \rightarrow A^{2+} + B + 2e^{-} \\ & \qquad \rightarrow A^{+} + B^{+} + 2e^{-} \\ & \qquad \rightarrow A + B^{2+} + 2e^{-}, \end{aligned} \qquad (A.12) $$

where $X_{ph}$ is an X-ray photon. When the $i$th species is molecular (A.12), the reaction rate $R_{i \rightarrow jl}^{Xpr}$ for the reaction that destroys the species $i$ has two indices $j$, and $l$, for the resulting species. Otherwise the reaction rate should simply read $R_{i \rightarrow jl}^{Xpr}$ as only one species is produced out of reactions (A.10) and (A.11). For simplicity, from now on, we use the molecular reaction rate as an example for X-ray primary ionization. The rate is calculated as

$$R_{i \leftarrow jl}^{Xpr} = \int_{E_i}^{\infty} \sigma_i(E) F(E, r) dE \quad [s^{-1}]. \qquad (A.13)$$

It depends on the X-ray radiation field $F(E, r)$ at the point $r$ where it is computed and on the cross section $\sigma_i$ of the species $i$. It does not depend on the local particle density. From Eqs. (A.3) and (A.6), we see that the contribution to the equilibrium equation for the species $i$ can be written as

$$F_i = F_{\text{oth}} - \sum_{jl} R_{i \leftarrow jl}^{\text{Xpr}} \cdot n_i + \sum_{jl} R_{j \rightarrow il}^{\text{Xpr}} \cdot n_j. \qquad (A.14)$$

The contribution to the Jacobi element of this reaction is then

$$J_{ik} = \frac{\partial F_{\text{oth}}}{\partial n_k} - \sum_{jl} R_{i \leftarrow jl}^{\text{Xpr}} \cdot \delta_{ik} + \sum_{jl} R_{j \rightarrow il}^{\text{Xpr}} \cdot \delta_{jk}. \qquad (A.15)$$

### A.2. Secondary ionization

We consider X-ray secondary ionization for all the atomic species A and only for a single molecule H₂:

$$A + e^{-} \rightarrow A^{+} + 2e^{-} \qquad (A.16)$$

$$H_2 + e^{-} \rightarrow H_2^{+} + 2e^{-} \qquad (A.17)$$

$$ \rightarrow H + H^{+} + 2e^{-}. \qquad (A.18) $$

The reaction rate for the X-ray secondary ionization of species $i$ is

$$R_{i \rightarrow j}^{\text{sec}} = r_i \frac{H_X n_{\langle \text{H} \rangle}}{W n_{\text{H}}} \quad [\text{s}^{-1}], \qquad (A.19)$$

where $r_i$ is a parameter that takes into account the geometrical cross section of species $i$ compared to the hydrogen cross section, $n_{\langle \text{H} \rangle}$ is the total hydrogen nuclei density, $H_X$ is the X-ray energy deposition (Maloney et al. 1996), and $W$ [eV] is the mean energy consumed per ion pair (Dalgarno et al. 1999).

W for hydrogen and all the other atomic elements except He is

$$W_{\text{H}} = W_{0,\text{H}} \left(1 + c_1 \left(\frac{n_{\text{el}}}{n_{\langle \text{H} \rangle}}\right)^{\alpha}\right) \left(1 + c_2 \frac{n_{\text{H}_2}}{n_{\text{H}}}\right) \quad [\text{eV}], \qquad (A.20)$$

where $W_{0,\text{H}}$ is the collisional ionization rate for H in a pure neutral atomic gas (13.6 eV), $n_{\text{el}}$ is the electron density, $n_{\text{H}_2}$

is the molecular hydrogen density, and $\alpha$, $c_1$ and $c_2$ are fitting parameters.

$W_{\text{H}}$ is then the energy needed to collisionally ionize hydrogen in a gas mixture with $n_{\text{el}}$, $n_{\text{H}}$ and $n_{\text{H}_2}$. The values for other atomic elements are scaled up, considering the geometrical factor $r_i = \sigma_i^{\text{coll}}/\sigma_{\text{H}}^{\text{coll}}$ (Ádámkovics et al. 2011). W has been calculated for He and H₂ as well; these cases will be treated separately further on.

### A.3. Hydrogen

The X-ray energy deposition is a quantity defined per unit of hydrogen nuclei as follows:

$$H_X = \int \sigma_{\text{tot}}(E) F(E, r) dE \quad [\text{erg} \langle \text{H} \rangle^{-1}], \qquad (A.21)$$

where $\sigma_{\text{tot}}$ is given by

$$\sigma_{\text{tot}}(E) = \sum_{i=1}^{N_{\text{sp}}} \sigma_i(E) n_i / n_{\langle \text{H} \rangle} \quad [\text{cm}^2], \qquad (A.22)$$

and $F(E, r)$ is the radiation field in unit of erg s⁻¹ cm⁻² eV⁻¹. For a given atomic species $i$, at the equilibrium, the contribution of the secondary ionization to the volumetric rate is

$$ \begin{aligned} F_i &= F_{\text{oth}} + F_i^{\text{sec}} \\ &= F_{\text{oth}} - R_{i \leftarrow j}^{\text{Xsec}} \cdot n_i + R_{j \rightarrow i}^{\text{Xsec}} \cdot n_j \\ &= F_{\text{oth}} - r_i \frac{H_X n_{\langle \text{H} \rangle}}{W_H n_{\text{H}}} n_i + r_j \frac{H_X n_{\langle \text{H} \rangle}}{W_H n_{\text{H}}} n_j. \end{aligned} \qquad (A.23) $$

Following Eq. (A.7), the derivative of $F_i^{\text{sec}}$ with respect to $n_k$ is

$$ \frac{dF_i^{\text{Xsec}}}{dn_k} = -r_i \times \left[ \frac{n_i}{W_H n_{\text{H}}} \frac{\partial (H_X n_{\langle \text{H} \rangle})}{\partial n_k} - \frac{H_X n_{\langle \text{H} \rangle} n_i}{W_{\text{H}}^2 n_{\text{H}}} \frac{\partial W_{\text{H}}}{\partial n_k} - \frac{H_X n_{\langle \text{H} \rangle} n_i}{W_H n_{\text{H}}^2} \delta_{\text{H},k} + \frac{H_X n_{\langle \text{H} \rangle}}{W_H n_{\text{H}}} \delta_{ik} \right]. \qquad (A.24) $$

Here we develop separately the partial derivative for all involved terms:

$$ \begin{aligned} \frac{\partial (H_X n_{\langle \text{H} \rangle})}{\partial n_k} &= \frac{\partial}{\partial n_k} \int \sum_j n_j \sigma_j(E) F(E, r) dE \\ &= \int \sigma_k(E) F(E, r) dE \\ &= D_K \end{aligned} \qquad (A.25) $$

$$ \begin{aligned} \frac{\partial W_{\text{H}}}{\partial n_k} &= \frac{\partial}{\partial n_k} \left[ W_{0,\text{H}} \left(1 + c_1 \left(\frac{n_{\text{el}}}{n_{\langle \text{H} \rangle}}\right)^{\alpha}\right) \left(1 + c_2 \frac{n_{\text{H}_2}}{n_{\text{H}}}\right) \right] \\ &= W_{0,\text{H}} c_1 \frac{\alpha q_k n_{\text{el}}^{\alpha-1}}{n_{\langle \text{H} \rangle}^{\alpha}} \left(1 + c_2 \frac{n_{\text{H}_2}}{n_{\text{H}}}\right) \\ &\quad - W_{0,\text{H}} c_1 \frac{\alpha n_{\text{el}}^{\alpha} Q_{k,\text{H}}}{n_{\langle \text{H} \rangle}^{\alpha+1}} \left(1 + c_2 \frac{n_{\text{H}_2}}{n_{\text{H}}}\right) \\ &\quad + W_{0,\text{H}} \left(1 + c_1 \frac{n_{\text{el}}^{\alpha}}{n_{\langle \text{H} \rangle}^{\alpha}}\right) \frac{c_2}{n_{\text{H}}} \delta_{\text{H}_2, k} \\ &\quad - W_{0,\text{H}} \left(1 + c_1 \frac{n_{\text{el}}^{\alpha}}{n_{\langle \text{H} \rangle}^{\alpha}}\right) c_2 \frac{n_{\text{H}_2}}{n_{\text{H}}^2} \delta_{\text{H}, k}. \end{aligned} $$
Substituting these results in Eq. (A.24) we obtain seven terms:

$$ \frac{dF_i^{\text{Xsec}}}{dn_k} = -r_i \frac{n_i}{W_H n_H} D_K \\ - r_i \frac{H_X n_{(H)}}{W_H n_H} \delta_{ik} \\ - r_i \frac{H_X}{W_H^2 n_H} n_i W_{0,H} c_1 \alpha \left( \frac{n_{el}}{n_{(H)}} \right)^\alpha \left( 1 + c_2 \frac{n_{H2}}{n_H} \right) Q(k, H) \\ - r_i \frac{H_X n_{(H)} n_i}{W_H^2 n_H^2} W_{0,H} \left( 1 + c_1 \left( \frac{n_{el}}{n_{(H)}} \right)^\alpha \right) c_2 \frac{n_{H2}}{n_H} \delta_{H,k} \\ + r_i \frac{H_X}{W_H^2 n_H} n_i W_{0,H} c_1 \alpha \left( \frac{n_{el}}{n_{(H)}} \right)^{\alpha-1} \left( 1 + c_2 \frac{n_{H2}}{n_H} \right) q_k \\ + r_i \frac{H_X n_{(H)} n_i}{W_H^2 n_H^2} W_{0,H} \left( 1 + c_1 \left( \frac{n_{el}}{n_{(H)}} \right)^\alpha \right) c_2 \delta_{H,k} \\ + r_i \frac{H_X n_{(H)} n_i}{W_H n_H^2} \delta_{H,k}. \qquad (A.26) $$

The seven terms can be interpreted as follows:

1. The first term represents the variation of the energy deposition $H_X$ when $n_k$ is increased. According to this term only, increasing $n_k$ brings more fast electrons in the gas phase via primary ionization of the species k. Hence more electrons are available for the secondary ionization of the species i.

2. The second term directly impacts $R_{i \leftarrow j}^{\text{Xsec}}$ if $i = k$ more particles are available for secondary ionization.

3. The third calculates the variations in W if the species considered contains hydrogen.

4. The fourth term comes from the dependency of $W_H$ on $n_H$ (the higher $n_H$ the closer is $W_H$ to $W_{0,H}$).

5. Analogous to the previous one, the fifth term shows that if the molecular hydrogen density increases (atoms are less likely to be ionized by fast electrons because) $W_H$ will be higher.

6. If the density of a positively charged particle increases, the electron density must increase at the same time as well. This favors Coulomb losses in energy of the incoming fast electrons over secondary ionization of species i. If $k = H^-$, this sixth term becomes negative because fewer electrons will be in the gas phase, favouring secondary ionization over Coulomb heating.

7. The seventh term comes directly from the secondary ionization rate. The higher the atomic hydrogen density, the lower is $R_{i \leftarrow j}^{\text{Xsec}}$.

### A.4. Molecular hydrogen and helium

The mean energy per ion pair for $H_2$ is

$$ W_{H2} = W_{0,H2} \left(1 + c_1 \left(\epsilon_{el}^*\right)^\alpha\right) \left(1 + c_2 \frac{n_H}{n_{H2}}\right), \qquad (A.27) $$

while for helium, it is

$$ W_{\text{He}} = W_{0,\text{He}} (1 + c_1 (\epsilon_{\text{el}})^\alpha) \qquad (\text{A.28}) $$

where $\epsilon_{el}^*$ is

$$ \epsilon_{el}^* = \frac{1.83 \epsilon_{el}}{1 + 0.83 \epsilon_{el}}, \qquad (A.29) $$

with $\epsilon_{el}$ is the electron fraction ($n_{el}/n_{(H)}$). Inserting these terms into (A.23) and generating the derivative with respect to $n_k$ gives the Jacobi elements for molecular hydrogen and helium respectively. Since they look very similar to Eq. (A.26), we will not write them out.

## References

Ádámkovics, M., Glassgold, A. E., & Meijerink, R. 2011, ApJ, 736, 143

Anicich, V. G. 1993, J. Phys. Chem. Ref. Data, 22, 1469

Aresu, G., Kamp, I., Meijerink, R., et al. 2011, A&A, 526, A163

Aresu, G., Meijerink, R., Kamp, I., et al. 2012, A&A, 547, A69

Badnell, N. R. 2006, ApJS, 167, 334

Brittain, S. D., Rettig, T. W., Simon, T., Gibb, E. L., & Liskowsky, J. 2010, ApJ, 708, 109

Cazaux, S., & Tielens, A. G. G. M. 2004, ApJ, 604, 222

Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368

D'Alessio, P., Canto, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411

Dalgarno, A., Yan, M., & Liu, W. 1999, ApJS, 125, 237

Draine, B. T. 1978, ApJS, 36, 595

Draine, B. T., & Bertoldi, F. 1996, ApJ, 468, 269

Dullemond, C. P., Dominik, C., & Natta, A. 2001, ApJ, 560, 957

Dullemond, C. P., Hollenbach, D., Kamp, I., & D'Alessio, P. 2007, Protostars and Planets V, 555

Ercolano, B., Drake, J. J., Raymond, J. C., & Clarke, C. C. 2008, ApJ, 688, 398

Fedele, D., Pascucci, I., Brittain, S., et al. 2011, ApJ, 732, 106

Fireman, E. L. 1974, ApJ, 187, 57

Glassgold, A. E., Najita, J., & Igea, J. 2004, ApJ, 615, 972

Glassgold, A. E., Najita, J. R., & Igea, J. 2007, ApJ, 656, 515

Glassgold, A. E., Meijerink, R., & Najita, J. R. 2009, ApJ, 701, 142

Gorti, U., & Hollenbach, D. 2008, ApJ, 683, 287

Güdel, M., Briggs, K. R., Arzner, K., et al. 2007, A&A, 468, 353

Güdel, M., Lahuis, F., Briggs, K. R., et al. 2010, A&A, 519, A113

Gullbring, E., Calvet, N., Muzerolle, J., & Hartmann, L. 2000, ApJ, 544, 927

Hartmann, L. 1998, Accretion Processes in Star Formation (Cambridge, UK, New York: Cambridge University Press)

Hashimoto, J., Tamura, M., Muto, T., et al. 2011, ApJ, 729, L17

Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35

Hogerheijde, M. R., Bergin, E. A., Brinch, C., et al. 2011, Science, 334, 338

Hollenbach, D. J., & Tielens, A. G. G. M. 1999, Rev. Mod. Phys., 71, 173

Jonkheid, B., Faas, F. G. A., van Zadelhoff, G.-J., & van Dishoeck, E. F. 2004, A&A, 428, 511

Kamp, I., & Dullemond, C. P. 2004, ApJ, 615, 991

Kamp, I., Woitke, P., Pinte, C., et al. 2011, A&A, 532, A85

Lahuis, F., van Dishoeck, E. F., Blake, G. A., et al. 2007, ApJ, 665, 492

Landini, M., & Fossi, B. C. M. 1991, A&AS, 91, 183

Leen, T. M., & Graff, M. M. 1988, ApJ, 325, 411

Lennon, M. A., Bell, K. L., Gilbody, H. B., et al. 1988, J. Phys. Chem. Ref. Data, 17, 1285

Maloney, P. R., Hollenbach, D. J., & Tielens, A. G. G. M. 1996, ApJ, 466, 561

Mathews, G. S., Dent, W. R. F., Williams, J. P., et al. 2010, A&A, 518, L127

Meeus, G., Montesinos, B., Mendigutia, I., et al. 2012, A&A, 544, A78

Meijerink, R., & Spaans, M. 2005, A&A, 436, 397

Meijerink, R., Glassgold, A. E., & Najita, J. R. 2008, ApJ, 676, 518

Merín, B., Brown, J. M., Oliveira, I., et al. 2010, ApJ, 718, 1200

Millar, T. J., Adams, N. G., Smith, D., Lindinger, W., & Villinger, H. 1986, MNRAS, 221, 673

Najita, J. R., Dopppmann, G. W., Bitner, M. A., et al. 2009, ApJ, 697, 957

Nomura, H., & Millar, T. J. 2005, A&A, 438, 923

Pinte, C., Woitke, P., Ménard, F., et al. 2010, A&A, 518, L126

Pontoppidan, K. M., Salyk, C., Blake, G. A., et al. 2010, ApJ, 720, 887

Pontoppidan, K. M., Blake, G. A., & Smette, A. 2011, ApJ, 733, 84

Riviere-Marichalar, P., Ménard, F., Thi, W. F., et al. 2012, A&A, 538, L3

Robitaille, T. P., Whitney, B. A., Indebetouw, R., Wood, K., & Denzmore, P. 2006, ApJS, 167, 256

Sturm, B., Bouwman, J., Henning, T., et al. 2010, A&A, 518, L129

Thalmann, C., Grady, C. A., Goto, M., et al. 2010, ApJ, 718, L87

Thi, W.-F., Woitke, P., & Kamp, I. 2011, MNRAS, 412, 711

van Zadelhoff, G.-J., Aikawa, Y., Hogerheijde, M. R., & van Dishoeck, E. F. 2003, A&A, 397, 789

Woitke, P., Kamp, I., & Thi, W.-F. 2009, A&A, 501, 383

Woitke, P., Riaz, B., Duchêne, G., et al. 2011, A&A, 534, A44