© ESO, 2014

1. Introduction

The asymptotic giant branch (AGB) is the last nuclear-burning phase for low- to intermediate-mass stars (from 0.8 M_⊙ to ~8 M_⊙), where the core of the star evolves into an inert and degenerate C-O core (Herwig 2005). AGB stars can be classified by the elemental C/O ratio, namely, C-rich AGB stars (C/O > 1), M-type AGB stars (C/O < 1) and S-type AGB stars (C/O ≈ 1). Our own Sun will become an M-type AGB star. Generally, the envelopes of C-rich stars contain CO plus C-bearing molecules, such as C₂, CN, HCN, and C₂H₂, while O-rich stars contain CO plus O-bearing molecules such as H₂O, TiO, and VO (Le Bertre 1997). AGB stars undergo considerable mass loss and eject dust and molecules into the surrounding regions, creating circumstellar envelopes (CSEs). These gas and dust envelopes eventually merge with the interstellar medium (ISM), enriching molecular clouds in which new stars may be born. Furthermore, CSEs are one of the richest sources for detecting new molecules. The study of AGB stars is of particular interest and importance for our understanding of molecular formation, destruction, and the recycling of material between star birth and star death. The schematic structure of the CSE of an AGB star is shown in Fig. 1.

IRC +10216 (CW Leonis), the brightest object in the sky at mid-infrared wavelengths outside the solar system, is the nearest C-rich AGB star and attracts intensive theoretical and observational study (e.g., Morris 1975; Bieging & Rieu 1988; Glassgold 1996; Mauron & Huggins 1999; Millar et al. 2000; Cernicharo et al. 2000; Woods et al. 2003; Cordiner & Millar 2009; De Beck et al. 2012; Agúndez et al. 2012; McElroy et al. 2013). More importantly, it is one of the richest molecular sources in the sky. To date, around 180 molecular species, not counting isotopologues, have been identified in the interstellar medium or CSEs¹ while more than 80 of them have been detected in IRC +10216. Among these, more than 30 nitrogen-bearing species have been identified, including HCN, CN, SiN, PN, HNC, MgCN, NaCN, NH₃, HC₇N, and C₅N⁻, amongst others (Wakelam et al. 2010). Even water, which was thought unlikely to form in C-rich AGB envelopes, has been detected here (Melnick et al. 2001; Decin et al. 2010). Some molecular detections in IRC +10216 were their first discoveries in astrophysical environments, for example, the first detection of the cyanide anion CN⁻ (Agúndez et al. 2010), and FeCN (Zack et al. 2011).

Most of the detected species in IRC +10216 are found in the cool, expanding outer CSE, where photodissociation processes dominate the destruction of molecules. The investigation of “key” reactions, which might affect the entire chemical network and therefore the comparison of simulations and observations, is an important topic. According to a systematic sensitivity study using a large chemical network, the reaction ${\mathrm{N}}_{\mathrm{2}}\mathrm{+}\mathrm{h}\mathit{\nu }\mathrm{\to }\mathrm{N}\mathrm{+}\mathrm{N}$(1)is one of the most significant reactions in the outer CSE, which directly affects the abundances of many N-bearing species, such as N, N₂, HC₂N, C₃N, C₃N⁻ and C₂N (Wakelam et al. 2010). Eq. (1) is the primary destruction route of N₂ in any region where UV photons are present. From an observational point of view, the direct detection of N₂ is very challenging because it has no permanent electric dipole and thus possesses no electric-dipole-allowed pure rotational spectrum. The only reported detection of molecular nitrogen is via its far-UV electronic transitions observed in the interstellar medium (Knauth et al. 2004). One can infer N₂ indirectly through the protonated ion, N₂H⁺ (Turner 1974; Herbst et al. 1977) or its deuterated form, N₂D⁺. However, neither N₂H⁺ nor N₂D⁺ have been identified in IRC +10216. From the simulation point of view, Eq. (1) has usually not been treated properly in models because these have not included N₂ self-shielding. The importance of self-shielding of molecules in CSEs was first noticed for CO some 30 years ago (Morris & Jura 1983). Even with an approximate treatment of CO self-shielding (the “one-band approximation”), much better agreement was obtained between simulations and observations.

We here employ the latest reported photodissociation rate and shielding functions for N₂ (Li et al. 2013; Heays et al. 2014) to investigate the effects in a chemical model of the CSE of IRC +10216. These are based on a concerted laboratory (e.g., Ajello et al. 1989; Helm et al. 1993; Sprengers et al. 2004; Stark et al. 2008; Lewis et al. 2008a; Heays et al. 2011) and theoretical (e.g., Spelsberg & Meyer 2001; Lewis et al. 2005; 2008b; Ndome et al. 2008) effort over the last two decades. An update was also made to the photodissociation of CO using the self-shielding functions from Visser et al. (2009) following a similarly large experimental effort over the past decades. While the absolute unshielded rates of both N₂ and CO are changed only at the level of ~30% with respect to the previous values (van Dishoeck 1988), it is important to realize that the uncertainty in the photorates is reduced by an order of magnitude. Also, we develop and employ a new fully spherically-symmetric (SS) model that computes the self-shielding of molecules in an isotropic interstellar radiation field, rather than the usual plane-parallel (PP) geometry.

The paper is organised as follows: the CSE model, the improvements in N₂ and CO photodissociation, as well as the theory and details for the evaluation of photodissociation rates for the SS model are described in Sect. 2. The results and discussion can be found in Sect. 3, followed by the concluding remarks in Sect. 4. The impact of these improvements on the molecular abundances in the CSE of IRC +10216 are discussed. Special attention is given to those species which have already been detected and those which may be detectable in the near future, e.g., using the Atacama Large Millimeter/submillimeter Array (ALMA). A full description of the SS model used for calculating the photodissociation rates of CO and N₂, together with the proposed numerical methods for implementing molecular shielding functions, are included in Appendices A and B.

2. Methods

2.1. CSE Model

The CSE model described in McElroy et al. (2013) was employed and extended in the present work. Specific details can be found in Millar et al. (2000) and Cordiner & Millar (2009). Improvements in the model are summarised in Sect. 2.2. The latest (fifth) release of the UMIST Database for Astrochemistry (UDfA), McElroy et al. (2013), hereafter RATE12, was adopted in all of the calculations. RATE12 contains 6173 gas-phase reactions involving 467 species. McElroy et al. (2013) tested this network in a dark cloud model and in a CSE model of IRC +10216 and compared the results with previous models and observations.

