© ESO, 2017

1. Introduction

Our Sun acts as a particle accelerator and produces energetic particles with energies ≥10 MeV (e.g. Mewaldt et al. 2007). Such particles are also called solar cosmic rays as their energies are comparable to Galactic cosmic rays. They are accelerated in highly violent events like flares and/or close to the solar surface due to shocks produced by coronal mass ejections (Reames 2015). Therefore the energetic particle flux is strongly correlated with the activity of the Sun (e.g. Mewaldt et al. 2005; Reedy 2012).

From X-ray observations of T Tauri stars we know that their X-ray luminosities can be up to 10⁴ times higher than the X-ray luminosity of the contemporary Sun (e.g. Feigelson & Montmerle 1999; Güdel et al. 2007). Such high X-ray luminosities are rather a result of enhanced flare activity of young stars than coronal effects (e.g Feigelson et al. 2002). Enhanced activity of T Tauri stars implies an increase of their stellar energetic particle (SP) flux. From simple scaling with the X-ray luminosity and considering that young stars produce more powerful flares, Feigelson et al. (2002) derived a typical SP flux for T Tauri stars ≈10⁵ times higher than for the contemporary Sun. Under these assumptions T Tauri stars might show on average a continuous proton flux of f_p(E_p ≥ 10 MeV) ≈ 10⁷ protons cm^-2 s^-1 at a distance of 1 au from the star.

Such a scenario is also likely for the young Sun. Measurements of decay products of short-lived radionuclides (SLR) like ¹⁰Be or ²⁶Al in meteorites indicate an over-abundance of SLRs in the early phases of our solar system (e.g. Meyer & Clayton 2000). One likely explanation for these abundance anomalies are spallation reactions of SPs with the dust in the protosolar nebula (e.g. Lee et al. 1998b; McKeegan et al. 2000; Gounelle et al. 2001, 2006). However, this would require a strongly enhanced SP flux of the young Sun by a factor ≳3 × 10⁵ compared to the contemporary Sun (McKeegan et al. 2000), consistent with the estimated SP flux of T Tauri stars derived by Feigelson et al. (2002).

More recently Ceccarelli et al. (2014) reported a first indirect measurement of SPs in the protostar OMC-2 FIR 4. They observed a low HCO⁺/N₂H⁺ abundance ratio between three and four that requires high H₂ ionization rates >10^-14 s^-1 throughout the protostellar envelope. They explain this high ionization rate by the presence of SPs. From this ionization rate they derived a particle flux of f_p(E_p ≥ 10 MeV) ≥ 3 − 9 × 10¹¹ protons cm^-2 s^-1 at 1 au distance from the star. Such a high flux would be more than sufficient to explain the over-abundance of SLRs in the solar nebula.

A further indication of SPs in young stars is the anti-correlation of X-ray fluxes with the crystalline mass fraction of the circumstellar dust found by Glauser et al. (2009). They argue that SPs are responsible for the amorphization of dust particles and that the correlation can be explained if the SP flux scales with the stellar X-ray luminosity. Trappitsch & Ciesla (2015) tested such a scenario by using detailed models of SP transport for a protoplanetary disk (i.e. SP flux as a function of height of the disk). They also considered the vertical “mixing” of dust particles. According to their models SP irradiation of the disk cannot explain the total SLR abundances in the solar nebula but might play a role for dust amorphization.

Like Galactic cosmic rays SPs not only interact with the solid component but also with the gas component of the disk. However, little is known about the impact of SPs on the chemical structure of disks. Turner & Drake (2009) investigated the relevance of SPs on the size of dead-zones in disks. Assuming similar enhancement factors as mentioned above they find that SPs can decrease the size of the dead zone depending on the disk model and other ionization sources.

We present a first approach to study the impact of SP ionization on the chemical structure of the disk. We assume a typical SP flux as proposed for T Tauri stars to study the impact of SP ionization on the common disk ionization tracers HCO⁺ and N₂H⁺. We use the radiation thermo-chemical disk model PRODIMO (PROtoplanetary DIsk MOdel, Woitke et al. 2009, 2016; Kamp et al. 2010; Thi et al. 2011) to model the thermal and chemical structure of the disk. We argue that spatially resolved radial intensity profiles of molecular ion emission of the disk allow to constrain the SP flux of T Tauri stars.

In Sect. 2 we describe our method to derive the ionization rate due to SPs and the disk model. Our results are presented in Sect. 3 where we show the impact of SPs on the common disk ionization tracers HCO⁺ and N₂H⁺. In Sect. 4 we discuss possibilities to constrain the SP flux via observations of molecular ion emission and future prospects for modelling of SP ionization in protoplanetary disks. We present a summary and our main conclusions in Sect. 5.

2. Method

To investigate the impact of SPs on the disk chemical structure we first needed to determine the SP flux and the particle energy distribution. With these particle spectra we were able to calculate the ionization rate throughout the disk. We applied this to a disk structure representative for disks around T Tauri stars. With the radiation thermo-chemical disk code PRODIMO we calculated the chemical abundances. We did this for a series of models in which we also considered other important high energy ionization sources like Galactic cosmic rays (CR) and X-rays.

2.1. Stellar energetic particle spectra

As the actual particle spectra and fluxes of young stars are unknown, we derived the spectra from the knowledge available from our Sun. The origin of solar energetic particles are most likely flares and/or shock waves driven by coronal-mass ejections (CME; Reames 2013, 2015). Flares act like point sources on the solar surface whereas the shock waves can fill half of the heliosphere at around 2 Solar radii (Reames 2015). Particle fluxes are not continuous but rather produced in events lasting from several hours to days (Feigelson et al. 2002; Mewaldt et al. 2005).

Based on observed X-ray luminosities of solar analogs in the Orion nebula Feigelson et al. (2002) estimated that SP fluxes in young stars are likely ≈10⁵ times higher than in the contemporary Sun (see also Glassgold et al. 2005). As T Tauri stars are very active Feigelson et al. (2002) argue that it is likely that X-ray flares with luminosities below the detection limit occur several times a day (the same argument holds for CMEs). In that case the X-ray flares and consequently also SP events overlap, resulting in an enhanced continuous SP flux.

Based on these arguments we have assumed here a continuous and enhanced SP flux for young T Tauri stars. This approximation is consistent with the assumption of powerful and overlapping flare and CME events of T Tauri stars.

In Mewaldt et al. (2005) measurements of five different solar particle events are reported. We used their fitting formulae (see their Eq. (2) and Table 5) and derive SP spectra (protons in this case) averaged over the duration of the observed events. These measurements are for particles with energies up to several 100 MeV. We extrapolated their results up to energies typical for Galactic cosmic rays of ≈10 GeV. This is consistent with the maximum energy E_max ≈ 30 GeV derived by Padovani et al. (2015, 2016) for particles accelerated on protostellar surfaces. The resulting spectra are shown in Fig. 1. The flux levels for the different events can vary by up to two orders of magnitude and there is also some variation in the shape of the spectra.

Reedy (2012) reported proton fluxes of the contemporary Sun for five solar cycles. Typical values for the cycle averaged fluxes at a distance of 1 au are f_p(E_p> 10 MeV) = 59 − 213 protons cm^-2 s^-1. For the Event1 spectrum in Fig. 1 we get f_p(E_p> 10 MeV) = 151 protons cm^-2 s^-1 at 1 au, very similar to the reported values of Reedy (2012).

Here we used only the Event1 spectrum, which we call the “active Sun” spectrum. We simply scaled the active Sun spectrum by a factor of ≈10⁵, as proposed by Feigelson et al. (2002), to get a typical “active T Tauri” spectrum. The resulting SP flux of f_p(E_p> 10 MeV) = 1.51 × 10⁷ protons cm^-2 s^-1 is consistent with the value of ≈10⁷ protons cm^-2 s^-1 derived by Feigelson et al. (2002) for young solar analogs in the Orion Nebula cluster. The two SP spectra are shown in Fig. 2, where we also show two cases of Galactic cosmic-ray spectra for comparison (see Sect. 2.3.2).

For comparison we also present models applying the same approach for the treatment of SPs as Turner & Drake (2009). For their model they assumed that SPs behave very similar to Galactic CRs (e.g. particle energies). The details of the Turner model are discussed in Appendix D.

It is not clear if SPs actually reach the disk (see Feigelson et al. 2002, for a discussion). However, as we are interested in the possible impact of SPs on the chemistry, we simply assumed that all SPs reach the disk. We further discuss this assumption in Sect. 4.3.

2.2. Stellar particle transport and ionization rate

Energetic particles hitting the disk interact with its gas and dust contents. Although the interaction with the solids is relevant for the production of SLR, here we are only interested in the interaction with the gas. Dust only plays a minor role in the actual attenuation of particles as only ≈1% of the mass in protoplanetary disk is in solids (see also Trappitsch & Ciesla 2015). Very similar to Galactic cosmic rays, SPs mainly ionize the gaseous medium (i.e. molecular hydrogen). Energetic particles interact multiple times and ionize many atoms and molecules on their way until they eventually have lost their energy completely. This complex process requires detailed particle transport models.

To model the transport of energetic particles through the disk gas, we used the continuous slowing down approximation, which assumes that particles lose an infinitesimal fraction of their energy during propagation (Takayanagi 1973). We used the results obtained by Padovani et al. (2009, 2013a) who compute the propagation of CRs in a 1D slab, taking all the relevant energy loss processes into account. They give a useful fitting formula for the ionization rate as a function of the column density of molecular hydrogen.

In order to apply the results of this 1D transport model, we assumed that SPs travel along straight lines (i.e. no scattering due to their high energies) and that they originate from a point source (the star). We also disregarded the effect of magnetic fields that could increase or decrease the ionization rate depending on their configuration (Padovani & Galli 2011; Padovani et al. 2013b, see also Sect. 4.3.2).

From the detailed 1D particle transport model we derived a simple fitting formulae for the two SP spectra considered here. The SP ionization rate ζ_SP for molecular hydrogen as a function of the total hydrogen column density N_{⟨ H ⟩} = N_H + 2 × N_H2 is given by ${\mathit{\zeta }}_{\mathrm{SP}}\mathrm{\left(}{\mathit{N}}_{\mathrm{⟨}\mathrm{H}\mathrm{⟩}}\mathrm{\right)}\mathrm{=}{\left[\frac{\mathrm{1}}{{\mathit{\zeta }}_{\mathrm{L}}{\left(\frac{{\mathit{N}}_{\mathrm{⟨}\mathrm{H}\mathrm{⟩}}}{{\mathrm{10}}^{\mathrm{20}}\hspace{0.17em}\mathrm{c}{\mathrm{m}}^{-2}}\right)}^{\mathit{a}}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{\zeta }}_{\mathrm{H}}{\left(\frac{{\mathit{N}}_{\mathrm{⟨}\mathrm{H}\mathrm{⟩}}}{{\mathrm{10}}^{\mathrm{20}}\hspace{0.17em}\mathrm{c}{\mathrm{m}}^{-2}}\right)}^{\mathit{b}}}\right]}^{-1}\mathrm{\left[}{\mathrm{s}}^{-1}\mathrm{\right]}\mathit{,}$(1)and for N_{⟨ H ⟩}>N_E by ${\mathit{\zeta }}_{\mathrm{SP}\mathit{,}\mathrm{E}}\mathrm{\left(}{\mathit{N}}_{\mathrm{⟨}\mathrm{H}\mathrm{⟩}}\mathrm{\right)}\mathrm{=}{\mathit{\zeta }}_{\mathrm{SP}}\mathrm{\left(}{\mathit{N}}_{\mathrm{⟨}\mathrm{H}\mathrm{⟩}}\mathrm{\right)}\mathrm{\times }\exp \left[\mathrm{-}\left(\frac{{\mathit{N}}_{\mathrm{⟨}\mathrm{H}\mathrm{⟩}}}{{\mathit{N}}_{\mathrm{E}}}\mathrm{-}\mathrm{1.0}\right)\right]\mathrm{\left[}{\mathrm{s}}^{-1}\mathrm{\right]}\mathit{.}$(2)The two power laws in Eq. (1) are a consequence of the shape of the SP input spectra (see Fig. 2). The two parts of Eq. (1) account for the ionization rate at low (ζ_L) and high (ζ_H) column densities. Equation (2) accounts for the exponential drop of the SP ionization rate starting at a certain column density given by N_E (i.e. similar to CRs). For the two SP spectra considered here N_E = 2.5 × 10²⁵ cm^-2. The other fitting parameters ζ_L, ζ_H, a and b are given in Table 1.

Equations (1) and (2) provide the unattenuated SP ionization rate at 1 au distance from the star (i.e. for a SP flux at 1 au). To account for geometric dilution we scaled ζ_SP(N_{⟨ H ⟩}) by 1 /r² at every point in the disk (r is the distance to the star in au). For the chemistry we simply added ζ_SP to the ionization rate for Galactic cosmic rays ζ_CR (see Sect. 2.3.2).

2.3. Other ionization sources

To investigate the impact of SPs on the disk ionization structure also other ionization sources common to T Tauri stars must be considered. Besides SPs our model includes stellar UV and X-ray radiation, interstellar UV radiation and Galactic cosmic rays (CRs). However, most relevant for our study are the high energy ionization sources capable of ionizing molecular hydrogen: SPs, X-rays and CRs.

2.3.1. X-rays

To model the stellar X-ray spectrum we used an approximation for an isothermal bremsstrahlung spectrum (Glassgold et al. 1997; Aresu et al. 2011) $\mathit{F}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{\propto }\frac{\mathrm{1}}{\mathit{E}}\exp \mathrm{\left(}\mathrm{-}\mathit{E}\mathit{/}\mathit{k}{\mathit{T}}_{\mathrm{X}}\mathrm{\right)}\mathit{,}$(3)where E is the photon energy (here in the range of 0.1 to 20 keV), k the Boltzmann constant and T_X is the plasma temperature. This spectrum is scaled to a given total X-ray luminosity L_X (0.3 ≤ E ≤ 10 keV, e.g. Güdel et al. 2010).

As a result of the activity of the stars (e.g. flares) X-ray radiation of young stars is variable. We account for this in a simple way by including a spectrum with a X-ray luminosity and temperature representative for a typical T Tauri star, and a spectrum which represents a flaring spectrum (more activity) with higher luminosity and a harder (hotter) radiation. However, we actually ignored the variability and assume time averaged X-ray fluxes. According to Ilgner & Nelson (2006) this is a reasonable assumption as typically the recombination timescales in the disk are longer than the flaring period and large parts of the disk (r ≳ 2 au in their model) respond to an enhanced average X-ray luminosity.

We also included the stellar X-ray properties used in Turner & Drake (2009) for the Turner models (Appendix D). The parameters for the various X-ray input spectra are given in Table 2 and the spectra are shown in Fig. 3.

To derive the X-ray ionization rate ζ_X we used X-ray radiative transfer including scattering and a detailed treatment of X-ray chemistry (Aresu et al. 2011; Meijerink et al. 2012). For more details on the new X-ray radiative transfer module in PRODIMO see Appendix A.

2.3.2. Galactic cosmic rays

Protoplanetary disks are exposed to Galactic cosmic rays (CR). Differently to SPs and X-rays, CRs are not of stellar origin and hit the disk isotropically. Cleeves et al. (2013) proposed that for T Tauri disks the actual CR ionization rate might be much lower compared to the interstellar medium (ISM) due to modulation of the impinging CRs by the heliosphere (“T-Tauriosphere”).

From modelling molecular ion observations of the TW Hya disk, Cleeves et al. (2015) derived an upper limit for the total H₂ ionization rate of ζ ≲ 10^-19 s^-1. This upper limit applies for all ionization sources including SLRs (e.g. Umebayashi & Nakano 2009). SLR ionization is a potentially important ionization source in the midplane of disks. However, similar to Cleeves et al. (2015) we modelled the low ionization rate scenario by reducing the CR ionization rate and did not explicitly treat SLR ionization in the models presented here.

We considered two different CR input spectra, the canonical local ISM CR spectrum (Webber 1998) and a modulated CR spectrum which accounts for the exclusion of CRs by the “T-Tauriosphere” (Cleeves et al. 2013, 2015). For simplicity we call these two spectra “ISM CR” and “low CR”, respectively. To calculate the CR ionization rate ζ_CR in the disk we applied the fitting formulae provided by Padovani et al. (2013a) and Cleeves et al. (2013): ${\mathit{\zeta }}_{\mathrm{CR}}\mathrm{\left(}{\mathit{N}}_{\mathrm{⟨}\mathrm{H}\mathrm{⟩}}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\zeta }}_{\mathrm{l}}\hspace{0.17em}{\mathit{\zeta }}_{\mathrm{h}}}{{\mathit{\zeta }}_{\mathrm{h}}\mathrm{[}{\mathit{N}}_{\mathrm{⟨}\mathrm{H}\mathrm{⟩}}\mathit{/}{\mathrm{10}}^{\mathrm{20}}\hspace{0.17em}\mathrm{c}{\mathrm{m}}^{-2}{\mathrm{]}}^{\mathit{a}}\mathrm{+}{\mathit{\zeta }}_{\mathrm{l}}\mathrm{[}\exp \mathrm{(}\mathrm{\Sigma }\mathit{/}{\mathrm{\Sigma }}_{\mathrm{0}}\mathrm{)}\mathrm{-}\mathrm{1}\mathrm{]}}\mathrm{·}$(4)For simplicity we have assumed that CRs enter the disk perpendicular to the disk surface. Therefore we used the disk vertical hydrogen column density N_{⟨ H ⟩ ,ver} and surface density Σ_ver for Eq. (4) to calculate ζ_CR at every point in the disk. The fitting parameters ζ_l, ζ_h, Σ₀ and a for the two CR input spectra are given in Table 3 (see Padovani et al. 2009, 2013a, for details). The typical resulting H₂ ionization rates in the disk are ζ_CR ≈ 2 × 10^-17 s^-1 for the ISM CR spectrum and ζ_CR ≈ 2 × 10^-19 s^-1 for the low CR spectrum (see Sects. 3.1 and 3.2).

2.4. Disk model

To model the disk we used the radiation thermo-chemical disk code PRODIMO (Woitke et al. 2009, 2016; Kamp et al. 2010; Thi et al. 2011). PRODIMO solves the wavelength dependent continuum radiative transfer which provides the disk dust temperature and the local radiation field. The gas temperature (the balance between heating and cooling) was determined consistently with the chemical abundances. The chemical network included 235 different species and 3143 chemical reactions (see Appendix B for more details).

We used a disk model representing the main properties of a disk around a typical T Tauri star. The stellar properties and the disk structure of this model are identical to the so-called reference model presented in Woitke et al. (2016). Here we only provide a brief overview of the disk model and refer the reader to Woitke et al. (2016) for details.

In Fig. 4 we show the gas number density, dust density, the local far-UV (FUV) radiation field, gas temperature, dust temperature and the vertical hydrogen column density for the reference model (model CI_XN, see Sect. 2.5 and Table 5). All relevant disk model parameters are given in Table 4.

Based on the similarity solution for viscous accretion disks, we used an axisymmetric flared gas density structure with a Gaussian vertical profile and a powerlaw with a tapered outer edge for the radial column density profile (e.g. Lynden-Bell & Pringle 1974; Andrews et al. 2009). The vertical scale height as a function of radius is expressed by a simple powerlaw. The disk has a total mass of 0.01 M_⊙ and extends from 0.07 au (the dust sublimation radius) to 620 au where the total vertical hydrogen column density reaches N_{⟨ H ⟩ ,ver} ≈ 10²⁰ cm^-2 (panel (f) in Fig. 4).

For the dust density distribution we assumed a dust to gas mass ratio of d/g = 0.01. There is observational evidence for dust growth and settling in protoplanetary disks (e.g. Williams & Cieza 2011; Dullemond & Dominik 2004). To account for dust growth we assumed a power law dust size distribution f(a) = a^-3.5 with a minimum and maximum grain radius of a_min = 0.05 μm and a_max = 3000 μm. For dust settling we applied the method of Dubrulle et al. (1995) with a turbulent mixing parameter of 10^-2.