The main assumptions and key parameters adopted in our model of IRC +10216 are summarised in Table 1. Parent species that are injected at the inner radius of the envelope are listed in Table 2. The simulations start at the inner envelope with the following set of parameters: the radius, molecular hydrogen number density, visual extinction, and kinetic temperature of the gas are r₁ = 1.0 × 10¹⁵ cm, n(H₂) = 1.3 × 10⁷ cm^-3, A_V = 13.8 mag, and T = 575 K, respectively. The outer radius of the envelope is set to r_f = 7 × 10¹⁷ cm, where the density has decreased to n(H₂) = 26 cm^-3, the temperature to T = 10 K, and the visual extinction to A_V = 0.02 mag. We pay particular attention to the chemistry in the outer CSE where photodissociation is an important process. For all photoreactions, species are destroyed by photons from all directions, i.e., the radiation field is isotropic. For IRC +10216, one has to only take into account photons that come from the interstellar medium. The star, with a temperature of ~2330 K (De Beck et al. 2012), is too cool to generate photons that could lead to the photodissociation of molecules.

2.2. What’s new

Our aim is to improve the treatment of N₂ and CO photodissociation. In the work of McElroy et al. (2013), they employed the following parameters and methods:

• (a)

  Unshielded photodissociation rates

  CO:2.0 × 10^-10 s^-1, from van Dishoeck (1988)

  N₂: 2.0 × 10^-10 s^-1, from van Dishoeck (1988)

• (b)

  Shielding functions

  CO: dust + self-shielding, from Morris & Jura (1983), neither taking into account all lines for CO self-shielding nor shielding from H₂

  N₂: dust shielding only

• (d)

  Method for implementing self-shielding functions

  The “one-band approximation” (Morris & Jura 1983).

In this work, we employ the most accurate molecular data to date and develop a new method for implementing molecular shielding functions in the SS model:

• (a)

  Unshielded photodissociation rates

  CO:2.6 × 10^-10 s^-1, 30% higher, from Visser et al. (2009)

  N₂: 1.65 × 10^-10 s^-1, 28% lower, from Li et al. (2013)

• (b)

  Shielding functions

  CO: dust + self- + H + H₂ shielding, from Visser et al. (2009)

  N₂: dust + self- + H + H₂ shielding, from Li et al. (2013)

• (c)

  Method for implementing self-shielding functions

  By iterating the calculations to get converged results for the column densities in spherical symmetry for an isotropic radiation field. Details of this method are described in Sect. 2.3 and Appendix B.

2.3. Photodissociation rate and shielding function

The definitions of the photodissociation rate, k, unshielded rate, k⁰, and shielding function, Θ, appropriate for a plane-parallel model (only considering photons from the normal direction) are the same as those described in our recent paper (Li et al. 2013). In this work, we consider photons from all directions in space (SS model) for photoreactions that occur in the outer CSE. Suppose an incident ray, inclined with angle ϑ to the outward normal direction, arrives at a specific radius, r_i. The total column density of H₂ integrated from this point along the ray to infinity can be written as (Jura & Morris 1981; Morris & Jura 1983) ${\mathit{N}}_{{\mathrm{H}}_{\mathrm{2}}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\right)}\mathrm{=}{\mathit{N}}_{{\mathrm{H}}_{\mathrm{2}}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,}\mathrm{0}\mathrm{\right)}\mathit{\vartheta }\mathit{/}\sin \mathit{\vartheta }{\mathrm{cm}}^{-2}\mathit{.}$(2)This expression diverges for ϑ close to zero and π. In practice, the value of ϑ/ sinϑ equals 1 if ϑ = 0. For the case of ϑ = π, its value is obtained using a linear extrapolation from the preceding data points on our angular grid.

In the CSE of IRC +10216, most hydrogen is locked in H₂ due to self-shielding. Thus, the corresponding visual extinction, A_V(r_i,ϑ), is ${\mathit{A}}_{\mathit{V}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\right)}\mathrm{=}\mathrm{2}{\mathit{N}}_{{\mathrm{H}}_{\mathrm{2}}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\right)}\mathit{/}\mathrm{(}\mathrm{1.87}\mathrm{\times }{\mathrm{10}}^{\mathrm{21}}\mathrm{)}\mathrm{mag}\mathit{,}$(3)with the conversion factor based on the observations of Bohlin et al. (1978) and Rachford et al. (2009). The photodissociation rates of most interstellar species are mainly attenuated by dust, which is characterised by dust shielding functions given by ${\mathrm{\Theta }}_{\mathrm{dust}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\right)}\mathrm{=}{\mathrm{e}}^{\mathrm{-}\mathit{\gamma }{\mathit{A}}_{\mathit{V}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\right)}}\mathit{,}$(4)where γ is the parameter used to estimate the increased dust extinction at ultraviolet wavelengths. In general, the contribution of a specific ray of interstellar photons to the overall photodissociation rate at r_i may be calculated using $\mathit{k}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\right)}\mathrm{=}{\mathit{k}}^{\mathrm{0}}{\mathrm{\Theta }}_{\mathrm{dust}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\right)}{\mathrm{s}}^{-1}\mathit{.}$(5)However, the photodissociation rates of N₂, CO, and H₂ are more complex due to self- and mutual-shielding by H, H₂, CO, and other molecules (which are wavelength and column density dependent) and also continuum shielding by dust. A convenient way to take these effects into account in models is by use of molecular shielding functions that are calculated from high accuracy wavelength-dependent absorption cross sections. There is no need to consider H₂ photodissociation in the current study because it is fully shielded everywhere in the CSE. For N₂ and CO, one must consider self-shielding and mutual-shielding by H₂. The photodissociation rate of N₂ is calculated by ${\mathit{k}}_{{\mathrm{N}}_{\mathrm{2}}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\right)}\mathrm{=}{\mathit{k}}_{{\mathrm{N}}_{\mathrm{2}}}^{\mathrm{0}}{\mathrm{\Theta }}_{\mathrm{dust}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\right)}{\mathrm{\Theta }}_{\mathrm{mol}}\mathrm{\left[}{\mathit{N}}_{{\mathrm{H}}_{\mathrm{2}}}\mathrm{\right(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\left)}\mathit{,}{\mathit{N}}_{{\mathrm{N}}_{\mathrm{2}}}\mathrm{\right(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\left)}\mathrm{\right]}{\mathrm{s}}^{-1}\mathit{,}$(6)where ${\mathit{k}}_{{\mathrm{N}}_{\mathrm{2}}}^{\mathrm{0}}$ is the unshielded photodissociation rate of N₂, and N_N2(r_i,ϑ) is the column density of N₂ at radius r_i, integrated from this point along angle ϑ to infinity. Θ_mol [N_H2(r_i,ϑ),N_N2(r_i,ϑ)] is the molecular shielding function that is responsible for the multiple shielding effects caused by N₂ (self-) and H₂. It can be obtained via a bi-linear interpolation of N₂ and H₂ shielding functions that depend only on their column densities.

Correspondingly, the CO photodissociation rate is calculated using ${\mathit{k}}_{\mathrm{CO}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\right)}\mathrm{=}{\mathit{k}}_{\mathrm{CO}}^{\mathrm{0}}{\mathrm{\Theta }}_{\mathrm{dust}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\right)}{\mathrm{\Theta }}_{\mathrm{mol}}\mathrm{\left[}{\mathit{N}}_{{\mathrm{H}}_{\mathrm{2}}}\mathrm{\right(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\left)}\mathit{,}{\mathit{N}}_{\mathrm{CO}}\mathrm{\right(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\left)}\mathrm{\right]}{\mathrm{s}}^{-1}\mathit{,}$(7)where ${\mathit{k}}_{\mathrm{CO}}^{\mathrm{0}}$ is the CO unshielded photodissociation rate, and N_CO(r_i,ϑ) is the column density of CO at radius r_i, integrated along angle ϑ to infinity. Atomic hydrogen can also shield N₂ and CO, although in practice this is unimportant for the case of CSEs because the majority of hydrogen is locked in H₂.

A full description of the SS model and the numerical method for the computation of N(r_i,ϑ) can be found in Appendix A. The implementation of N₂ and CO shielding functions in the SS model is not straightforward. This is because before starting the evolution of the chemistry, one needs to know in advance the abundances of N₂ and CO at each radius to evaluate their column densities, and to compute the molecular shielding functions. However, these data are generated from the output of the model. One approximate solution is to assume that the N₂ and CO number densities are a constant proportion of H₂, which falls as r^-2 (Jura & Morris 1981), then work out the corresponding column densities. This assumption is inaccurate in the photon-dominated regions where N₂ and CO densities are no longer constant with respect to H₂. In this study, we propose a new procedure which gives more accurate and convergent results by iterating the computation of the chemistry. This method is described in Appendix B.

The photodissociation of N₂ and CO becomes considerable at the edge of the CSE, where the visual extinction is less than 1.0 mag and the temperature is around 10 K, as can be seen in Fig. 2. When calculating the N₂ shielding functions, the excitation temperature for both H₂ and N₂ are chosen to be 10 K. The Doppler widths of N₂, H₂, and H are 0.2, 3, and 5 km s^-1, respectively. The column density of atomic hydrogen is taken to be 1.0 × 10¹⁴ cm^-2. For CO shielding, the excitation temperature for H₂ and CO are chosen to be 11 and 5 K, respectively. At the low densities in the outer CSE, CO is subthermally excited. The Doppler widths of CO, H₂, and H are 0.3, 3, and 5 km s^-1. Electronic tables of the N₂ and CO shielding functions are available online². Note that these shielding functions were calculated by considering (H₂, H, and self-) shielding at various temperatures for a few typical column densities. Shielding functions for other temperatures and column densities can be obtained by interpolation. In practice, the uncertainties in the shielding functions over a temperature variation of a few degrees are negligible, see Li et al. (2013). The parameter γ used to account for the dust extinction at ultraviolet wavelengths for N₂ and CO was taken to be 3.9 and 3.5, respectively.

Finally, the radial photodissociation rate in a spherical geometry with isotropic incident radiation is evaluated by $\mathit{k}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}\mathit{k}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\right)}\sin \mathit{\vartheta }\hspace{0.17em}\mathrm{d}\mathit{\vartheta }{\mathrm{s}}^{-1}\mathit{.}$(8)

3. Results

3.1. Improvements in N₂ and CO photodissociation

Figure 3 compares the calculated column densities and photodissociation rates for the spherically-symmetric and plane-parallel models. For the PP case, we consider photons from the radial direction only but the full shielding (dust + self- + H + H₂) was included. The envelope was illuminated by a Draine (1978) field, extended to wavelengths longer than 2000 Å according to van Dishoeck & Black (1988) and with a scaling factor, χ. In the SS model, photons from all directions in space were considered with χ = 1 and the angular weighting scheme is given by Eq. (8). For the PP model, radiation from one side was taken into account explicitly with χ = 0.5 and radiation from the far side assumed completely attenuated. Hence, in the absence of the envelope the same unattenuated Draine (1978) radiation field is recovered as for the SS model. Indeed, for χ = 0.5, the N₂ photodissociation rates in the SS and PP cases are very similar (middle panel Fig. 3), whereas when assuming χ = 1.0, the PP rate is up to a factor of two higher at the outer edge. As a result, the PP model with χ = 1.0 underestimates the N₂ column density at this location, but the χ = 0.5 model results are very close to those of the SS model. The bottom panel of Fig. 3 shows that the PP and SS models start to deviate more in the inner envelope, where the PP model gives a much higher photodissociation rate. However, the absolute photodissociation rates are so small here that they are negligible in the chemistry. Overall, we conclude that for the CSE chemistry of IRC +10216, the PP model with a scaling factor of χ = 0.5 provides excellent agreement with the full SS model.

In the SS model the abundance of N₂ is higher than that in the PP model in the outer envelope. The distribution of CO in the SS and PP models is similar to that for N₂, also indicating a higher abundance of CO at the edge of the cloud for the SS case.

Including or excluding different shielding effects in the model is another crucial factor in the investigation of N₂ and CO photodissociation. Shielding effects, which act like a “smoke screen”, help molecules survive in photon dominated regions. N₂ and CO can be significantly shielded by both dust (wavelength independent) and molecules (wavelength dependent, line-by-line shielding). Fig. 4 presents the various shielding effects on N₂ photodissociation. Self-shielding is the most significant shielding effect, whereas (H₂ + H) shielding contributes the least. Shielding from dust largely depends on the properties of the dust (van Dishoeck et al. 2006). It is necessary to include all sources of shielding for N₂ because the combined effect has a strong influence on the abundance and distribution of N₂.