The irradiation of the disk by the star is important for the temperature and the chemical composition of the disk. For the photospheric emission of the star we use PHOENIX stellar atmosphere models (Brott & Hauschildt 2005). We considered a 0.7 M_⊙ star with an effective temperature of 4000 K and a luminosity of 1L_⊙. In addition to the photospheric emission, T Tauri stars commonly show far ultra-violet (FUV) excess (e.g. France et al. 2014) due to accretion shocks and strong X-ray emission (e.g. Güdel & Nazé 2009). For the excess FUV emission, we used a simple power law spectrum with a total integrated FUV luminosity of L_FUV = 0.01 L_∗, in the wavelength interval [91.2 nm, 250 nm]. The details for the stellar X-ray properties are discussed in Sect. 2.3.1 above.

We use the molecules HCO⁺ and N₂H⁺ primarily to study the impact of SP ionization. To verify if our model gives reasonable results concerning HCO⁺ and N₂H⁺ abundances, we compared the modelled fluxes for the J = 3 − 2 transition of HCO⁺ and N₂H⁺ with the observational sample of Öberg et al. (2010, 2011a) finding a good agreement (for details see Appendix E).

2.5. Model series

It is likely that the different ionization sources are correlated. As already discussed the SP flux of young stars is actually derived from their stellar X-ray properties (Lee et al. 1998b; Feigelson et al. 2002). Also the CR ionization rate might be anti-correlated with the activity of the star (Cleeves et al. 2013). However, these possible correlations are not well understood. We therefore ran a series of full disk models using the X-ray, SP and CR spectra (described above) as inputs and also included the Turner SP model (Appendix D). We do not discuss any model with the active Sun SP spectrum as in this case SPs do not have a significant impact on the disk chemical structure (see Sect. 3.1). An overview of all presented models is given in Table 5.

3. Results

3.1. Ionization rates as a function of column density

Before we discuss our results for the full disk model, we compare the SP, X-ray and CR ionization rates as a function of the total hydrogen column density N_{⟨ H ⟩} (N_{⟨ H ⟩} = N_H + 2N_H2). Figure 5 shows such a comparison for our different input spectra discussed in Sects. 2.2 and 2.3.

From visual inspection of Fig. 5 it becomes clear that for a SP flux comparable to our Sun (active Sun spectrum) SP ionization cannot compete with X-ray ionization assuming typical T Tauri X-ray luminosities. However, for the active T Tauri SP spectrum SP ionization becomes comparable to X-ray ionization or even dominates for N_{⟨ H ⟩} ≲ 10²⁴ − 10²⁵ cm^-2. For N_H ≲ 10²³ cm^-2ζ_SP is determined by the particles with E_p ≲ 5 × 10⁷ eV whereas higher energy particles dominate for N_{⟨ H ⟩} ≳ 10²³ cm^-2. The kink at N_{⟨ H ⟩} ≈ 2 × 10²⁵ cm^-2 is caused by the rapid attenuation of the SPs at high column densities. At such high column densities even the most energetic particles have lost most of their energy and the ionization rate drops exponentially.

For X-rays, Fig. 5 shows the differences between the normal and high X-ray spectrum. The X-ray ionization rates are higher for the high X-ray spectrum due to the higher X-ray luminosity. Additionally, the harder X-ray photons can penetrate to deeper layers but are also more efficiently scattered than lower energy X-ray photons. Compared to the normal X-ray case, the X-ray ionization rate increases by several orders of magnitude for N_{⟨ H ⟩} ≳ 10²⁴ cm^-2 for the high X-ray case.

Galactic cosmic rays are the most energetic ionization source. The peak in the particle energy distribution is around 10⁸ − 10⁹ eV (see Fig. 2). As a consequence the CR ionization rate ζ_CR stays mostly constant and only decrease for N_{⟨ H ⟩} ≳ 10²⁵ cm^-2. Only for such high column densities CR particle absorption becomes efficient.

In the Turner model it is implicitly assumed that SPs have the same energy distribution than Galactic CRs (see Appendix D). As a consequence the SP ionization rate in the Turner model is simply a scaled up version of the CR ionization rate. The slight differences to our model in CR attenuation is caused by the different methods used to calculate the SP/CR ionization rates; Turner & Drake (2009) use the fitting formulae of Umebayashi & Nakano (2009). Compared to our active T Tauri SP spectrum ζ_SP in the Turner model is larger for N_{⟨ H ⟩}> 10²⁵ cm^-2 but significantly lower at low column densities.

3.2. Disk ionization rates

In Fig. 6 we show the ionization rates as a function of the vertical hydrogen column density N_{⟨ H ⟩ ,ver} at two different radii of the disk. Shown on the figure are the models CI_XN_SP (ISM CR, normal X-rays; left panel) and CL_XH_SP (low CR, high X-rays; right panel).

CRs are only significantly attenuated for N_{⟨ H ⟩ ,ver}> 10²⁵ cm^-2 and r ≲ 1 au. For most of the disk, CRs provide a nearly constant ionization rate of ζ_CR ≈ 2 × 10^-17 s^-1 for the ISM like and ζ_CR ≈ 2 × 10^-19 s^-1 for the low CR spectrum.

X-rays are strongly attenuated as a function of height and radius (i.e. geometric dilution). However, due to scattering X-rays can become the dominant midplane ionization source for large regions of the disk. For ISM-like CRs, CR ionization is the dominant midplane ionization source even in the high X-ray models. In the low CRs models X-rays are the dominant midplane ionization source for r ≲ 100 au in the normal X-ray model and for all radii in the high X-ray model.

Differently to X-rays, SPs are not scattered towards the midplane. Due to their high energies they propagate along straight lines (provided that the SPs are not shielded by magnetic fields, see Sect. 4.3.2). As SPs are of stellar origin they penetrate the disk only along radial rays. The radial column densities close to the midplane of the disk are N_{⟨ H ⟩ ,rad} ≫ 10²⁵ cm^-2 and therefore SPs are already strongly attenuated at the inner rim of the disk. From Fig. 6 we see that the SP ionization rate ζ_SP drops below 10^-19 s^-1 for N_{⟨ H ⟩ ,ver} ≳ 10²⁴ cm^-2 at r = 1 au and for N_{⟨ H ⟩ ,ver} ≳ 10²² cm^-2 at r = 100 au. However, at higher layers (N_{⟨ H ⟩ ,ver} ≲ 10²² − 10²³ cm^-2) SPs are the dominant ionization source even in the high X-ray models. Expressed in radial column densities: SPs are the dominant ionization source in disk regions with N_{⟨ H ⟩ ,rad} ≲ 10²⁴ − 10²⁵ cm^-2.

In Fig. 7 we show the dominant ionization source at every point in the disk for all SP models. An ionization source is dominant at a certain point in the disk if its value is higher than the sum of the two other ionization sources. The first two columns in Fig. 7 show our models, the last column shows the Turner models. In our models, SPs are the dominant ionization source in the upper layers of the disk (above the white solid contour line for N_{⟨ H ⟩ ,rad} = 10²⁵ cm^-2), whereas in the midplane always CRs or X-rays dominate.

For the Turner model the picture is quite different. In their model SPs can also penetrate the disk vertically (Appendix D) and reach higher vertical column densities before they are completely attenuated (Fig. 5). As a consequence SPs can become the dominant ionization source in the midplane of the disk (e.g. for the low CR case). In the upper layers always X-rays dominate as ζ_SP<ζ_X for low column densities. In the Turner model ζ_SP ≲ 10^-13 s^-1 for N_{⟨ H ⟩ ,rad}< 10²⁵ cm^-2 which is several orders of magnitudes lower than in our models. The reason for this is that in the Turner model SPs are simply a scaled version of ISM like CRs. The high ζ_SP values in our model in the upper layers of the disk are caused by the high number of particles with energies E_p ≲ 10⁸ eV, which are missing in the Turner model.

3.3. Impact on HCO⁺ and N₂H⁺

The molecules HCO⁺ and N₂H⁺ are the two most observed molecular ions in disks (e.g. Dutrey et al. 2014) and are commonly used to trace the ionization structure of disks (e.g. Dutrey et al. 2007; Öberg et al. 2011b; Cleeves et al. 2015). Also Ceccarelli et al. (2014) used these two molecules to trace SPs in a protostellar envelope.

The main formation path of HCO⁺ and N₂H⁺ is the ion-neutral reaction of H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$ with their parent molecules CO and N₂, respectively. H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$ is created by ionization of H₂ by CRs, X-rays and in our model additionally by SPs. The main destruction pathway for HCO⁺ and N₂H⁺ is via dissociative recombination with free electrons.

The chemistry of HCO⁺ and N₂H⁺ is linked to the CO freeze-out. To form HCO⁺, gas phase CO is required, whereas N₂H⁺ is efficiently destroyed by CO (e.g. Aikawa et al. 2015). Consequently the N₂H⁺ abundance peaks in regions where CO is depleted and N₂, the precursor of N₂H⁺, is still in the gas phase. The result of this chemical interaction is a vertically layered chemical structure for HCO⁺ and N₂H⁺ (see Fig. 8). For further details on the HCO⁺ and N₂H⁺ chemistry see Appendix B, where we also list the main formation and destruction pathways for HCO⁺ and N₂H⁺ (Table B.1).

3.3.1. Abundance structure

In the following we describe details of the molecular abundance structure that are relevant for the presentation of our results for our reference model CI_XN. The abundance ϵ of a molecule X is given by ϵ(X) = n_X/n_{⟨ H ⟩}, where n_X is the number density of the respective molecule and n_{⟨ H ⟩} = n_H + 2 n_H2 is the total hydrogen number density. Figure 8 shows the resulting abundance structure for CO, HCO⁺ and N₂H⁺ for the CI_XN model.

We define the location of the CO ice line where the CO gas phase abundance is equal to the CO ice-phase abundance (white dashed line in Fig. 8). The CO ice line is located at dust temperatures in the range T_d ≈ 23 − 32 K (density dependence of the adsorption/desorption equilibrium; e.g. Furuya & Aikawa 2014). The radial CO ice line in the midplane (z = 0 au) is at r ≈ 12 au and T_d ≈ 32 K. At r ≈ 50 au the vertical CO ice line is at z ≈ 8.5 au (z/r ≈ 0.17) and T_d ≈ 26 K. Inside/above the CO ice line ϵ(CO) ≈ 10^-4. Outside/below the CO ice line ϵ(CO) rapidly drops to values ≲10^-6. In regions where non-thermal desorption processes are efficient (r ≳ 150 au) ϵ(CO) ≈ 10^-6 down to the midplane. For the regions inside/above and outside/below the CO ice line we use the terms warm and cold molecular layer, respectively.