Figure 5 shows the fractional abundances of N₂ and N, calculated using various models. The location of the transition zone from N₂ to N is shifted outwards by a factor of 5 when molecular shielding is taken into account in the PP (χ = 1.0) model. Further changes happen when we use the SS model and PP (χ = 0.5) model.

The photodissociation of CO has an impact on both carbon and nitrogen chemistry. In particular, reactions with C⁺ are important destruction mechanisms for quite a few N-bearing species, for example C_nN⁻ (n = 1,3,5). As shown in Fig. 6, the conversion of CO to C and C⁺ is affected when CO self-shielding is calculated using the updated shielding functions. The chemistry of other C-containing molecules is also affected by the photodissociation of CO, but are beyond the scope of this work. The SS model with full shielding predicts that CO is relatively abundant up to a radius of ~ 6 × 10¹⁷ cm, which is in better agreement with observations than the PP model (which is described further in Sect. 3.3.2).

3.2. Impacts of N₂ and CO photodissociation on the CSE chemistry

The photodissociation of N₂ and CO significantly affects the chemistry in the outer CSE. This influence has two observable aspects: one leads to changes in the radial distribution of the species’ peak abundances. These distributions will become increasingly observable after exploiting the increased spatial resolution and sensitivity of new facilities such as ALMA. The other effect leads to changes in the peak abundances, which are directly reflected in total column densities. In the following subsections, we discuss fractional abundances of the most interesting (detectable or promisingly detectable) species which are sensitive to the photodissociation of N₂ and CO. Special attention is given to N-bearing species. Unless stated elsewhere, comparisons are made between the SS model including the updated N₂ and CO shielding and a recreation of the model of McElroy et al. (2013). Additional comparisons between observations and simulations from this study and McElroy et al. (2013) with regards to the 46 detected species in IRC +10216 are discussed in Sect. 3.3.

The dominant formation and destruction mechanisms of most species vary with physical conditions in the CSE. Here, we mainly focus on the chemistry of the outer CSE where photodissociation processes dominate.

3.2.1. CN⁻, C₃N⁻ and C₅N⁻

When N₂ photodissociation is treated in an appropriate manner, a higher abundance of N₂ is present in the outer CSE with a corresponding decrease in N. The most direct effects are on those N-bearing species whose formation processes require atomic N. C_nN⁻ (n = 1,3,5) have already been observed in the CSE of IRC +10216 (Agúndez et al. 2010; Thaddeus et al. 2008; Cernicharo et al. 2008). The following discussion of C_nN⁻ chemistry is based on the adopted reaction network from Walsh et al. (2009) which was incorporated into McElroy et al. (2013) and our model. Many of the adopted rate coefficients come from Eichelberger et al. (2007).

The chemistry of the molecular anions, CN⁻ and C₃N⁻, in the CSE of IRC +10216 has been described in Kumar et al. (2013). As shown in Fig 7, C_nN⁻ (n = 1,3,5) are very sensitive to N₂ photodissociation and are significantly decreased in the new treatment. This is because their formation processes are strongly affected by the abundance of atomic N in the outer CSE. As an example, at 1.0 × 10¹⁷ cm CN⁻ and C₃N⁻ are mainly formed through: $\mathrm{N}\mathrm{+}\mathrm{C}\begin{array}{c}\mathrm{-}\\ \mathit{n}\end{array}\mathrm{\to }\mathrm{CN}\mathrm{-}\mathrm{+}\mathrm{C}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{(}\mathit{n}\mathrm{=}\mathrm{5}\mathit{,}\mathrm{6}\mathit{,}\mathrm{7}\mathit{,}\mathrm{8}\mathit{,}\mathrm{9}\mathit{,}\mathrm{10}\mathrm{)}\mathit{,}$(9)$\mathrm{N}\mathrm{+}\mathrm{C}\begin{array}{c}\mathrm{-}\\ \mathit{n}\end{array}\mathrm{\to }\mathrm{C}\mathrm{3}\mathrm{N}\mathrm{-}\mathrm{+}\mathrm{C}\mathit{n}\mathrm{-}\mathrm{3}\mathrm{(}\mathit{n}\mathrm{=}\mathrm{6}\mathit{,}\mathrm{7}\mathit{,}\mathrm{8}\mathit{,}\mathrm{9}\mathit{,}\mathrm{10}\mathrm{)}\mathit{.}$(10)The abundances of these anions decrease with the decreasing abundance of N. As shown in Fig. 5, the fractional abundance of N drops from 4 × 10^-4 (McElroy et al. (2013) model) to 3 × 10^-5 (SS model) at 10¹⁷ cm, leading to ~1 order of magnitude drop in CN⁻ and C₃N⁻. Interestingly, C₅N⁻ has a different formation route. It comes primarily from C₅N, via $\mathrm{C}\mathrm{5}\mathrm{N}\mathrm{+}\mathrm{e}\mathrm{-}\mathrm{\to }\mathrm{C}\mathrm{5}\mathrm{N}\mathrm{-}\mathrm{+}\mathit{h\nu }\mathit{.}$(11)The radiative electron attachment of C₅N in Eq. (11) is very efficient compared with C₃N, which leads to the higher abundance of C₅N⁻ (see, e.g., Herbst & Osamura 2008).

3.2.2. CN, C₃N and C₅N

The chemistry of C_nN (n = 1, 3, 5) are also sensitive to N₂ photodissociation, but not as sensitive as C_nN⁻ (n = 1, 3, 5). This is because they are only partly affected by the available abundance of atomic N. Take CN as an example, it is mostly formed by the following reactions, $\begin{array}{ccc}& & \\ & & \\ & & \mathrm{N}\mathrm{+}\mathrm{CH}\mathrm{\to }\mathrm{CN}\mathrm{+}\mathrm{H}\mathit{.}\end{array}$The photodissociation of HCN dominates CN formation inside a radius of ~4 × 10¹⁶ cm. Further out, the reactions involving atomic N play a role, although these reactions do not affect its peak abundance, see Fig. 8. The cases of C₃N and C₅N are similar to CN.

There is a shift in the peak abundances of CN, C₃N⁻ and C₅N⁻ from the outside of the envelope inwards towards the star. In the case of C₅N⁻, this shift is from 10¹⁷ to 4 × 10¹⁶ cm. This difference is caused by the change in location of the N₂ to N transition zone when full shielding is included.