There are two main reservoirs for HCO⁺, one in the warm molecular layer above the CO ice line and one in the outer disk (r ≳ 150 au) below the CO ice line where non-thermal desorption becomes efficient. In the warm molecular layer, the ionization fraction ϵ(e⁻) ≈ 10^-7 is dominated by sulphur as it is ionized by UV radiation (e.g. Teague et al. 2015, see also Sect. 4.2.3). Those free electrons efficiently destroy molecular ions via dissociative recombination. This causes a dip in the vertical HCO⁺ abundance structure within the warm molecular layer with ϵ(HCO⁺) ≈ 10^-12 − 10^-11, whereas at the top and the bottom of the warm molecular layer ϵ(HCO⁺) reaches values of ≈10^-10 − 10^-9. The peak in the top layer is mainly caused by the high X-ray ionization rate for H₂ (ζ_X ≳ 10^-12 s^-1). At the bottom of the warm molecular layer more HCO⁺ survives. This region is already sufficiently shielded from UV radiation and the free electron abundance drops rapidly. In the second reservoir, below the CO ice line where non-thermal desorption is efficient ϵ(HCO⁺) ≈ 10^-11 − 10^-10.

The main N₂H⁺ reservoir resides in the cold molecular layer just below the CO ice line with ϵ(N₂H⁺) ≳ 10^-11. The lower boundary of this layer with ϵ(N₂H⁺) < 10^-11 is reached at T_d ≈ 16 K where ϵ(N₂) ≲ 10^-6 due to freeze-out. Radially this layer extends from the inner midplane CO ice line out to r ≈ 250 − 300 au. Close to the midplane ϵ(N₂H⁺) ≲ 10^-12 for r ≳ 150 au due to non-thermal desorption of ices. There is also a thin N₂H⁺ layer at the top of the warm molecular layer with ϵ(N₂H⁺) ≈ 10^-12 extending from the inner radius of the disk out to r ≈ 100 au. In this layer the X-ray ionization rate is high enough to compensate for the destruction of N₂H⁺ by CO.

The detailed appearance of this layered structure is especially sensitive to the dust temperature and therefore also to dust properties (e.g. dust size distribution). The above described abundance structure for CO, HCO⁺ and N₂H⁺ is consistent with the model of Aikawa et al. (2015) that includes millimetre sized dust particles with a dust size distribution similar to what is used here (for details see Appendix B).

3.3.2. Vertical column densities

To study the impact of SP ionization quantitatively we compare vertical column densities of HCO⁺ and N₂H⁺ for models with and without SPs. In Fig. 9 we show the vertical column densities N_HCO⁺ and N_N2H⁺ as a function of the disk radius r for all models listed in Table 5. The left column in Fig. 9 shows the models with ISM CRs, the right column the models with low CRs. At first we discuss the models without SPs and compare them to other theoretical models.

The N_HCO⁺ profile shows a dip around r ≈ 50 − 100 au in the ISM CR models CI_XN and CI_XH (high X-rays). This dip is also seen in the models of Cleeves et al. (2014). Unlike Cleeves et al. (2014), in our model this dip is not predominantly due to the erosion of CO by reactions with He⁺ (“sink effect” e.g. Aikawa et al. 1997, 2015; Bergin et al. 2014; Furuya & Aikawa 2014), but mainly due to the interplay of CO freeze-out and non-thermal desorption in the outer disk. The CO sink effect is also active in our model but less efficient (see Appendix B.2.2).

The lack of HCO⁺ in the disk midplane at r ≈ 50 − 100 au due to CO freeze-out is also visible in the HCO⁺ abundance structure shown in Fig. 8. Non-thermal desorption in the midplane produces CO abundances ≳10^-7 for r ≳ 150 au and consequently also a slight increase in N_HCO⁺. In the low CR models the ionization rate is too low to produce a significant amount of HCO⁺ in the cold molecular layer and the dip in the profile vanishes.

N₂H⁺ traces the distribution of gas phase CO as it is efficiently destroyed by CO (Qi et al. 2013a,b, 2015). This is also seen in our model. The sharp transition in the N₂H⁺ column density at r ≈ 30 au traces the onset of CO freeze-out (Fig. 9). We note, however, that the actual midplane CO ice-line is at ≈ 12 au (see Appendix C for details). In the ISM CR models N_N2H⁺ is dominated by CR ionization as N₂H⁺ mainly resides in the cold molecular layer. In the warm N₂H⁺ layer X-ray ionization dominates. However, due to the lower densities in the warm molecular layer, this layer only contributes significantly to the column density within the radial CO ice line, and the impact of X-rays is only visible there (compare models CI_XH and CI_XN in Fig. 9). The high X-ray luminosity decreases the contrast between the peak of N_N2H⁺ close to the radial CO ice line and N_N2H⁺ inside the CO ice line by about a factor of five.

For the low CR case N_N2H⁺ drops by more than an order of magnitude compared to the ISM CR case. Such a strong impact of CR ionization on N_N2H⁺ is also reported by Aikawa et al. (2015) and Cleeves et al. (2014). Higher X-ray luminosities can compensate for low CR ionization only to some extent. In the high X-ray model CL_XH, N_{N2H +} is lower by a factor of five compared to the ISM CR models.

3.3.3. Impact of SPs

The column densities N_HCO⁺ and N_N2H⁺ for models with SPs are shown in Fig. 9 (blue solid and dashed lines). The solid orange line in Fig. 9 shows the results for the Turner model. We discuss the Turner model separately in Sect. 3.3.4. We define a change in the column densities by at least a factor of three compared to the reference model as significant. This is indicated by the grey area around the column density profiles of the reference models CI_XN and CL_XN.