3.2.3. C₇N and C₉N

As shown in Fig. 9, the peak abundances of C₇N and C₉N are located at radii between 5 × 10¹⁶ and 1 × 10¹⁷ cm. Formation mechanisms for C₇N and C₉N are thought to be very similar to those of C_nN (n = 1, 3, 5). They have a formation route via: $\begin{array}{ccc}& & \mathrm{N}\mathrm{+}\mathrm{C}\begin{array}{c}\mathrm{-}\\ \mathrm{7}\end{array}\mathrm{\to }\mathrm{C}\mathrm{7}\mathrm{N}\mathrm{+}\mathrm{e}\mathrm{-}\mathit{,}\\ & & \mathrm{N}\mathrm{+}\mathrm{C}\begin{array}{c}\mathrm{-}\\ \mathrm{9}\end{array}\mathrm{\to }\mathrm{C}\mathrm{9}\mathrm{N}\mathrm{+}\mathrm{e}\mathrm{-}\mathit{.}\end{array}$These larger C_nN molecules have not yet been identified in IRC +10216, but they are also sensitive to N₂ photodissociation and may be useful in indirect investigations of the chemistry and abundance of N₂. Their predicted peak abundances are around one order of magnitude lower than that of C₅N.

3.2.4. HC_nN (n = 3, 5, 7, 9)

The cyanopolyynes, HC_nN (n = 3, 5, 7, 9), are also influenced by N₂ photodissociation, especially in the outer shell where photodissociation processes are important, as shown in Fig. 10. The fractional abundances of HC_nN (n = 3, 5, 7, 9) are all affected when the full shielding of N₂ is considered. The most significantly affected species is HC₅N. The changes in their total column densities are rather small because the peak abundances, which contribute most to the column densities, are unchanged. However, the radial distribution of HC₅N may be spatially resolvable by sensitive interferometers.

The chemistry of large organic molecules is more complex than for simple species. At the photon-dominated radius ~10¹⁷ cm, the main destruction process for the parent species, HCN, is photodissociation, which gives rise to its daughter molecule CN. Another important process for the destruction of HCN is by C⁺, via $\mathrm{HCN}\mathrm{+}\mathrm{C}\mathrm{+}\mathrm{\to }\mathrm{CNC}\mathrm{+}\mathrm{+}\mathrm{H}\mathit{.}$(17)However, there are quite a few reformation reactions at this location. Amongst the most important are two reactions relevant to the present study, $\begin{array}{ccc}& & \mathrm{N}\mathrm{+}\mathrm{CH}\mathrm{2}\mathrm{\to }\mathrm{HCN}\mathrm{+}\mathrm{H}\mathit{,}\\ & & \mathrm{H}\mathrm{+}\mathrm{CN}\mathrm{-}\mathrm{\to }\mathrm{HCN}\mathrm{+}\mathrm{e}\mathrm{-}\mathit{.}\end{array}$The decrease of HCN in the outer CSE, shown in Fig. 10, is largely due to the enhanced shielding of N₂, through the decreased abundance of N and, indirectly, CN⁻.

In our calculations, HC_nN are synthesised by the reactions $\mathrm{HC}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{H}\mathrm{+}\mathrm{CN}\mathrm{\to }\mathrm{HC}\mathit{n}\mathrm{N}\mathrm{+}\mathrm{H}\mathrm{(}\mathit{n}\mathrm{=}\mathrm{3}\mathit{,}\mathrm{7}\mathit{,}\mathrm{9}\mathrm{)}\mathit{.}$(20)The changes in the radial distributions for HC_nN (n = 3,7,9) are very similar to those of CN, see Fig. 8. HC₅N is an exception, and has the largest change among all HC_nN species in the photodissociation region. This is because the main formation process for HC₅N is via associative electron detachment, $\mathrm{H}\mathrm{+}\mathrm{C}\mathrm{5}\mathrm{N}\mathrm{-}\mathrm{\to }\mathrm{HC}\mathrm{5}\mathrm{N}\mathrm{+}\mathrm{e}\mathrm{-}\mathit{.}$(21)As discussed previously, the abundance and distribution of C₅N⁻ is strongly affected by the inclusion of N₂ shielding, thus, the abundance of HC₅N is also affected.

3.2.5. H₂CN and CH₂CN

The H₂CN and CH₂CN peak abundances are also significantly decreased when full N₂ shielding is included, as shown in Fig. 11. Their main formation processes in the outer CSE are $\mathrm{N}\mathrm{+}\mathrm{CH}\mathrm{3}\mathrm{\to }\mathrm{H}\mathrm{2}\mathrm{CN}\mathrm{+}\mathrm{H}$(22)and $\mathrm{N}\mathrm{+}\mathrm{C}\mathrm{2}\mathrm{H}\mathrm{3}\mathrm{\to }\mathrm{CH}\mathrm{2}\mathrm{CN}\mathrm{+}\mathrm{H}\mathit{.}$(23)There has already been a detection of CH₂CN in IRC +10216 (Agúndez et al. 2008), whereas H₂CN has only been detected in dark clouds (Ohishi et al. 1994). Since the column density computed here is close to the value observed in dark clouds, this species could also be detectable in IRC +10216. Because the chemistry of the two species are relatively simple and largely depend on the availability of N, we may be able to use them to constrain the fractional abundance of N₂ in space by comparing observations with models.

3.2.6. PN

The chemistry of P-bearing and Si-bearing species in C-rich CSEs are not completely understood (McElroy et al. 2013). The most sensitive P-bearing and Si-bearing species to N₂ shielding are HNSi, SiN, and PN, see Fig. 12. Here we concentrate only on the chemistry of PN, which is mainly formed through the following route at a radius of 10¹⁷ cm: $\mathrm{N}\mathrm{+}\mathrm{PH}\mathrm{\to }\mathrm{PN}\mathrm{+}\mathrm{H}\mathit{.}$(24)This process is largely reliant on atomic nitrogen, which is significantly decreased at the PN peak radius when updated shielding is used.

3.3. Comparison with observations

3.3.1. Column density

The radially-integrated column densities of most species are controlled by their peak abundances. After investigation, the abundances of quite a few of the 467 species in the network are affected when we use an appropriate method to calculated the photodissociation of N₂ and CO. A summary of the changes in column densities is an efficient way to extract the effects of these updates on the chemistry in the CSE. Here, we concentrate on those 46 species which have been observed in IRC +10216. The calculated and observed column densities are listed in Table 3. Of these, the calculated values for 14 species vary by more than 30%, and the values for 6 species vary by more than 80%, when compared with results from McElroy et al. (2013). The most affected N-bearing species with respect to the previous model are CN⁻ (−91%), SiN (−85%), and CH₂CN (−93%).

All of the species listed in Table 3 contain nitrogen and/or carbon and their formation and destruction may be influenced by N₂ and CO photodissociation. The important reactions for some of these species have been discussed in the previous subsections. The most significantly altered column densities include CN⁻ and C₃N⁻, whose formation is sensitive to the available N, which is reduced by the increased shielding of N₂ in our new calculations. A factor-of-10 decrease for SiN and CH₂CN column densities can be similarly explained.

Some of these changes reduce the agreement of the CSE chemistry model with observed column densities. There are many possible reasons for the discrepancies. First of all, some of the assumptions in the physical model of the CSE, unrelated to photodissociation (see Table 1), may be too simple. Further refinement of these assumptions may improve the agreement. For example, in this study, we ignored the presence of gas and dust shells in the expanding envelope. According to the study of Cordiner & Millar (2009), the inclusion of these shells gives a significant improvement in both modelled column densities and spatial distributions. A study of photodetachment as a destruction mechanism for the N-bearing anions, CN⁻ and C₃N⁻, in the CSE of IRC +10216, concluded that the inclusion of shells with enhanced density, similar to the model of Cordiner & Millar (2009), increases the column densities of the anions by about 20% (Kumar et al. 2013).

One of the challenges in CSE simulations is the uncertainty in the abundance of parent species (see Table 2). In particular the modelled column densities of Si-bearing and P-bearing species in Table 3 are several orders of magnitude lower than that suggested by observations. Possibly, the currently modelled Si-bearing and P-bearing parent species should have a higher abundance or there exist other Si- and P-containing parent species. On the other hand, the uncertainties of some observations are considerable. Take SiC₂ as an example. The cited column density for this molecule in Table 3 is 2.0 × 10¹⁴ cm^-2 (Thaddeus et al. 1984). However, according to a more recent study by Cernicharo et al. (2010) employing high quality Herschel/HIFI spectra an average column density of ~6.4 × 10¹⁵ cm^-2 was found. Our computed value, 2.6 × 10¹⁵ cm^-2, is in better agreement with the latter observed value.

3.3.2. Radius of peak abundance

The radial location of peak molecular abundance is another observable parameter and is of equal importance to the column densities. Many daughter species have a clear radial peak in their modelled abundances which we may directly compare with spatially resolved observations. A selection of such observations have been reduced to their peak radii (without considering detailed excitation or radiative transfer) and compared with our model calculations, and are summarised in Table 3 and Fig. 13.

The model predicts that most peak radii fall within shells between 1 × 10¹⁶ (4′′) and 2 × 10¹⁷ cm (89′′) from the central star. Here, the distance between the Earth and IRC +10216 is assumed to be 150 pc, based on the study by De Beck et al. (2012). Different values for this parameter have been employed in other work which may affect intercomparison, e.g., Huggins et al. (1988) who assumed 200 pc, and Cordiner & Millar (2009) who employed 130 pc.

The outer radius of CO using the SS model is slightly larger, ~7 × 10¹⁷ cm (~310′′), than the previous value, see Fig. 6. Emission from CO in IRC +10216 has been extensively investigated. A good summary can be found in De Beck et al. (2012), who also obtained various physical parameters of the structure of IRC +10216. However, conclusions from these observations are not always the same. For instance, Huggins et al. (1988) detected CO emission out to ~210′′ in the J = 1−0 line and 180′′ in the J = 2−1 line. Later, the J = 1−0 line was detected out to 190′′ (Fong et al. 2003, 2006). Recently, Cernicharo et al. (priv. comm.) observed weak J = 2−1 emission out to a radius of ~300′′, in excellent agreement with the present simulations.

Figure 13 compares the simulated and observed peak abundance radii for a few selected species for which interferometric maps are available. The observed peak radii, R_obs, were deduced from direct measurement of the peak radius and averaged over those maps for which multiple rotational transitions were available. These observations were done with high angular resolution telescopes, such as the Plateau de Bure Interfrometer (PdBI; Guélin et al. 1998, 1999; Guélin 2011) and the Very Large Array (VLA; Dinh-V-Trung & Lim 2008). In our models, it is assumed that the peak emission occurs at the peak abundance position, i.e., no excitation or radiative transfer of the model lines is performed.

With these assumptions, it is seen in Fig. 13 that the overall agreement between our simulations and observations is within a factor of two. Interestingly, the simulations predict a peak abundance at 18′′ for the parent species SiC₂ (see Table 3). This is very close to the value estimated from observations (Lucas et al. 1995; Guélin 2011), 16′′. No other parent species exhibits a similar maximum.

None of the peak radii of observationally-mapped species were found to be significantly altered by the introduction of a more appropriate photodissociation treatment, even if their column densities were substantially changed. However, large differences were deduced for some unmapped carbon chain molecules and anions, specifically C₃N⁻, C₅N⁻, C₅N, C₇N, and C₉N. The peak radii of these species all shift inwards after applying the new photodissociation shielding functions. For the case of C₅N⁻, the peak radius retreats by a factor of 2.5, from 1.0 × 10¹⁷ (44′′) to 4.5 × 10¹⁶ (20′′)  cm.

Our new model predicts that the CN⁻ peaks at a somewhat larger radius than its corresponding neutral CN, whereas the peaks of C₃N⁻ and C₅N⁻ match those of C₃N and C₅N, respectively, as shown in Table 3, Figs. 7, and 8. This conclusion differs from some previous models, e.g., Kumar et al. (2013) and references therein. Similar differences were found between the modelled peak radii of the neutral and anionic forms of C_nH (n = 2,4,6; see Table 3). Interferometric maps of these species could serve as a test of the model predictions.

3.4. N₂ abundance in IRC +10216