To better quantify the impact of SPs compared to the competing H₂ ionization sources , X-rays and CRs, we introduce the weighted column density ${\mathit{N}}_{\mathit{\zeta }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathit{n}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{\times }{\mathit{f}}_{\mathit{\zeta }}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{z}\hspace{0.17em}\hspace{0.17em}\mathrm{\left[}\mathrm{c}{\mathrm{m}}^{-2}\mathrm{\right]}\mathit{.}$(5)N_ζ is the weighted column density for a particular H₂ ionization source ζ, n is the number density of a particular molecule in units of cm^-3 and ${\mathit{f}}_{\mathit{\zeta }}\mathrm{=}\frac{\mathit{\zeta }}{{\mathit{\zeta }}_{\mathrm{X}}\mathrm{+}{\mathit{\zeta }}_{\mathrm{SP}}\mathrm{+}{\mathit{\zeta }}_{\mathrm{CR}}}\mathrm{·}$(6)N_ζ represents the fraction of the column density dominated by a particular ionization source ζ.

In Fig. 10 we show N_ζ for HCO⁺ and N₂H⁺ normalized to the total column density of the respective molecule as a function of radius. The region where more than 50% of the total column density of the molecules arise from disk regions above the CO ice line (i.e. from the warm molecular layer) is roughly indicated on each plot.

SP ionization has a significant impact on the N_HCO⁺ profile in all our models. In the ISM CR model (CI_XN_SP) N_HCO⁺ increases by a factor ≈3 for 50 ≲ r ≲ 100 au and the dip in the profile seen in the models without SPs (CI_XN) vanishes (top left panel in Fig. 9). In the low CR models this region increases to 25 ≲ r ≲ 150 au and N_HCO⁺ reaches values up to an order of magnitude higher compared to the CL_XN model (top right panel in Fig. 9). In the models with high X-rays, SP ionization still has an significant impact and N_HCO⁺ increases by up to a factor of three (dashed lines in Fig. 9).

This situation is also clearly visible in Fig. 10 where we show N_ζ as a function of radius (Eq. (5)). Although X-rays can be the dominant ionization source close to the star, in all models N_HCO⁺ is dominated by SP ionization for 50 ≲ r ≲ 100 − 200 au. Figure 10 also shows that in the low CR models N_HCO⁺ is mainly built up in the warm molecular layer for r ≲ 200 au. In this region (indicated in Fig. 10) the warm molecular layer contributes more to the total column density than the cold molecular layer.

For N₂H⁺ the picture is more complex. In the ISM CR models SPs have only very little impact on N_N2H⁺. Only in the inner 30 au, within the radial CO ice line, the N₂H⁺ profile is significantly affected. The reason for this is the high SP ionization rate ζ_SP ≳ 10^-12 s^-1 in the warm molecular layer of the disk close to the star. In this region the abundance ratio of HCO⁺/N₂H⁺ drops from >10³, in the models without SPs, to around 10 to 100 in models with SPs. These high ratios can be explained by the efficient destruction of N₂H⁺ by CO. However, due to the high ζ_SP this destruction path becomes less important, and the molecular ion abundance are mainly determined by the balance between ionization and recombination.

Ceccarelli et al. (2014) reported a very low measured HCO⁺/N₂H⁺ ratio between three and four in the Class 0 source OMC-2 FIR 4. They explain this low ratio by the high ionization rates due to SPs (ζ_SP> 10^-14). In our disk model ζ_SP in the warm molecular layer is comparable, but the HCO⁺/N₂H⁺ ratio is much larger than four. This higher ratio is due to the higher densities of 10⁸ − 10¹⁰ cm^-3 and the stronger UV field in the warm molecular layer, compared to the physical conditions in OMC-2 FIR 4. This is in agreement with the chemical models presented in Ceccarelli et al. (2014).

The impact of SPs on the N₂H⁺ abundance in the warm layer can extend out to r ≈ 200 au (similar to HCO⁺). However, for r> 30 au this layer does not significantly contribute to the total N₂H⁺ column density as N_N2H⁺ is dominated by the high density layer below the vertical CO ice line.

Beyond the radial CO ice line N_N2H⁺ is dominated by CRs in the ISM CR models and is not affected by X-rays nor SPs (bottom left panel in Fig. 9). Compared to the ISM CR models, N_N2H⁺ is reduced by a factor of a few around the peak and by more than an order of magnitude at larger radii in the low CR models. As a consequence the profile is also steeper. Although higher X-rays (CL_XH model) and also SPs can to some extent compensate low CR ionization rates, the N_N2H⁺ profile is still steeper and N_N2H⁺ is lower by a factor of between approximately two and six than in the ISM CR models for r ≳ 30 au.

From Fig. 10 we see that only in the CL_XN_SP model SPs dominate N_N2H⁺ for r ≳ 70 au. Actually the importance of SPs increases with r in this model. In the CL_XN_SP model SPs are the dominant ionization source in regions with N_{⟨ H ⟩ ,rad} ≲ 10²⁵ cm^-2. The N_{⟨ H ⟩ ,rad} = 10²⁵ cm^-2 iso-contour is below the vertical CO ice line for r> 70 au and the layer between the CO ice line and N_{⟨ H ⟩ ,rad} = 10²⁵ cm^-2 becomes thicker with radius (see Figs. 7 and 8). This also explains the change in the slope of N_N2H⁺ compared to models without SPs.

In the high X-ray model, CL_XH_SP, the picture is quite different. N_N2H⁺ is now dominated by X-rays for r ≳ 30 au. X-rays are efficiently scattered towards the midplane and therefore ζ_X>ζ_SP in the cold N₂H⁺ layer (see also Fig. 7). X-rays affect the cold N₂H⁺ layer at all radii therefore the slope of N_N2H⁺ is steeper than in the model where SP dominate (compare the blue solid line with the red dashed line in Fig. 9).

Our results show that HCO⁺ is always significantly affected by SP ionization but N₂H⁺ only in models with low CRs and normal X-rays (for r> 30 au). As SPs can only reach the upper layers of the cold molecular layer, N₂H⁺ is less sensitive to SP ionization than HCO⁺.

3.3.4. Impact of SPs in the Turner model

In the Turner model SP ionization is just a scaled up version of CR ionization where SPs can also penetrate the disk vertically (see Appendix D). The results of the Turner models are also shown in Figs. 9 and 10.

In the ISM CR models there is no significant impact on HCO⁺ and N₂H⁺ by SP ionization. ζ_SP is significantly lower at low column densities compared to our models (see Sect. 3.1). Therefore X-rays are the dominant ionization source in the warm molecular layer and HCO⁺ is not significantly affected by SP ionization. The slight increase in the HCO⁺ column density is mainly due to the higher X-ray luminosity in the Turner model compared to our reference model with normal X-rays. Similar to our models CRs dominate in the cold molecular layer.

In the Turner model with low CRs, SPs become the dominant ionization source in the cold molecular layer as they also penetrate the disk vertically. In this layer ζ_SP reaches values of ≈10^-17 s^-1 at r ≈ 100 au. However, in the Turner model SPs also cannot compensate for a low CR ionization rate as ζ_SP ∝ 1 /r² (geometric dilution).

The impact of SPs in the Turner model is rather limited and restricted to the cold molecular layer, in strong contrast to our models. The differences are mainly due to the assumptions concerning the SP transport. In the Turner model SPs hit the surface of the disk and penetrate the disk vertically, whereas we assume that SPs travel only along radial rays. However, both approaches are an approximation of a likely more complex picture of SP transport in disks. We discuss this in more detail in Sect. 4.3.

4. Discussion

4.1. Constraining the SP flux of T Tauri stars

SP ionization has a significant impact on the column densities of HCO⁺ and N₂H⁺ in all our models. To actually constrain the SP flux from observations it is necessary to disentangle the contribution of SP ionization from the competing ionization sources CRs and X-rays.

From a chemical point of view all three ionization sources act the same way, they ionize molecular hydrogen and drive the molecular ion chemistry. However, they also show distinct differences in how they irradiate and penetrate the disk. CRs act like a background source and irradiate the disk isotropically, whereas X-rays and SPs originate from the star and act like a point source. A further difference is their energy distribution. Because of their high energies, CRs and SPs tend to move on straight lines whereas (hard) X-rays also experience scattering during their interaction with the disk.

Those differences in their irradiation properties and their energy distribution allow to disentangle their impact on the ion chemistry at different locations in the disk. Our models show that the stellar ionization sources are more effective closer to the star and at the surface layers of the disk (Fig. 10). X-rays can also become an important ionization source in the midplane of the disk but not for the whole disk as, roughly speaking, the ionization rate of stellar ionization sources is ∝1 /r². CR ionization affects the whole disk, but is, in contrast to the stellar ionization sources, more important for the outer disk and the cold layers and the midplane of the disk (this argument also holds for SLR ionization).

So far we have only shown models using the commonly proposed SP flux of f_p(E_p ≥ 10 MeV) ≈ 10⁷ protons cm^-2 s^-1 (e.g. Feigelson et al. 2002). However, this value should be only seen as an order of magnitude estimate (see Sect. 2.2).

In Fig. 11 we show the column densities for HCO⁺ and N₂H⁺ for models with a SP flux a factor of ten higher and lower with respect to the reference value. Also shown are the reference model for the low CR case CL_XN and the model with the reference SP flux CL_XN_SP. For these models we use the low CR ionization rate (ζ_CR ≈ 2 × 10^-19 s^-1) and the normal X-ray luminosity (L_X = 10³⁰ erg s^-1). Further we show a model with L_X = 3 × 10³¹ erg s^-1, to illustrated the impact of (very) strong X-ray emission (e.g. ζ_X ≳ 10^-17 s^-1 in the disk midplane).

In the low SP flux model (CL_XN_SPL) the impact on HCO⁺ and N₂H⁺ is quite limited. For such a case it is still possible to define upper limits for the SP flux in disks.

In the high SP flux model (CL_XN_SPH) the SP ionization rate reaches values of ζ_SP ≳ 10^-13 s^-1 and ζ_SP ≈ 10^-17 s^-1 in the warm and cold molecular layer respectively. The N₂H⁺ column density profile beyond the radial CO ice line is comparable to the profile for models with ISM CRs (compare with Fig. 9). Also the model with very high X-rays (CL_X31) and no SPs shows a similar profile for N₂H⁺. However, as seen from Fig. 11 the corresponding HCO⁺ profiles differ significantly.

Comparing the very high X-ray model (CL_X31) to the reference model with SPs (CL_XN_SP) shows that the corresponding HCO⁺ profiles are similar but the N₂H⁺ profiles differ by an order of magnitude. This shows again that it is indeed possible to distinguish between the different ionization sources by simultaneous modelling of HCO⁺ and N₂H⁺ column density profiles.

To trace this interplay of ionization sources, spatially resolved observations of molecular ion lines tracing different vertical layers of the disk are required. With modern (sub)millimetre interferometers like the Atacama Large Millimeter Array (ALMA), the NOrthern Extended Millimeter Array (NOEMA) and the Submillimeter Array (SMA) such observations with a spatial resolution of tens of au are already possible (e.g. Qi et al. 2013b; Cleeves et al. 2015; ALMA Partnership et al. 2015; Yen et al. 2016) and will become available on a regular basis in the near future. Here, we use HCO⁺ and N₂H⁺ as the tracers of the warm and cold molecular layer respectively, but also other molecules like DCO⁺, which traces similar regions as N₂H⁺ (Teague et al. 2015; Mathews et al. 2013) can be used.

Complementary, far-infrared lines of HCO⁺ and N₂H⁺, as used by Ceccarelli et al. (2014) to trace SP ionization in a protostellar envelope, are good tracers of molecular ion emission in the warm inner region of the disk. However, a more detailed analysis with proper modelling of line emission is required to identify the best observational tracers of SP ionization. We will present such an analysis in a follow-up paper.

4.2. Chemical implications

Besides the H₂ ionization rates there are other “chemical parameters”, which have an impact on the molecular ion abundance in disks. In the following we discuss the dependence of our results on the location of the CO ice line, depletion of CO and the assumed initial metal abundances. Those chemical properties of disks are not well constrained from observations and/or can vary between different targets.

4.2.1. Location of the CO ice line

Recent ALMA observations provide direct constrains on the location of the CO ice line. However, these results depend on the method or more precisely the molecule used to trace the CO ice line (see Qi et al. 2013b; Schwarz et al. 2016; Nomura et al. 2016, for TW Hya). Further, due to complex chemical processes like the CO sink effect it is possible that the actual location of the CO ice line does not only depend on the CO freeze-out temperature (Aikawa et al. 2015, Sect. 3.3.2).

To investigate the dependence of our results on the location of the CO ice line we artificially moved the CO ice line in our model by adapting the binding energy for CO. We consider two cases: E_B(CO) = 950 K and E_B(CO) = 1350 K (i.e 200 K lower and higher compared to our reference model). In both cases we kept the ratio of E_B(N₂) /E_B(CO) = 0.67 constant (see Appendix B.2.2). As a consequence also the N₂H⁺ layer moves accordingly to the CO ice line (see Sect. 3.3.1). For E_B(CO) = 950 K the CO ice line moves to T_d ≈ 20 − 24 K, (i.e. deeper into the disk) and for E_B(CO) = 1350 K to T_d ≈ 25 − 36 K (i.e. higher up in the disk).

For a CO ice line deeper in the disk the contribution of CR ionization to the total column of HCO⁺ and N₂H⁺ increases in the ISM CR models. In the ISM CR models SP ionization is not significant anymore (i.e. N_HCO⁺ increases by less than a factor two). However, in the low CR models the impact of SPs remains significant.

A CO ice line higher up in the disk has the opposite effect. The total column densities of the molecular ions are now dominated by layers higher up in the disk which can efficiently be ionized by SPs. As a consequence the relative contribution of SP ionization to N_HCO⁺ and N_N2H⁺ increases.

In summary for a CO ice line location deeper in the disk SP ionization becomes less important, for a CO ice line higher up in the disk SP ionization becomes more important. However, in both cases the interplay of the different ionization sources is qualitatively speaking similar to what is shown in Fig. 10.

4.2.2. CO depletion

There is observational evidence for CO depletion in protoplanetary disks (Dutrey et al. 1997; Bruderer et al. 2012; Favre et al. 2013; Kama et al. 2016a; Schwarz et al. 2016; McClure et al. 2016). The best constraint case is TW Hya. Using spatially resolved ALMA spectral line observations of several CO isotopologues Schwarz et al. (2016) derived a uniform CO abundance of ≈10^-6 in the warm molecular layer, two order of magnitudes lower than the canonical value of ≈10^-4. However, the degree of CO depletion seems to vary from source to source. Using Herschel HD J = 1 − 0 line observations McClure et al. (2016) derived CO depletions of a factor of approximately five and up to ≈100 for DM Tau and GM Aur, respectively.

The cause of CO depletion in disks is not yet clear. Although freeze-out of CO certainly contributes to depletion it is unlikely that it is the only process acting. Several other mechanisms that can at least partly explain CO depletion are proposed:

• the destruction of CO by He⁺ and the subsequent conversion of atomic carbon to more complex carbon bearing molecules with higher freeze-out temperatures (Aikawa et al. 1996; Bergin et al. 2014; Helling et al. 2014; Furuya & Aikawa 2014);

• depletion of CO in layers above the CO ice line (up to T ≈ 30 K) due to conversion of CO to CO₂ on the surfaces of dust grains (Reboussin et al. 2015);

• CO isotopologue selective photodissociation, which affects CO isotopologue line emission and therefore the derived CO depletion factors;

• carbon and/or oxygen depletion in the warm disk atmosphere due to settling and mixing of ice coated dust grains (Du et al. 2015; Kama et al. 2016b).

It is as yet unclear which of the proposed mechanisms is the most efficient one; none of them can be excluded with certainty. It is also possible that all of these processes are at work. So far the impact of CO (carbon and/or oxygen depletion) on molecular ion emission was not studied in detail. For the modelling of HCO⁺ and N₂H⁺ line emission of TW Hya, Cleeves et al. (2015) reduced the initial atomic carbon abundance by two orders of magnitude to match C¹⁸O line observation. However, the impact of C and CO depletion on HCO⁺ and N₂H⁺ was not discussed in detail.

To simulate CO depletion we simply reduce the total carbon and oxygen element abundances by one order of magnitude throughout the disk. This results in a CO abundance of ≈10^-5 in the warm molecular layer. We applied this “artificial” CO depletion to all models listed in Table 5; all other parameters of the models are fixed.

The CO depletion models show a factor of ≈5 − 10 lower N_HCO⁺ for r ≳ 50 au compared to the non depleted models. For r< 30 auN_N2H⁺ increases by more than an order of magnitude. N_N2H⁺ beyond the CO ice line is not affected as N_N2H⁺ resides within the CO freeze-out zone where gas phase CO is anyway depleted.

In the CO depletion models SP ionization is slightly more efficient for N_N2H⁺ as the contribution of the warm N₂H⁺ layer to N_N2H⁺ increases. Due to the lower CO abundance in the warm molecular layer the N₂H⁺ abundance increases as the destruction pathway via CO is less efficient. For HCO⁺ the opposite is true. The HCO⁺ abundance decreases by roughly an order of magnitude in the warm molecular layer, whereas in the CO freeze-out zone the impact is smaller (i.e. CO is frozen-out anyway). Relatively speaking the contribution of the cold HCO⁺ to N_HCO⁺ increases in the CO depletion models. Therefore the impact of SPs on N_HCO⁺ is less significant whereas X-rays and CRs become more important. Although there are some differences, the impact of SPs on N_N2H⁺ and N_HCO⁺ is qualitatively very similar to the non-depleted models. In particular the main trends derived from Fig. 10 are also seen in the CO depletion models.

CO depletion is certainly more complex than modelled here. A more thorough study of the impact of CO depletion on molecular ion abundances is certainly desirable and possibly provides new constraints on CO gas phase depletion in disks. However, this is beyond the scope of this paper.

4.2.3. Metal abundances

Heavy metals such as sulphur play an important role in the molecular ion disk chemistry (Teague et al. 2015, Rab et al., in prep.). Here we refer with the term metals to the elements Na, Mg, Si, S and Fe. Metal ionization due to UV radiation can produce a large number of free electrons. Those free electrons destroy molecular ions like HCO⁺ and N₂H⁺ via dissociative recombination. Dissociative recombination is more efficient than radiative recombination of metals. As a consequence a high abundance of metals significantly reduces the abundance of molecular ions (e.g. Mitchell et al. 1978; Graedel et al. 1982).

In the ISM and in disks most of the metals are likely locked up in refractory grains and are therefore depleted compared to Solar abundances. We use metal abundances similar to the commonly used “low metal” abundances (Graedel et al. 1982; Lee et al. 1998a). These low metal abundances are depleted by a factor of ≈100 − 1000 compared to Solar abundances (i.e. low metal sulphur abundance ϵ(S) ≈ 10^-7). However, the actual gas phase abundance of metals is difficult to constrain from observations and a stronger degree of depletion already prior to disk formation is possible (e.g. Maret & Bergin 2007; Maret et al. 2013).

To investigate the dependence of our results on the metal abundances we deplete the initial gas phase metal abundance by an additional factor of ten compared to the low metal abundances (i.e. ϵ(S) ≈ 10^-8). This means that a larger fraction of metals is locked-up in refractory dust grains and cannot be released back into the gas phase.

Decreasing the metal abundances increases ϵ(HCO⁺) in the warm molecular layer by up to an order of magnitude (i.e. the gap in the vertical ϵ(HCO⁺) profile nearly vanishes; Sect. 3.3.1). In the cold molecular layer the metal abundances are not as important since the metals are frozen-out on dust grains anyway. However, in regions where non-thermal desorption processes are efficient (i.e. metal ices are released back into the gas-phase), ϵ(HCO⁺) but also ϵ(N₂H⁺) are higher by a factor of a few in the lower metal abundance models. The HCO⁺ column density is higher by a factor of between approximately two and three at all radii in the strong metal depletion models compared to the models with the reference abundances. The N₂H⁺ column density is only affected for r ≳ 150 au (higher by a factor of between approximately two and three), where non-thermal desorption of metals is efficient.

Lower gas phase metal abundances lead to an increase of the molecular ion abundances in the warm molecular layer and the contribution of this layer to the total column densities increases. However, in the strong metal depletion models the interplay and the relative contributions to the column densities of the different ionization sources is nearly identical to the reference model grid (i.e. Fig. 10 does not change significantly). Our arguments concerning the impact of SPs on HCO⁺ and N₂H⁺ are therefore also valid for the case of strong metal depletion.

4.3. Future prospects

Our model results show that SPs can indeed become an important ionization source in T Tauri disks. However, our model should be seen as a first attempt at a comprehensive modelling of SP ionization in protoplanetary disks. In the following we discuss further important aspects like variability, non-stellar origin of SPs, the importance of magnetic fields and future prospects for SP modelling in protoplanetary disks.

4.3.1. Flares and variability

In our models we assumed continuous (i.e. time-averaged) particle fluxes and X-ray luminosities. Although this is a reasonable assumption (see Sect. 2.2) it is likely that the disk is also hit by singular powerful X-ray and/or SP flares.

From the X-ray COUP survey of the ONC cloud, Feigelson et al. (2002) and Wolk et al. (2005) derived a median X-ray flare luminosity of ≈6 × 10³⁰ erg s^-1, comparable to the X-ray luminosity in our high X-ray models, and peak flare luminosities up to 10³¹ erg s^-1 . The duration of such flares can last from hours up to three days with a typical frequency of roughly one powerful flare per week.

Ilgner & Nelson (2006) argue that for such flare properties the disk ion chemistry responds to time-averaged ionization rates. However, in case the duration between two flares is longer than the recombination time scale in the disk, such singular and strong flares would produce asymmetric features in molecular ion emission (i.e. a singular flare only affects a certain fraction of the disk). The spatially resolved HCO⁺J = 3 − 2 SMA observations of TW Hya show indeed such an asymmetric structure (Cleeves et al. 2015). However, as discussed by Cleeves et al. (2015) these features could also have a different origin like spiral arms or an hidden planet locally heating the disk.

If such asymmetric structures are caused by stellar flares they would provide complementary constraints on the X-ray/SP activity of the star. Multi-epoch data of spatially resolved molecular ion emission is required to prove the flare scenario (i.e. the features should disappear quickly). The current (sub)mm interferometers like ALMA, NOEMA and SMA provide the required spatial resolution and might even allow for monitoring of disks in molecular ion lines on a daily/weekly basis in the future.

Modelling of such observations does not necessarily require 3D chemical models. Assuming that the disk physical structure is not affected by flares, radial cuts through the disk can be modelled with 2D (time-dependent) thermo-chemical models as presented here.

4.3.2. Magnetic fields

In our model we neglect the impact of magnetic fields on the SP transport. Stellar and disk magnetic fields are in particular relevant for the question if particles actually hit the disk (see Feigelson et al. 2002, for a discussion).

Magnetic fields can either drag the particles away from the disk (e.g. like in a wind) but could also funnel the particles and concentrate the ionizing flux in particular regions of the disk. In the first scenario SPs can still have an impact on the upper layers of the disk but certainly not on the disk midplane. In the second scenario particles are likely focused on regions close to the inner radius of the disk and their impact on the outer disk will become less significant. Also the trajectory of the particles will be affected, and they might penetrate the disk also vertically. Our models, where particles are transported only radially, are closer to the wind scenario whereas the Turner model would represent an extreme case for magnetically focused particles. However, to qualitatively estimate the impact of magnetic fields more complex SP transport models are required.

It is possible to consider magnetic field effects in high energy particle transport models (Desch et al. 2004; Padovani & Galli 2011; Padovani et al. 2013b). In principle such methods can also be applied to disk models. The main challenge though is to determine the structure of the star and disk magnetic fields. However, as argued by Ceccarelli et al. (2014) identifying distinct observational signatures of SP ionization in disks and/or envelopes would allow to derive constrains for the magnetic field structure. This is certainly challenging, but with the availability of spatially resolved observational data and interpretation of such data with (improved) thermo-chemical disk models this might be feasible in the near future.

4.3.3. Non-stellar origin of high energy particles

Besides the stellar surface the close environment of young stars offers also alternative particle acceleration sites. X-ray flares and particles can be produced close to the inner disk in the so called reconnection ring where the stellar and disk magnetic field interact (X-wind model Shu et al. 1994, 1997). More recently Padovani et al. (2015, 2016) proposed jet shocks as alternative acceleration sites for high energy particles.

In the X-wind model the particle source and also the X-ray emitting source is located closer to the disk and slightly above the disk midplane. The typically assumed source location in this “lamppost” scenario is r_L ≈ 0.05 au and z_L ≈ 0.05 au (≈10 R_⊙) (Lee et al. 1998b; Igea & Glassgold 1999). For X-rays Ercolano et al. (2009) found that the height of the emitting source has relatively little impact concerning X-ray radiative transfer in disks. Moving the emitting source closer to the disk certainly has an impact on the very inner disk but for e.g. 10 au distance from the star the stellar particle flux increases only by about 1% compared to the stellar origin. The height of the SP emitting source has some impact on where in the disk SPs see a column density N_{⟨ H ⟩} ≳ 10²⁵ cm^-2. For particles accelerated at a height z_L ≈ 0.05 au and moving along a ray parallel to the midplane this happens at r ≈ 0.5 au. The consequences are that SPs can penetrate into slightly deeper layers of the disks but they still cannot penetrate to the disk midplane. Therefore our conclusions on the impact of SPs on HCO⁺ and N₂H⁺ remain valid for the X-wind scenario.

In the jet shock scenario, the emitting source would be located far above the disk. Padovani et al. (2016) considered a particle emitting source located at 1.8 × 10³ au above the star and calculated the resulting ionization rates for a 2D disk structure. Depending on their particle acceleration model they found ionization rates up to ζ ≈ 10^-14 s^-1 at the surface layers of the disk. Although X-ray ionization rates in the upper layer of the disk are typically higher such an irradiation scenario could have a significant impact on the ionization of the outer disk (r ≳ 50 au). In the outer disk, the jet accelerated particles can penetrate the disk vertically and therefore can also reach the disk midplane (similar to the Turner model, Sect. 3.3.4). This might become important if Galactic CRs are efficiently attenuated (i.e. in the low CR case).

The scenarios described above for non-stellar particle sources are certainly worth being investigated in detail. For the future we plan to extend our model to a proper treatment of such non-stellar emitting sources.

5. Summary and conclusions

In this work we investigated the impact of stellar energetic particle (SP) ionization on disk chemistry with a focus on the common disk ionization tracers HCO⁺ and N₂H⁺. We assumed a typical SP flux of f_p(E_p ≈ 10 MeV) ≈ 10⁷ protons cm^-2 s^-1 (at 1 au) as commonly proposed for T Tauri stars and a particle energy distribution derived from measurements of solar particle events. Based on a detailed particle transport model we derived an easy to use formula (see Sect. 2.2) to calculate the SP ionization rate in the disk as a function of hydrogen column density and radius, assuming that the particles can penetrate the disk only radially. With a small grid of models considering varying properties of the competing high energy disk ionization sources, X-rays and Galactic cosmic rays, we studied the interplay of the different ionization sources and identified possible observational tracers of SP ionization. Our main conclusions are the following:

• SPs cannot penetrate the disk mid-plane. At hydrogen column densitiesN_{⟨ H ⟩} ≳ 10²⁵ cm^-2 even the most energetic particles are attenuated (stopped) and the SP ionization rate drops rapidly. As the radial hydrogen column densities for full T Tauri disks are typically N_{⟨ H ⟩} ≫ 10²⁵ cm^-2 the midplane SP ionization rate is ζ_SP ≪ 10^-20 s^-1 already at a distance of 1 au from the star.

• For the assumed SP flux (see above), SPs become the dominant H₂ ionization source in the warm molecular layer of the disk above the CO ice line, provided that SPs are not shielded by magnetic fields. This is even true for enhanced X-ray luminosities (i.e. L_X = 5 × 10³⁰ erg s^-1).

• SP ionization can increase the HCO⁺ and N₂H⁺ column densities by factors of between approximately three and ten for disk radii r ≲ 200 au. The impact is more significant in models with low CR ionization rates (i.e. ζ_CR ≈ 10^-19 s^-1).

• SP ionization becomes insignificant for an SP flux one order of magnitude lower than the proposed value for T Tauri stars. In such a case H₂ ionization is solely dominated by X-rays and CRs.

• As SPs cannot penetrate the deep layers of the disk, X-rays and/or CRs usually remain the dominant H₂ ionization source in the cold disk layers (i.e. below the CO ice line). Therefore HCO⁺, which traces the warm molecular layer, is more sensitive to SP ionization than N₂H⁺ that resides in the cold molecular layer.

• Simultaneous modelling of spatially resolved radial intensity profiles of molecular ions tracing different vertical layers of the disk allows to disentangle the contributions of the competing high energy ionization sources to the total H₂ ionization rate. Consequently such observations allow to constrain the SP flux in disks. Such a method is likely to be model dependent and ancillary observations constraining the vertical chemical structure of disks are required.

We have shown that stellar energetic particles can be an important ionization agent for disk chemistry. Modelling of spatially resolved observations of molecular ions with a model such as presented here allows to put first constraints on the stellar particle flux in disks around T Tauri stars.

Further model improvements concerning the stellar energetic particle transport (i.e. magnetic fields) are required to answer the question to what extent stellar particles reach the disk. Additionally non-stellar origins (i.e. jets) of high energy particles should be considered. With such models and spatially resolved molecular ion observations it will be possible to put stringent constraints on stellar energetic particle fluxes of T Tauri stars and to infer properties of the stellar and disk magnetic fields.

---

Acknowledgments

The authors thank the anonymous referee for useful suggestions and comments. We want to thank C. Dominik for very helpful discussions on stellar particle transport in disks. The research leading to these results has received funding from the European Union Seventh Framework Programme FP7-2011 under grant agreement No. 284405. R.CH. acknowledges funding by the Austrian Science Fund (FWF): project number P24790. M.P. acknowledges funding from the European Unions Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 664931. The computational results presented have been achieved using the Vienna Scientific Cluster (VSC). This publication was supported by the Austrian Science Fund (FWF). This research has made use of NASA’s Astrophysics Data System. All figures were made with the free Python module Matplotlib (Hunter 2007).

References

Appendix A: X-ray radiative transfer

Here we have used the newly developed X-ray radiative transfer module of PRODIMO. The details of this module will be presented in a forthcoming paper Rab et al. (in prep.). Here we only briefly describe the implementation and present results concerning the X-ray ionization rate for comparison with other models.

The main difference between the new model and the existing implementation of Aresu et al. (2011), Meijerink et al. (2012) is the treatment of X-ray scattering. We used the absorption and scattering cross-section from the publicly available xraylib library¹ (Schoonjans et al. 2011). For the X-ray radiative transfer we use the same numerical approach (e.g. use of wavelength bands) as is used for the dust radiative transfer in PRODIMO (Woitke et al. 2009). We assumed isotropic scattering but applied a correction factor (anisotropic factor g) to account for anisotropic Compton scattering. According to Cleeves et al. (2016), neglecting anisotropic scattering has only a limited impact on the X-ray ionization rate of about a factor of two. We consider this possible deviation as not significant for the results presented here. We further disregarded the dust in the X-ray radiative transfer, since with the assumed gas to dust ratio of 100 and dust settling X-ray photons mainly interact with the gas (Bethell & Bergin 2011; Glassgold et al. 2012).

The X-ray chemistry in PRODIMO is presented in detail in Meijerink et al. (2012). The interaction of X-rays with the gas changes the species abundances and consequently also the X-ray opacities. We therefore iterated between the X-ray radiative transfer and the chemistry until convergence is reached.

In Figs. A.1 and A.2 we show the X-ray ionization rate ζ_X for our disk model (Sect. 2.4) and for both X-ray spectra (Sect. 2.3.1). Figure A.1 shows ζ_X as a function of vertical column density at 1 au whereas Fig. A.2 shows ζ_X in the midplane of the disk as a function of radius. These results are in good quantitative agreement with Ercolano & Glassgold (2013) and Cleeves et al. (2015) considering that different disk models and different implementations of X-ray radiative transfer (e.g. cross-sections) are applied.

Appendix B: Chemistry

For the chemical reaction network we use the UMIST 2012 database for gas phase chemistry (McElroy et al. 2013) for a selection of 235 species. Additionally the network includes X-ray chemistry (Aresu et al. 2011; Meijerink et al. 2012), PAH (polycyclic aromatic hydrocarbons) chemistry (Thi et al. 2014), adsorption and desorption (thermal, CR and photo desorption) for ice species and H₂ formation on grains (Woitke et al. 2009). In total the network includes 3143 chemical reactions. For the ice species we use the binding energies from the UMIST 2012 release but updated a couple of values for oxygen bearing species (see Sect. B.2). Further details on the network can be found in Kamp et al. (2017). The most relevant gas phase reactions for HCO⁺ and N₂H⁺ with their rate coefficients are listed in Table B.1.

To test the robustness of our chemical model with respect to HCO⁺ and N₂H⁺ we performed chemical tests with time-dependent chemistry, varying binding energies and the KIDA chemical network (Kinetic Database for Astrochemistry, Wakelam et al. 2012, 2015). The different test models are described in Table B.2 and discussed in the following sections. The results for the HCO⁺ and N₂H⁺ column densities are shown in Figs. B.1 and B.2.

Appendix B.1: N₂ shielding

Dust shielding of N₂ photodissociation is strongly reduced in protoplanetary disks due to dust evolution (e.g. dust growth) and shielding by H₂ becomes important (Li et al. 2013). To account for this we implemented the H, H₂ shielding and self-shielding functions for N₂ of Li et al. (2013) in our chemistry model. The data are taken from the Leiden photodissociation database² .

We find that N₂ shielding is very important for the abundance structure of N₂H⁺. In models without N₂ shielding the vertical column density of N₂H⁺ is reduced by one order of magnitude throughout the disk. Further the warm N₂H⁺ layer above the vertical CO ice-line vanishes.

Appendix B.2: Binding energies

Appendix B.2.1: Oxygen binding energy

Several laboratory experiments (e.g. Ward et al. 2012; He et al. 2014; Minissale et al. 2016a,b) reported binding energies for oxygen in the range of E_B(O) = 1500 − 1800 K. This is significantly higher than the value of E_B(O) = 800 K listed in the UMIST 2012 database. We therefore updated the binding energies for oxygen and several other oxygen bearing species with the values listed in Minissale et al. (2016b). The new values with their references are given in Table B.3.

The higher binding energy for oxygen has a significant impact on N₂H⁺. Using the UMIST 2012 binding energy reduces the N₂H⁺ column density by about a factor of three for r ≳ 50 au (see model EBUMIST in Fig. B.1). Similar to CO, oxygen also destroys N₂H⁺ via the reaction O + N₂H⁺ → N₂ + OH⁺. Due to this reaction the N₂H⁺ abundance is strongly reduced in a thin layer below the vertical CO ice line and above the vertical oxygen ice line. With the higher binding energy for oxygen this thin layer vanishes as oxygen freezes-out at higher temperatures than CO.

For HCO⁺ the dip in the column density profile at r ≈ 50 au alone is significantly affected. In this region HCO⁺ is efficiently destroyed by water. Removing oxygen and OH from the gas phase reduces also the water abundance near the vertical CO ice line, consequently the HCO⁺ abundance increases.

Appendix B.2.2: N₂ binding energy

In the UMIST 2012 database E_B(N₂) = 0.67 × E_B(CO) (see Table B.3). This ratio is close to the value of 0.65 used to explain observations of prestellar cores showing that N₂ freezes out at lower temperatures than CO (e.g. Bergin & Langer 1997; Bergin et al. 2001; Ceccarelli & Dominik 2005). Also disk observations of N₂H⁺ indicate that E_B(N₂) <E_B(CO) (Qi et al. 2003, 2013b). However, the ratio of 0.65 is in contradiction with laboratory experiments which show E_B(N₂) ≈ 0.9 × E_B(CO) (Öberg et al. 2005; Bisschop et al. 2006; Fayolle et al. 2016). One explanation for this discrepancy might be differential freeze-out but it is unclear if this is sufficient (Bisschop et al. 2006; Fayolle et al. 2016).

In our models the N₂H⁺ column density drops by about one order of magnitude if we set E_B(N₂) = 0.9 × E_B(CO) (model EBN2 in Fig. B.1). Contrary to our model Aikawa et al. (2015) find that in their disk model the N₂H⁺ abundance structure is not strongly affected by setting E_B(N₂) = E_B(CO). They argue that due to the sink effect, CO is depleted anyway and therefore the N₂H⁺ abundance is not strongly affected by setting E_B(N₂) = E_B(CO). We have not included dust surface chemistry in our model therefore the CO sink effect is less efficient. Aikawa et al. (2015) argue that the sink effect depends on various parameters (turbulence in the disk, efficiency of the conversion of CO to less volatile species) therefore they also studied a no-sink model. In this no-sink model the N₂H⁺ abundance is sensitive to the binding energy of N₂, very similar to our models.

The main difference between the sink and no-sink model is the location of the CO ice line and the N₂H⁺ layer. In the sink model CO is also depleted at temperatures higher than its sublimation temperature, consequently the CO ice line moves to higher and warmer layers of the disk. In the no-sink model the CO ice line is given by the adsorption/desorption equilibrium for CO which depends on the chosen CO binding energy. However, with the exception of the exact location of the CO ice line and the N₂H⁺ layer our models with E_B(N₂) = 0.67 × E_B(CO) are in good agreement with the full chemical network/sink model of Aikawa et al. (2015). We discuss the impact of the CO ice line location on our results in Sect. 4.2.1.

Appendix B.3: Comparison to the KIDA chemical network

The reaction rates for nitrogen chemistry are not as well known as for carbon/oxygen chemistry (e.g. Hily-Blant et al. 2010; Le Gal et al. 2014). Wakelam et al. (2013) reviewed a large number of important reactions for nitrogen chemistry. The new derived reactions rates are in included in the latest KIDA gas phase chemistry database release (Wakelam et al. 2015).

We ran models using the KIDA gas phase chemistry database instead of the UMIST 2012 database (see Table B.2). We used the KIDA 2011 (Wakelam et al. 2012) and the KIDA 2014 (Wakelam et al. 2015) releases. The additional chemistry included in PRODIMO (e.g. X-ray chemistry) remains the same.

Figure B.1 shows the resulting HCO⁺ and N₂H⁺column densities for the KIDA2011 and KIDA2014 model in comparison to the reference model (UMIST 2012). In the KIDA2011 model the N₂H⁺ column density is about a factor of three higher compared to the reference and KIDA2014 models. The reason are updated rate coefficients for the dissociative recombination reactions of N₂H⁺ with electrons. The reaction rates at 20 K for N₂H⁺ + e⁻ → N₂ + H (N₂H⁺ + e⁻ → NH + N) are a about a factor of five (three) higher in the KIDA 2014/UMIST 2012 releases than in the KIDA 2011 release. This explains the higher N₂H⁺ abundance in the KIDA 2011 model. However, we find a good agreement for the KIDA 2014 and the UMIST 2012 release for both molecules HCO⁺ and N₂H⁺.

Appendix C: Steady-state versus time-dependent chemistry

For our models we have assumed that the chemistry reaches a steady-state within typical lifetimes of disks (couple of million years). To verify this assumption we run time-dependent chemistry models for the reference models CI_XN and CL_XN. We find that for the species considered in this work HCO⁺ and N₂H⁺ the chemistry reaches steady-state within ≈1 Myr in a large fraction of the disk. This is in agreement with the models of Aikawa et al. (2015) although they use a different, in particular larger chemical network.

In Fig. B.2 we compare the vertical column densities N_HCO⁺ and N_N2H⁺ of the steady-state model to the results of the time-dependent models at different times. For r ≳ 30 au and t ≳ 1 Myr to column densities of the steady-state and the time-dependent model are nearly identical. However, the steady-state model underestimates N_HCO⁺ and N_N2H⁺ for r ≲ 30 au. The differences are caused by the “sink effect” for CO and N₂ (see Sect. 3.3 and Aikawa et al. 2015). This erosion of CO and N₂ and other neutral molecules in the disk midplane by the reaction with He⁺ is a slow chemical process (e.g. Furuya & Aikawa 2014; Bergin et al. 2014; Helling et al. 2014) and is therefore over-estimated in the steady-state models.

As N₂H⁺ resides in deeper layers than HCO⁺ the deviations in the column densities are more pronounced. As a consequence the N₂H⁺ column density does not exactly trace the CO ice line in the midplane of the disk. The midplane CO ice line is located at r ≈ 12 au in the time-dependent model, whereas in the steady-state model N_N2H⁺ indicates a CO ice line at r ≈ 30 au. We note that in the low CR model CL_XN the sink effect is less efficient due to the lower midplane ionization rate, therefore the deviations in the steady-state models for N_HCO⁺ and N_N2H⁺ are smaller (N_N2H⁺ indicates a CO ice line at r ≈ 20 au).

Our comparison shows that for the bulk of N_HCO⁺ and N_N2H⁺ the assumption of steady-state chemistry is well justified. In particular our conclusions concerning the impact of SP ionization are not affected by the artifacts in the steady-state models. A steady-state model requires about a factor of ten less computational time than the time-dependent models. The use of the steady-state models allows us to study other important aspects like the impact of CO depletion (Sect. 4.2.2) or the comparison of different chemical networks.

Appendix D: The Turner model

In Turner & Drake (2009) the SP ionization rate is calculated for a mean solar mass nebulae disk model. To calculate the SP ionization rate at the surface of the disk they scaled their CR ionization rate of ζ_CR = 5 × 10^-18 s^-1 by a factor of 10⁴(r/ au)^-2. To account for the attenuation of SPs as a function of hydrogen column density they applied the same equation as they use for Galactic cosmic rays (their Eq. (2)).

To compare our results with the approach of Turner & Drake (2009) we implemented their method in our model. We also used their X-ray input spectrum and applied our X-ray radiative transfer to calculate the X-ray ionization rate. In Fig. D.1 we show the resulting X-ray, SP and CR ionization rates as a function of radius at a vertical column density of ≈8 g cm^-2 (≈3 × 10²⁴ cm^-2) for the Turner model (CI_T) and for our CI_XN_SP model. A comparison of Fig. D.1 with Fig. 1 of Turner & Drake (2009) shows that the CI_T model reproduces their SP ionization rates. We note that in our disk model (i.e lower disk mass) only for r< 10 au a vertical column density of >3 × 10²⁴ cm^-2 is reached; therefore only the inner 10 au are shown.

Figure D.1 clearly shows the differences of the two approaches. In our models SPs can only penetrate the disk radially and therefore cannot reach the disk midplane (see Sect. 3.2). In the Turner model the particles can also penetrate the disk vertically down to the midplane as they use the same equation for SP attenuation as is used for Galactic CRs. We cannot reproduce their results if we allow only for radial transport of particles. The differences in the X-ray ionization rate are due to the different X-ray spectrum (L_X = 2 × 10³⁰ erg s^-1,T_X = 5.8 × 10⁷ K, see Table 2) used by Turner & Drake (2009).

Appendix E: Comparison to observations

We compared the HCO⁺J = 3 − 2 and N₂H⁺J = 3 − 2 line fluxes from our reference model to the observations of Öberg et al. (2010, 2011a). From their sample we choose the targets where both lines were detected (the six targets shown in Fig. 6 of Öberg et al. 2011a). The observed integrated line fluxes are in the range of 3.1 − 17.8 Jy km s^-1 for HCO⁺J = 3 − 2 and in the range of 0.4 − 2.9 Jy km s^-1 for N₂H⁺J = 3 − 2 (all line fluxes are scaled to a distance of 140 pc). Excluding IM Lup, an extremely large and massive disk, from the sample gives a much narrower range of 3.1 − 5.4 Jy km s^-1 and 0.4 − 1.4 Jy km s^-1 for HCO⁺J = 3 − 2 and N₂H⁺J = 3 − 2, respectively.

To calculate the line fluxes for our model we use the line transfer module of PRODIMO (Woitke et al. 2011) and the molecular data from the Leiden Atomic and Molecular Database (Schöier et al. 2005; Botschwina et al. 1993; Flower 1999). For our reference model CI_XN we find fluxes of 2.4 Jy km s^-1 and 0.78 Jy km s^-1 for HCO⁺J = 3 − 2 and N₂H⁺J = 3 − 2, respectively. As the sample of Öberg et al. (2011a) is probably biased towards large and massive disks we also calculated the fluxes for a disk four times more massive. We find fluxes of 4.3 Jy km s^-1 and 1.18 Jy km s^-1 for HCO⁺J = 3 − 2 and N₂H⁺J = 3 − 2, respectively. These results are well within the range of the observations.

All Tables

All Figures