There is some uncertainty in the initial abundance of N₂ (a parent species) in the CSE of IRC +10216. In previous models, its adopted abundance relative to H₂ lay between 1 × 10^-5 and 2 × 10^-4 (e.g., Nejad & Millar 1987; Millar et al. 2000; MacKay & Charnley 2001; Agúndez & Cernicharo 2006). In the present work, we employ the same abundance as that in the model of McElroy et al. (2013), 2 × 10^-4, which gives an overall good agreement between the modelled and observed N-bearing species’ abundances in Table 3. However, the observational investigation of Milam et al. (2008) deduced an N₂ abundance of ~1 × 10^-7. They used the observation of PN to infer an N₂ abundance by assuming approximate equality of the abundance ratios, PN/N₂≃ [P]/[N] and [P]/[N] ≃ HCP/HCN ≃ 0.001.

To test the sensitivity of daughter species to the initial N₂ abundance we recomputed their chemical abundances after lowering the parent N₂ abundance to 1 × 10^-5, a factor of 20 lower than in our fiducial models. The results for four species, CN⁻, C₅N⁻, C₅N⁻, and N₂H⁺ are shown in Fig. 14.

N₂H⁺ has not been observed in IRC +10216, but has been detected in other interstellar and circumstellar environments (e.g., Liu et al. 2013; Qi et al. 2013; Codella et al. 2013). This species is thought to be a good tracer for N₂, from which it forms directly via the reaction, $\mathrm{N}\mathrm{2}\mathrm{+}\mathrm{H}\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}\mathrm{\to }\mathrm{N}\mathrm{2}{\mathrm{H}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\mathit{.}$(25)The utility of N₂H⁺ as an N₂ tracer is supported by the calculated abundances plotted in Fig. 14, which are directly proportional to the assumed N₂ abundance. The predicted column density of N₂H⁺ assuming initial N₂ abundances of 2 × 10^-4 and 1 × 10^-5 are 9.4 × 10¹⁰ and 5.3 × 10⁹ cm^-2, respectively. The peak abundance of N₂H⁺ is located at a radius of 5 × 10¹⁶ cm (22′′). Future deeper observations of IRC+10216 with ALMA may determine the column density of N₂H⁺ and allow for a new estimate of the N₂ abundance. Note that the abundances of neither N₂H⁺ nor HCO⁺, formed through a similar reaction between CO and H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$, are significantly affected by the new model treatment.

The abundances of C_nN⁻ anions are also not significantly affected by the altered parent N₂ abundance. These molecules are strongly dependent on the availability of atomic nitrogen (Sect. 3.2.1) which in the region of their peak abundance forms primarily from the photodissociation of HCN. Thus, HCN is the important N-bearing parent species affecting their column densities. Note that the main source of atomic N further out in the envelope is via the photodissociation of N₂.

4. Concluding remarks

In this work, the effect of self-shielding of N₂ and CO, and mutual-shielding by H₂, in the outer envelope of the C-rich AGB star, IRC +10216, was studied for the first time. This was performed using the latest available data for N₂ and CO photodissociation, and using an extended spherically-symmetric model for molecular self-shielding. The impact of these two improvements on the chemical evolution of the expanding envelope was investigated, with special attention paid to nitrogen chemistry and detectable species.

Key points from this study are:

(i): N₂ and CO are more abundant at the edge of the CSE than predicted by previous models due to molecular shielding from photodissociation.

(ii): The photodissociation of N₂ and CO affects the chemistry of most species in the outer CSE. Following its improved treatment here, the transition zones of N₂→ N and CO → C⁺→ C shift towards the outer edge of the envelope. The abundances of some species formed from N are reduced in the outer CSE. This induces large changes in the predicted column densities (factor of 10) of some species (e.g., C_nN and C_nN⁻ carbon chains) and the radii of their peak abundances.

(iii): Predictions are made for column densities and peak radii of molecules whose abundances are sensitive to N₂ and CO photodissociation and which could conceivably be detected. These may be verified by observations using high resolution and high sensitivity telescopes, especially ALMA and PdBI.

(iv): The abundances obtained from models which treat the full angularly-resolved spherically-symmetric radiation field incident on the outer CSE are very similar to those found in plane-parallel models in which the radiation is incident normally but with the intensity reduced by half. The proposed iterative method for implementing molecular shielding functions is most efficient when assuming that the molecule is fully shielded at the beginning of the calculation.

Other astrochemical models simulating the CSEs of AGB stars or other environments may show a similar sensitivity using a more realistic treatment of self-shielding and spherically-symmetric radiation as was found in this model of IRC +10216.

The photodissociation processes for N₂ and CO (both the unshielded rate as well as the depth dependence) are now very well understood, so any major discrepancies between observations and models are likely due to other assumptions in the model.

---

Acknowledgments

The authors would like to thank the anonymous referee for his/her valuable comments and constructive suggestions. X. Li is delighted to thank TJM for his hospitality during a two-week visit to Queen’s University Belfast (QUB) in 2013, and thank Dr. Markus Schmalzl and Prof. Xander Tielens for some stimulating discussions. Astrochemistry in Leiden is supported by the Netherlands Research School for Astronomy (NOVA), by a Spinoza grant and grant 648.000.002 from the Netherlands Organisation for Scientific Research (NWO), and by the European Community’s Seventh Framework Programme FP7/2007-2013 under grant agreements 291141 (CHEMPLAN) and 238258 (LASSIE). Astrophysics at QUB is supported by a grant from the STFC. C. W. acknowledges support from NWO (program number 639.041.335).

References

Appendix A: Spherically-symmetric model

In this appendix we describe the SS model employed in this work for the computation of the angle-dependent column densities of species at each radius. The SS model is shown in Fig. A.1. The star is located at point S and ejects dust and molecules to various shells with radii r_i (not equidistant). Since the star itself is too cool to generate photons that can induce photodissociation of molecules (in our case, IRC +10216, ~2330 K), one only needs to consider photons from the interstellar radiation field. For a specific position P at radius r_i, suppose ϑ is the angle between the direction of the incoming ray and the radius vector from the central star. Generally, the column density, N, of a species is defined as $\mathrm{\equiv }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathit{n}\hspace{0.17em}\mathrm{d}\mathit{z}{\mathrm{cm}}^{-2}\mathit{,}$(A.1)where n is the number density in units of cm^-3 and z is the integral path length in units of cm. In the present study, the total column density, N(r_i,ϑ), of a species integrated along a ray from point P to infinity is given by two cases:

• (1)

  0≤ϑ≤π/ 2, for example, $\overline{\mathrm{PA}}$, yielding $\begin{array}{ccc}& & \mathit{N}\overline{\mathrm{PA}}\mathrm{\left(}\mathit{r}\mathit{i}\mathit{,\vartheta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\left(\sqrt{{\mathit{r}}_{\mathit{i}\mathrm{+}\mathrm{1}}^{\mathrm{2}}\mathrm{-}{\mathit{r}}_{\mathit{i}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\vartheta }}\mathrm{-}{\mathit{r}}_{\mathit{i}}\cos \mathit{\vartheta }\right)\left[\mathit{n}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathrm{+}\mathit{n}\mathrm{(}\mathit{i}\mathrm{+}\mathrm{1}\mathrm{)}\right]\\ & & \mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\sum _{\mathit{j}\mathrm{=}\mathit{i}\mathrm{+}\mathrm{1}}^{\mathit{n}\mathrm{-}\mathrm{1}}\left(\sqrt{{\mathit{r}}_{\mathit{j}\mathrm{+}\mathrm{1}}^{\mathrm{2}}\mathrm{-}{\mathit{r}}_{\mathit{i}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\vartheta }}\mathrm{-}\sqrt{{\mathit{r}}_{\mathit{j}}^{\mathrm{2}}\mathrm{-}{\mathit{r}}_{\mathit{i}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\vartheta }}\right)\\ & & \mathrm{\left[}\mathit{n}\mathrm{\right(}\mathit{j}\mathrm{)}\mathrm{+}\mathit{n}\mathrm{(}\mathit{j}\mathrm{+}\mathrm{1}\mathrm{\left)}\mathrm{\right]}{\mathrm{cm}}^{-2}\mathit{,}\end{array}$(A.2)where n(i) is the number density of the molecule at radius r_i.

• (2)

  π/ 2<ϑπ, for instance $\overline{\mathrm{BP}}$. Since we assume the model is spherically symmetric, one may simplify the calculations by finding a point C that satisfies $\overline{\mathrm{SC}}$⊥$\overline{\mathrm{BD}}$, then get the column density along $\overline{\mathrm{BP}}$ with $\mathit{N}\overline{\mathrm{PB}}\mathrm{(}\mathit{r}\mathit{i}\mathit{,}{\mathit{\vartheta }}^{\mathrm{)}}\mathrm{=}\mathrm{2}\mathit{N}\overline{\mathrm{CB}}\mathrm{(}\mathit{r}\mathit{i}\mathit{,\pi }\mathit{/}{\mathrm{2}}^{\mathrm{)}}\mathrm{-}\mathit{N}\overline{\mathrm{PD}}\mathrm{(}\mathit{r}\mathit{i}\mathit{,\pi }\mathrm{-}{\mathit{\vartheta }}^{\mathrm{)}}{\mathrm{cm}}^{-2}\mathit{.}$(A.3)Both $\mathit{N}\overline{\mathrm{CB}}{}^{\mathrm{(}}\mathit{r}{}_{\mathit{i}}\mathit{,\pi }\mathit{/}{\mathrm{2}}^{\mathrm{)}}$ and $\mathit{N}\overline{\mathrm{PD}}{}^{\mathrm{(}}\mathit{r}{}_{\mathit{i}}\mathit{,\pi }\mathrm{-}{\mathit{\vartheta }}^{\mathrm{)}}$ can be conveniently calculated by case (1). Then, it is necessary to interpolate an abundance between the radially-gridded values in our model for points near the star, such as C.

Figure A.2 shows N(r_i,ϑ) for molecular N₂. According to our calculations, the major contribution of photons that induce photodissociation of N₂ come from the area within the dashed line (N(r_i,ϑ) < 1.0 × 10^-16 cm^-2), where r_i> 1.6 × 10¹⁷ cm and ϑ< 3π/ 4. Deeper in the envelope, the fractional abundance of N₂ stays constant because it is fully shielded.

Appendix B: Implementing the molecular shielding functions

In this Appendix we introduce the new procedure employed in the current work for the implementation of molecular shielding functions. Here we take N₂ as an example for the demonstration. The method has five steps.

• Step 1

  Set up the initial condition. Suppose N₂ is fully shielded everywhere in the CSE, then we can set k^initial(r_i) ≃ 0 for each radius, r_i, which means the corresponding shielding functions from dust and molecules are zero, i.e., ${\mathrm{\Theta }}_{\mathrm{dust}}^{\mathrm{initial}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\right)}\mathrm{=}\mathrm{0}$, and ${\mathrm{\Theta }}_{\mathrm{mol}}^{\mathrm{initial}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\right)}\mathrm{=}\mathrm{0}$. In all cases, we consider 201 incident rays as a function of ϑ between 0 and π, where ϑ = 0 gives the normal direction.

• Step 2

  Loop 1: compute the CSE chemistry for the first time, employing k^initial(r_i) prepared from Step 1. The calculated N₂ abundance generates column densities, N(r_i,ϑ), and corresponding dust and molecular shielding functions, ${\mathrm{\Theta }}_{\mathrm{dust}}^{\mathrm{loop}\mathrm{1}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\right)}$, and ${\mathrm{\Theta }}_{\mathrm{mol}}^{\mathrm{loop}\mathrm{1}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,\vartheta }\mathrm{\right)}$, and new photodissociation rates k^loop1(r_i).

• Step 3

  Loop 2: repeat the simulations of the CSE chemistry, but employing k^loop1(r_i) produced in Step 2 for each radius. Using the same method, generate another set of photodissociation rates, k^loop2(r_i).

• Step 4

  Compare the fractional abundances of N₂ at each radius for Loop 1 and Loop 2, to see if they are the same. If no changes are found, go to Step 5, otherwise go to Step 3, continuing the calculations.

• Step 5

  Stop the calculation. The outputs at this step are the final results.

The simulations can also approach convergence by assuming N₂ is unshielded in Step 1, namely k^initial(r_i) = k⁰ at each radius. In this case, the fractional abundance of N₂ will steeply decrease after starting the evolution of the chemical network, then approach the converged results quickly, as illustrated in Fig. B.1. The same converged results can be obtained from the two possible assumptions in Step 1. In practice, the fractional abundances of both N₂ and CO (as well as all other species) reach their converged abundances after ~15 loops. In practice, starting the calculation assuming full shielding in Loop 1 is most efficient.

The iterative method proposed here may be used in other models, where output variables are also required as inputs.

All Tables

All Figures