© ESO, 2012

1. Introduction

Ion-neutral reactions are the most important driving processes for gas-phase chemistry. Therefore it is important to understand the mechanisms by which chemical species in the interstellar medium become ionized, in order to have a more accurate picture of the chemistry in various interstellar sources. Near the edge of dense clouds and throughout diffuse clouds, UV photons can provide a powerful ionizing force upon the medium, especially if there is a sufficiently strong source of radiation nearby. These photons do not penetrate very far into dense clouds, decreasing exponentially with the column density. Other ionizing agents, like X-rays, will penetrate farther into dense clouds, but deep within the object, high-energy cosmic rays are the dominant ionizing force.

The recent detection of unexpectedly large abundances of H, however, in an assortment of diffuse clouds has raised the old question as to whether the high ionization rate needed is caused by a high flux of cosmic rays of  <1 GeV (McCall et al. 2003; Indriolo et al. 2007). Such low energy cosmic rays would not be expected to penetrate deeply into dense clouds, so that a column-dependent ionization rate due to cosmic rays might exist in denser sources. This question is best explored by examining the influence of cosmic ray ionization at different depths into a single object, and the Horsehead nebula is an ideal candidate for such an investigation.

The Horsehead nebula, also called Barnard 33, is a dark nebula of size about 5′ in the bright nebula IC 434. It is illuminated by σ   Ori from a distance of about 30′ (Anthony-Twarog 1982). The radiation field incident on the cloud is most commonly taken to be χ = 60 in Draine units (Draine 1978; Habart et al. 2005), and the geometry of the cloud is described as nearly “edge-on”, meaning that the line between σ   Ori and the Horsehead nebula is nearly perpendicular to the line of sight. This makes the Horsehead nebula ideal for observing column-dependent variables in a single source. It has an ambient magnetic field of  <6  μG (Zaritsky et al. 1986) and a steep density gradient ranging from 10² to 10⁵ cm^-3, and contains a pre-stellar core as well as at least one other dense region near the “throat” that will be able to be studied in greater detail by the Atacama Large Millimeter Array (ALMA) (Ward-Thompson et al. 2006).

High abundances of small carbon-bearing molecules were observed by Teyssier et al. (2004) and by Pety et al. (2005), with higher abundances of certain molecules (CCH, c-C₃H₂, C₄H) observed near the edge than at the center. This led Pety et al. to posit that polycyclic aromatic hydrocarbons (PAH’s) near the edge of the cloud are being destroyed by incident radiation, and that the products of their destruction are these small hydrocarbons. A number of other molecules have been detected, including the ions HCO⁺ and HOC⁺ (Goicoechea et al. 2009b), the carbon-bearing neutrals HCO (Gerin et al. 2009), l-C₃H, and c-C₃H (Teyssier et al. 2004), and the sulfur-bearing species CS and HCS⁺ (Goicoechea et al. 2006), although, except for HCO⁺, little information of their column dependence is available.

Chemical modeling of the Horsehead nebula was discussed by Winnewisser & Herbst (1993). Teyssier et al. (2004) provided the first detailed chemical PDR model of the Horsehead nebula, using the Meudon PDR code (Le Bourlot et al. 1993; Le Petit et al. 2002). A year later, Pety et al. (2005) modeled the Horsehead nebula with the same code, comparing the results with observations at three different lines of sight, and incorporating PAH’s into the model. Habart et al. (2005) determined a column-dependent temperature via thermal balance. None of these models is able to account for the high abundances of small hydrocarbons at the edge, or the HC₃N abundance.

Deuterium fractionation of HCO⁺ ( [DCO⁺] / [HCO⁺]  ~ 0.02) has been observed in the Cloud region at A_V ≈ 10, and used to constrain its temperature to about 20 K (Pety et al. 2007). Pety et al. (2007) also used this code to better understand deuterium fractionation at the Horsehead edge. Neutral atomic oxygen has also been detected (Goicoechea et al. 2009a), with hopes for Herschel’s heightened resolution to provide abundances of atomic oxygen for different regions in the cloud.

The effect of a higher sulfur abundance was considered by Goicoechea et al. (2006), using the Meudon code (Le Petit et al. 2006). The abundance of the negative ion C₆H ⁻  was calculated by Millar et al. (2007), although negative ions have not yet been observed in this region. Morata & Herbst (2008) developed a time-dependent PDR code, and first applied it to the Horsehead nebula, with mixed success. This is the code we make use of in this paper, in tandem with the Meudon PDR code, which we use to determine some of the physical conditions.

Compiègne et al. (2007) and Goicoechea et al. (2009b) have performed some recent modeling of the Horsehead region; Compiègne et al. (2007) explored the dust emission. Goicoechea et al. (2009b) self-consistently modeled the observed spatial distribution and line intensities with detailed depth-dependent predictions coupled with a nonlocal radiative transfer calculation for H¹³CO⁺, DCO⁺ and HOC⁺. They compared their model results with the Gerin et al. (2009) observations of HCO⁺ in order to constrain the electron fraction. Goicoechea et al. determined a very steep relative electron abundance of n_e/n_H ~ 10^-4 − 10^-8 (where n_H = n(H) + 2n(H₂)) at A_V ≈ 0.6 − 2.0 from the cloud edge, based on a faint emission line attributed to HCO⁺ near the edge.

In this paper, we report our investigation on the effect of a column-dependent cosmic ray ionization rate ζ(N_H) on model results for molecular abundances and their spatial variation in the Horsehead nebula. This is offered as a partial explanation of the high abundances of small hydrocarbons at the edge of the Horsehead nebula. In Sect. 2, we discuss the determination of three different ζ(N_H) functions, including a discussion of the role played by the magnetic field. In Sect. 3, we provide a detailed description of the PDR model used, and compare our calculated abundances with observational values using two different sets of elemental abundances. We also provide predicted abundances for observable species. In Sect. 4, we discuss the implications of these results, and a better determination of ζ(N_H) from single sources after the advent of ALMA.

2. The determination of ζ(N_H)

The cosmic ray ionization of the interstellar medium is caused primarily by relativistic protons, alpha particles, and electrons. This ionization rate, labeled ζ, is typically represented as a per-second rate at which cosmic rays ionize atomic hydrogen. Given the process where CR represents ionizing cosmic rays, ζ is defined by the kinetic equation where the brackets signify concentration. The ionization rate of other species, such as H₂ and He, is usually determined in chemical networks by multiplying ζ by a constant. Even near the edge of dense clouds, the majority of hydrogen is molecular in nature, so it is important to note that, to a good approximation, ζ_H2 ≈ 2ζ (Glassgold & Langer 1974).

In the last decade, results from diffuse sources (McCall et al. 2003; Le Petit et al. 2004; Indriolo et al. 2007), including recent observations with Herschel (Gerin et al. 2010; Neufeld et al. 2010), have most often indicated that in these environments ζ is more than an order of magnitude higher than the generally accepted value of 10^-17 s^-1. Earlier values for ζ ranging from 10^-17 − 10^-15 s^-1 had been proposed (Spitzer & Tomasko 1968; Hartquist et al. 1979; Dalgarno 2006). Table 1 contains a limited historical overview of some of the values of ζ utilized in previous models. These actually refer to molecular rather than atomic hydrogen.

The observations indicating a high ζ, along with this wide range of values, led us to initiate a calculation of column-dependent functions of ζ. At the same time, Padovani et al. (2009) undertook similar calculations. They used the ionization and energy loss cross sections for collisions between cosmic rays and atomic and molecular hydrogen (Cravens & Dalgarno 1978) as well as Helium to follow the flux-spectra of cosmic rays through a cloud and, from the flux spectra as a function of column, obtained the column-dependent cosmic-ray ionization rate for a number of initial flux-spectra. Here we report similar calculations but with a Monte Carlo approach in which we also include magnetic field effects.

2.1. Initial spectrum

We begin by considering the form of the initial cosmic ray flux-spectrum j(E) (cm^-2 s^-1 sr^-1 GeV^-1 per nucleon), as a function of energy. The spectrum has only been directly observed within our solar system, where the solar wind would have depleted the low energy cosmic rays (Parker 1958).

Different cosmic ray spectra have been proposed based on assumptions about the origin of the cosmic rays. Supernova shocks are currently the favored explanation for the origin of cosmic rays (Biermann et al. 2010; Axford 1981). The spectrum due to the supernova blast alone imposes a low-energy cutoff at about 100MeV because of energy loss due to debris and strong magnetic field effects (Hayakawa et al. 1961; Ip & Axford 1985). It is suspected that shocks in the debris may re-accelerate some of the thermalized cosmic rays (Ip & Axford 1985; Indriolo et al. 2009).

Shock models favor a steep power law for low-energy cosmic rays, with a new cutoff at 1 MeV, below which most cosmic rays would again lose a significant fraction of their energy into the debris, and would either be reabsorbed into the remnant, or would travel too slowly to propagate throughout the galaxy. Alternate theories for cosmic ray acceleration exist, but these also predict similar spectra for low-energy cosmic rays (Butt 2009).

Fig. 1 Three different cosmic ray flux spectra, taken from Hayakawa et al. (1961) (dashed line), Spitzer & Tomasko (1968) (dotted line), and Nath & Biermann (1994) (solid line).

Comparison between measurements of the cosmic ray flux and theoretical cosmic ray spectra have been very useful. Basic statistics, “leaky-box” models, convection methods, and Monte Carlo methods have been applied to better constrain cosmic ray spectra, often by examining the elemental composition of the cosmic rays themselves. Strong et al. (2007) contains an excellent review of these methods. Webber (1998) incorporated the newest results from Voyager into their Monte Carlo model, in order to determine the low energy spectrum better.

Nevertheless, Voyager is still in a region where solar winds have a substantial effect. In fact, the farther the Voyager satellite travels, the steeper the low energy spectrum becomes (Webber 1998). Indeed, as recently as Putze et al. (2011), statistical, Monte Carlo, and “leaky box” models have been unable to constrain the low energy cosmic ray spectrum, due to a lack of direct measurement of low energy cosmic ray protons outside the solar influence.

Indriolo et al. (2009) list many of the proposed cosmic ray spectra. We consider three representative spectra (Hayakawa et al. 1961; Spitzer & Tomasko 1968; Nath & Biermann 1994), which are shown in Fig. 1. These three spectra span the range of low energy cosmic rays. The spectrum from Spitzer & Tomasko (1968) is based on solar system measurements of the low energy cosmic ray flux, and contains the minimum low-energy cosmic ray spectrum. Nath & Biermann (1994) assume that the power-law for the cosmic ray spectrum at 1 GeV continues down until a hard cut-off at 1 MeV. Theirs is the highest published estimate of the low energy cosmic ray flux. We chose to use these three spectra in order to provide the full range of impact that different low energy cosmic ray flux spectra have on the ionization rate. The spectrum of Nath & Biermann (1994) increases the most steeply towards lower energies, that of Spitzer & Tomasko (1968) actually decreases towards lower energies, and that of Hayakawa et al. (1961) lies in the middle.

2.2. Cross sections

We calculate the loss of energy by considering 10⁴ cosmic ray protons with energies, E (in eV unless otherwise noted), distributed according to the three spectra selected above. The particles stream into a cloud of a number density n (cm^-3). At each distance increment, the particles are each assigned a random number, which is compared with the probability of an ionizing or other inelastic collision over an incremental distance, determined by cross-sections, σ (cm²). For ionizing collisions by protons we use the form of σ_i from Spitzer & Tomasko (1968). The cross section (cm²) for ionization of a hydrogen atom as a function of E and the rest-energy of the proton (E_P) is given by (1)The first term is dominant for “low” energies (500 keV  < E < 50 MeV), so for E < 50 MeV, σ_i,H ∝ 1/E, down to E ≈ 500 keV, when Eq. (1) ceases to be accurate. This cross-section is also used below for determining the cosmic ray ionization rate of atomic hydrogen. For molecular hydrogen, we simply multiply the cross section by a factor of 2.

Inelastic collisions are considered for atomic and molecular hydrogen only, and we use the cross-sections from Cravens et al. (1975), accounting for rotational, vibrational and electronic excitation as well as dissociation reactions. For the ionization of helium, the differential cross-section from Dalgarno et al. (1999) is integrated to yield a total cross section: (2)where A(E) ( ∝ 1/E for E ≲ 100 MeV;  ∝ log (E) for E ≳ 1 GeV) and ϵ₀ are parameters fit to the measurements of Shah et al. (1987).

2.3. Energy loss

The energy loss calculation assumes a column great enough that the cosmic rays will collide with gaseous atoms and molecules many times. Since our model is one-dimensional, we do not consider the effects of elastic collisions on the exclusion of low energy cosmic rays from molecular clouds.

Because there are many collisions, we are justified in utilizing the average energy lost by a cosmic ray in an ionizing collision, . This is equal to the ionization energy plus the average energy of the ejected electron. For molecular hydrogen, this is determined from the differential cross-section by Cravens & Dalgarno (1978) and Dalgarno et al. (1999) to be (3)where E (eV) is the energy of the cosmic ray before the ionizing event. The energy losses from other types of inelastic collisions with molecular hydrogen, as well as ionizing and other inelastic collisions with H and He are taken from the detailed forms in Cravens et al. (1975).

This energy loss is subtracted from the initial energy of the cosmic ray, and becomes the new energy. At each increment, a new flux-spectrum, j(E,N_H), is calculated, and new random numbers are assigned to the cosmic rays. Because of the energy-dependence of the σ functions, lower energy cosmic rays have more ionizing collisions. In the case of the spectrum of Nath & Biermann (1994), cosmic rays with E < 50 MeV contribute 99% to the value of ζ (see Sect. 2.5). To complicate matters, however, there is energy loss from magnetic effects in addition to the loss from collisions. Magnetic energy loss is assigned based on interactions with Alfvén waves, as discussed below, using a static magnetic field of 3 μGauss.

2.4. Magnetic field effects

Magnetic fields play an important role in the transport of cosmic rays. The Lorentz force is the largest magnetic force acting on cosmic rays, and affects energy loss by increasing the path length cosmic rays travel, as they spiral along the magnetic field lines. This resulting increase in path length is not, however, the primary source of energy loss. Rather, the dominant magnetic field effect on cosmic rays is due to irregularities in the magnetic field.

Because of the neutralization of low-energy cosmic rays, there will be far fewer cosmic rays at the center of the cloud than at the edge. Since cosmic rays are overwhelmingly positively charged, these losses introduce a charge imbalance in the cloud. Electrons are attracted to the edge, and their motion generates magnetic field irregularities moving from the center to the edge of the cloud with velocity υ_A = B(4πρ) ^− 1/2. These irregularities, called Alfvén waves, are the dominant source of energy loss, as discussed in Skilling & Strong (1976). Hartquist et al. (1978, 1979) first applied the work of Skilling & Strong (1976) to calculate cosmic ray ionization rates, and proposed different values of ζ, depending on the object.

Following Skilling & Strong (1976), we determine the charge imbalance using the Monte Carlo simulation with B = 0, and consider it in terms of a characteristic column density, λ(E) (cm^-2), determined by the simulation, at which the number of cosmic rays will be depleted by a factor of e. This means that N_H must be  ≳ λ(E) for cosmic rays of energy E to be significantly affected by magnetic field irregularities. This function will appear later in the analysis.

Alfvén waves are driven by the charge imbalance, and are damped by the friction between ions and the surrounding gas, as discussed by McIvor (1975) in terms of the collision rate Γ (s^-1) between ions and neutrals (Dalgarno & Dickinson 1968). The larger Γ is, the less effect the waves have. The static magnetic field enhances the damping by absorbing smaller irregularities. However, the larger static magnetic field also increases υ_A, and thus the frequency of collisions between the cosmic rays and the Alfvén waves.

This mechanism for cosmic ray energy loss by Alfvén wave effects is important for N_H < 10²⁴ cm^-2 when B ≲ 6 mG and n_H ≲ 10⁹ cm^-3. At a given column density, cosmic ray energy is substantially affected by Alfvén waves for energies less than the energy E₀. The static magnetic field outside denser regions is assumed to be much smaller than the field inside these regions. Because the difference between the magnetic field inside and outside the cloud significantly dampens the Alfvén waves for mid to high energy cosmic rays, E₀ cannot be greater than 50 MeV (Cesarsky & Volk 1978). E₀ is dependent on various physical parameters of the source in question. For typical cold and dense interstellar conditions, n(HI) = 1 cm^-3, n_H = 10⁴ cm^-3, T = 20 K, and B = 3 μG. Under these conditions, the use of j₀(E) from Nath & Biermann (1994) leads to E₀ = 1 MeV at N_H = 10¹⁹ cm^-2, while for N_H > 10²¹ cm^-2, E₀ = 50 MeV.

Integrating over energies up to this cutoff value, we can obtain the magnetohydrodynamic solution for j_IC(E,N_H), the “In-Cloud” cosmic ray flux-spectrum at a given N_H, to be (Skilling & Strong 1976) In this expression, j(E,N_H) is the spectrum determined using the Monte Carlo simulation in the absence of magnetic field effects, the magnetic energy density U_M = B²/2μ₀ (erg/cm³), Ω₀ is the gyromagnetic frequency (s^-1), the Compton-Getting factor α (Gleeson & Axford 1968) is where j₀ is the initial cosmic ray flux-spectrum, γ = (1 − υ²/c²) ^− 1/2 and where υ₀ is the velocity of a cosmic ray of energy E₀.

Given a steep initial j₀(E), the approximate effect of the Alfvén waves and Lorentz Force is to shift the cosmic ray spectrum, and thus the ionization rate (see next section), from ζ(N_H) to ζ(5N_H), so that the ionization rate decreases more strongly with column.

For cosmic ray flux-spectra that are not very steep (m < 2 for j ∝ E ^− m), the shift is less extreme. Of course, for the full description of the relationship of ζ to the column density, Eq. (4) must be calculated for E < E₀.

2.5. The column-dependent ionization rate

The value of ζ(N_H) is calculated by integrating the product of the flux-spectrum from Eqs. (4)–(5), j_IC(E,N_H), and σ_i,H(E) from Eq. (1), as a function of “depth” N_H into a cloud, with various correction factors: (6)The factor of 5/3 (Spitzer & Tomasko 1968; Dalgarno et al. 1999) takes into account the additional ionization caused by secondary electrons, while the factor of 1.8 accounts for ionization due to α particles (He ^+ 2). These particles are the second most important source of ionizing cosmic rays (ζ_α ≈ 0.8ζ_p). By comparison, relativistic electrons, the third most important ionizing source, have little effect: ζ_e ~ ζ_p/100.

Three different functions for ζ(N_H) have been calculated, based on the cosmic ray flux-spectra in Fig. 1, which are chosen to be widely divergent below E ≈ 500 MeV to account for the uncertainty in the low-energy region (Hayakawa et al. 1961; Spitzer & Tomasko 1968; Nath & Biermann 1994). Analytical expressions for ζ(N_H), used in the models below, are: and are valid for 10²⁴ cm^-2 cm^-2. These analytical expressions do not seem to change significantly for 100 cm^-3  < n < 10⁶ cm^-3 and 5 K  < T < 1000 K, beyond which the effects of the density and temperature on magnetic field effects becomes significant.

The results are depicted in Fig. 2, in terms of the visual extinction between the cloud and the UV source (A_V). We determine A_V ≈ 4.3 × 10^-22N_H, using Q efficiencies from Laor & Draine (1993) with a grain distribution (in terms of the “radius” of the grain, a) of n ∝ a^-3.5 with r_min = 5 nm and r_max = 1 μm. With these assumptions, the analytical expressions for ζ(A_V) are: These expressions are later referred to as “mid-range” and “high-range” values, respectively.

The wide range of the ionization rate demonstrates the importance of low-energy cosmic rays, especially at low N_H or A_V. Other calculations of ζ(N_H) have been performed, either for high column densities (>10²⁴ cm^-2) where low energy cosmic rays do not penetrate (Umebayashi & Nakano 1981; Finocchi & Gail 1997), without consideration of the magnetic field (Padovani et al. 2009), or in regions where there are no ionization losses (Padoan & Scalo 2005). Recently, Padovani & Galli (2011) have incorporated the effect of magnetic mirroring, whereas we have treated the effects of Alfvén waves on cosmic ray streaming. The results in this paper suggest that Alfvén waves may have a more substantial effect on ζ, with factor of  ~10 impact on ζ at certain N_H for Alfvén waves versus a factor of  ~2–4 impact on ζ from magnetic mirroring. Ultimately, a robust magnetohydrodynamics simulation of cosmic ray transport would be necessary to determine what magnetic field effects have the most significant impact on cosmic ray penetration.

In our study below, we determine the effect of four different functions for ζ on the chemistry in the Horsehead nebula. We consider the ζ(N_H) functions based on flux spectra from Nath & Biermann (1994) and Hayakawa et al. (1961), as well as constant values for ζ of 10^-15 s^-1, and 10^-17 s^-1, the latter of which is effectively the ζ(N_H) derived from the spectrum of Spitzer & Tomasko (1968).

Fig. 2 The results of the one-dimensional Monte Carlo model for ζ described in Sect. 2 in terms of AV. The solid red, dashed green, and dotted blue lines derive from the flux-spectra of Nath & Biermann (1994), Hayakawa et al. (1961), and Spitzer & Tomasko (1968), respectively. These lines fit the averaged result of dozens of iterations of the Monte Carlo model. The results from a single Monte Carlo run using the flux-spectrum of Hayakawa et al. (1961) are included (pink dotted) in order to show error.

3. Modeling the Horsehead nebula

We have used the PDR model of Morata & Herbst (2008) with the OSU 03/2008 gas-phase network¹. This network is a purely gas-phase one that treats the PDR as a semi-infinite series of slabs with the radiation source impinging on one edge. It does not account for freeze-out or any surface chemistry, aside from a simple approximation for H₂ formation on grains and selected ion recombination processes. Radiative transfer and self-shielding of H₂ and CO (Draine & Bertoldi 1996; Lee et al. 1996b) are calculated in progression, starting with the slab at the edge. The chemistry is solved with a time-dependent gas-phase kinetics model for each slab. This model, like our model for cosmic rays, is one-dimensional (1D). Because cosmic rays are thought to stream in from all sides, the effects of the geometry are mostly lost in this model. However, even with cosmic rays streaming in from all angles, low energy cosmic rays will dominate at the edge, and will be absent from the center. The average value of ζ at a slab near the edge will be close to the value determined from the 1D Monte Carlo model. Because the majority of slabs near the center will not have low-energy cosmic rays, the average ζ near the center also be close to the 1D value. Thus the average value of ζ at different slabs of the cloud will be close to the 1D values we use for ζ found in Fig. 2.

Following Pety et al. (2005), we compare, when possible, observations with model results for three regions at different optical extinctions (A_V) from the edge of the Horsehead PDR. These are the IR-edge (A_V = 1.56 ± 0.73), IR-Peak (A_V = 4.55 ± 1.7), and the Cloud (A_V = 11.7 ± 4.1). The error bars in A_V are based both on the beam size of the observations and uncertainty in the density profile of the cloud, as discussed in the next section. We determine the error in fractional abundance by taking the ratio between the observed column density of the species and the error in that column density, both from Pety et al. (2005).

3.1. Physical conditions and initial chemical abundances

The density profiles used are taken from Habart et al. (2005). The temperature profile is calculated from thermal balance (Le Petit et al. 2006). Cosmic rays heat the interstellar medium through the thermalization of secondary electrons and photons (Field et al. 1969; Glassgold & Langer 1973). Thermal heating by cosmic rays begins to dominate at A_V > 3, but the thermal impact of different cosmic ray ionization rates is not very significant until ζ > 10^-16 s^-1. Since even the highest ζ(N_H) drops to  ≈ 10^-16 s^-1 at the Cloud region, the temperature difference here between the high ζ(N_H) and ζ = 10^-17 s^-1 is only about 4 K. The density and temperature profiles are shown in Fig. 3.

The gas density increases with spatial distance into the nebula as a power law with an exponent β (Habart et al. 2005), which in terms of column density can be written as: (11)

where β ≥ 1 is a dimensionless constant used to parameterize the steepness of the number density, n_H,0 = 2 × 10⁵ cm^-3, x₀ = 0.02 pc is a length scale, and N_H,0 = (1 + β)^-11.23 × 10²² cm^-2 is the column density at a depth of x₀. For our analysis, we show the results for β = 1, and discuss results for both β = 1 and β = 4. The steeper density gradient impacts the UV photon flux and the resulting thermal balance. There are different total densities for the IR-edge and IR-peak regions. The difference in UV penetration, temperature and density at different values of A_V noticeably impacts the chemistry. The cosmic ray ionization, however, is not significantly altered by the density gradient, because for the ranges of density of 100 to 10⁵ cm^-3, ζ is column-dependent, but not density dependent.

Other densities and density profiles have been proposed. Pety et al. (2005) used several uniform number densities and profiles, while Goicoechea et al. (2009b) proposed a slowly changing piecewise function for the density, with three sections instead of two, reaching 2 × 10⁵ cm^-3 at A_V ≈ 5 instead of A_V ≈ 1.0, as used here. Until the number density is better determined, significant uncertainties in the extinction at a given angular depth will persist.

Fig. 3 The temperature (dashed line) and density (solid line) profiles as functions of visual extinction with ζH,Nath. The density profile is in the form of Habart et al. (2005), our Eq. (11), with β = 1. The temperature is from thermal balance (Le Petit et al. 2006). At AV = 10, ζ ≈ 10-16 s-1, which raises the temperature by  ≈ 4 K at the center compared to a ζ of 10-17 s-1.

The UV radiation field impinging on the Horsehead surface has been a topic of much discussion and uncertainty (Anthony-Twarog 1982; Zhou et al. 1993; Abergel et al. 2003). Values of χ = 30 to χ = 100 in Draine units (Draine 1978) have been proposed. We use χ = 60, because this is the most commonly used value for the Horsehead PDR. The external UV field is important to the chemistry only for the IR-edge. For the IR-peak and the Cloud regions, cosmic rays are the primary ionizing and photochemical agent.

Fig. 4 Fractional abundances of C2H, c-C3H2, C4H, and HC3N as functions of AV. The boxes represent observations with error bars, and the lines are the model results for ζ = 10-17 s-1 (blue dashed), ζ = 10-15 s-1 (pink dotted), and, from Fig. 2, the mid-range ζ(NH) (green dashed) and high-range ζ(NH) (red solid).

The initial chemical abundances used for the Horsehead PDR (Lee et al. 1996a; Morata & Herbst 2008) are listed in Table 2 and represent abundances for a dark cloud prior to the onset of a nearby star. These abundances comprise observed values for small (less-than-six-atom) species in TMC-1, as well as calculated early-time values from Smith et al. (2004) for atoms and small molecules that have not been observed, based on so-called “low-metal” elemental abundances.

In addition to these initial abundances, we also investigated cases with much higher elemental abundances of sulfur, based primarily on the analysis of CS and HCS⁺ by Goicoechea et al. (2006), who place the total elemental sulfur abundance with respect to n_H at 3.5 × 10^-6. On the other hand, Teyssier et al. (2004) used a value of  [S]  ~ 10^-7, similar to the low-metal value used in this part of the paper. To determine the effect of raising the sulfur abundance, we utilized elemental abundances for sulfur, relative to hydrogen, of 10^-6 and 10^-5, starting primarily from the neutral atomic form.

The abundances are calculated from time t = 0 to steady state (t = 5 × 10⁶ yr). Since the age of the Horsehead nebula is not well-determined, values from 10⁴ − 10⁶ yr have been considered (Morata & Herbst 2008). We focus only on the time of 10⁵ yr, because in general the calculated results are closest to observational values at this time. We also use this time because it is a reasonable age for a molecular cloud, given its size and velocity gradient (see Pound et al. 2003). Time-dependence was investigated by Morata & Herbst (2008) albeit with a different density profile from what is used here. They found that at times between 10⁵ yr and steady-state, the abundances of carbon chain species in the Cloud region become sharply lower, as is found in standard cold dark clouds. They also investigated times as early as 10⁴ yr, at which time the abundance profiles are flatter. Our calculations for carbon chain species have reached steady state by 10⁴ years for A_V < 5. For A_V > 5, our calculations confirm their findings.

Fig. 5 Relative abundances of HCO+, HCO and the ionization fraction as functions of AV. The boxes represent observations with error bars, and the lines are the model results for ζ = 10-17 s-1 (green dashed), ζ = 10-15 s-1 (pink dotted), and, from Fig. 2, the mid-range ζ(NH) ( blue dotted) and high-range ζ(NH) (red solid).

In Figs. 4 to 6, we show the calculated abundances of various molecules as continuous functions of visual extinction with observed values in boxes to delineate the uncertainties in both abundance and A_V. The calculated abundances are plotted with two fixed values of ζ: 10^-17 s^-1 and 10^-15 s^-1, as well as with two column-dependent ionization rates depicted in Fig. 2: the mid-range ζ(N_H) (dashed green line), and the high-range ζ(N_H) (solid red line). The fixed value of ζ = 10^-17 s^-1 is equivalent to the lowest-range ζ(N_H) in Fig. 2, which utilizes only high-energy protons. Neither of the two fixed values for ζ is likely to be physically reasonable; the low value can pertain to the inner Cloud region but is less likely to pertain to a region near the edge, where at least some low-energy cosmic rays exist, while the high value is more likely to pertain only to the edge of the PDR. Unless specified, the low elemental abundance of sulfur is utilized.

3.2. Results: C₂H, c-C₃H₂ and C₄H

Hydrocarbons are not direct tracers of ζ; nevertheless, an enhanced ζ at the surface of the Horsehead nebula may help to explain the high abundances of these small hydrocarbons at the edge. C₂H, c-C₃H₂ and C₄H are formed by a complex network of reactions, linked at least partially to the cosmic ray ionization rate via several sequence of reactions based on C and C⁺. The sequence involving neutral atomic C starts with the reactions: and CH⁺ initiates a series of chemical reactions that eventually results in C₂H, c-C₃H₂ and C₄H via recombination with electrons. The C⁺ ion is produced in three ways depending on physical conditions: at low extinction (A_V < 2.5), it is formed principally by photoionization, and can reach a fractional abundance as high as 10^-4, whereas at high extinction (A_V > 4.5) it is formed less efficiently by the reaction between He⁺ and CO. In the middle region (2 < A_V < 5), secondary photons from cosmic rays form a large amount of the C⁺. Once produced, it can radiatively associate with H₂ to form the CH ion, which initiates a series of reactions similar to those initiated by CH⁺ (Herbst & Millar 2008). Because of these alternate pathways, small hydrocarbons may not be as sensitive to ζ very close to the edge or deep within the Horsehead PDR. Regardless, our robust chemical network allows us to explore in detail the effect of a column-dependent ζ on the Horsehead nebula.

The model abundances for C₂H, c-C₃H₂, and C₄H vs. A_V can be found in Fig. 4, where observed abundances with estimated uncertainties are plotted as boxes for the three regions studied: the IR-edge, the IR-peak, and the Cloud. For C₂H, our use of temperature and density profiles seems to account for the observed abundance at the IR-edge, regardless of the value of ζ, probably because C₂H formation is so dependent on photon effects at the edge. The results diverge for the IR-peak and Cloud, where the high-range ζ(N_H) and ζ = 10^-15 s^-1 seem to do better than the other two choices of ζ. In the IR-peak, the abundances obtained with the high-range ζ(N_H) and ζ = 10^-15 s^-1 come within a factor of  ≈ 5 of the observed value, and are closer still for the Cloud region.

Fig. 6 Relative abundances of HCS+ and CS as a function of AV. The boxes are the observations with error bars, and the lines are the model results with  [S]  = 10-5 (red),  [S]  = 10-6 (green) and  [S]  = 7.2 × 10-8 (blue), all using the high column-dependent ζ from Fig. 2. The solid lines use a rate for S + H+ → S+ + H of 1.3 × 10-9 cm3 s-1 (Prasad & Huntress 1980) and the dashed lines use a rate of 1 × 10-14 cm3 s-1.

For c-C₃H₂, and for C₄H, none of the four plots comes particularly close to the observed values at the center of the IR-edge, although the curves obtained with the high-range ζ(N_H) and ζ = 10^-15 s^-1 graze the lower portion of the observation box for C₄H. This discrepancy suggests that, though a high surface ζ is important, there are likely other factors that must be taken into account, such as PAH fragmentation (Pety et al. 2005). For the IR-peak region, the high-range ζ(N_H) and ζ = 10^-15 s^-1 models lead to results that graze portions of the observational boxes for both species, with the others models exhibiting much too low an abundance. Finally, for the Cloud region, the high-range ζ(N_H) and ζ = 10^-15 s^-1 models do quite well for C₄H, and c-C₃H₂. while the lower ionization models show reasonable agreement only for the latter.

It would appear that, on balance, the results obtained with the high constant ζ and the high-range column-dependent ζ are closer to observation in most instances for these three hydrocarbons. To further distinguish between these two sets of results, we focus on the abundance ratios between IR-peak and Cloud regions for the three carbon-chain species. The ratios are taken at the visual extinctions where the models agree best with the observations, and are listed in Table 3. The reason for taking these ratios is that we can better compare results between a fixed and a column-dependent ionization rate in this manner. These ratios are examined only as a way to distinguish between a constant and a column-dependent ζ, and their use beyond this function is severely limited. For example, the C₂H emission attributed to the Cloud region may be from the FUV illuminated surface (for an analogous example involving HCO, see Gerin et al. 2009). It is likely that the observed ratios will change and will be far better constrained when the Horsehead nebula is explored at higher angular resolution.

For C₂H and c-C₃H₂, the ratios are much closer to observation for the column-dependent ζ(N_H) than for ζ = 10^-15 s^-1. In both of these cases, the ratios from the ζ(N_H) model are within a factor of 2 of the observed ratios. For ζ = 10^-15 s^-1, model ratios disagree by a factor of 5–7. In the case of C₄H, the ratio from the constant ζ agrees slightly better with observations than for ζ(N_H), although the ratios of both models are close to observation. Also, examining the C₄H abundances from Fig. 4, it is evident that, within the Cloud region, the results from ζ(N_H) are much closer to observation than the results from ζ = 10^-15 s^-1. In summary, as well as being unphysical, the results from a model with a constant ζ = 10^-15 s^-1 do not agree as closely as the results from a model with a column-dependent ζ(N_H).

For β = 4 at t = 10⁵ yr, the results are not significantly changed for the IR-edge, and the model underestimates the small hydrocarbon abundances for the IR-peak region by about an order of magnitude. The reasons for this seem to be as involved as the hydrocarbon chemistry. The most significant factor is that the production of these hydrocarbons at a higher density requires a higher ζ, and for β = 4, the density is much higher at the IR peak than for β = 1. Also, with the steeper density gradient, photons are more effective at ionizing and dissociating at the IR-Edge, but fall off more abruptly at higher A_V. The difference in C⁺ formation by photons between β = 1 and β = 4 density profiles is a factor of three, and only present at A_V < 2.5.

It should finally be mentioned that thermal balance from photons depends somewhat on the density, and so the temperature profiles with β = 1 and β = 4 are different. At A_V = 0.001, the gas-phase temperature for β = 1 is about 300 K, whereas for β = 4, T ≈ 600 K. The gas-phase temperatures for the two density profiles converge at A_V = 1, and this undoubtedly has some impact on the chemistry. It should be emphasized that ζ(N_H) is similar for the steep and gradual density gradients.

3.3. Results: HC₃N, HCO⁺, HCO and the electron fraction

Only one line of the carbon-chain species HC₃N has been detected, and this with a very large beam-size (Teyssier et al. 2004). We follow Teyssier’s tabulated value for A_V, and treat the emission as originating in the IR-peak, though there is some uncertainty about the origin of this emission. The four models all under-produce the observed abundance of HC₃N by a little less than an order of magnitude or more, as can be seen in Fig. 4, with the models with the high-range ζ(N_H) and the fixed ζ = 10^-15 s^-1 coming closest.

Cyanoacetylene (HC₃N) is not as dependent as the other species on cosmic ray ionization for much of the range of visual extinction. Two reactions primarily lead to its formation: At the edge, the first reaction is directly related to ζ through C₂H₂, but the second reaction involves C₃H₂N⁺, the formation of which is not strongly dependent on ζ. In the Cloud region, the situation is reversed: C₂H₂ is less dependent on ζ, and C₃H₂N⁺ is then closely linked with cosmic ray ionization. Because of the two channels for HC₃N we expect less dependence on A_V except in the middle range: 1 < A_V < 5. The results, shown in Fig. 4, roughly bear this out. Interestingly, both the observed and calculated abundances for HC₃N are much lower than the initial value, which is taken from the TMC-1 abundance. The discrepancy with the Cloud value, over three orders of magnitude, is especially large and very different from the analogous cases for the hydrocarbons in Fig. 4.

Figure 5 contains the observations and model results for HCO⁺ and HCO. Since HCO⁺ is optically thick, the carbon-13 isotopologue was used for observations. H¹³CO⁺ was observed in emission at  ≈ 40′′ from the PDR edge (Gerin et al. 2009), corresponding to an A_V ≈ 10, which is essentially the Cloud region. Following the analyses of Gerin et al. (2009) and Goicoechea et al. (2009b), we determined the abundance of HCO⁺ from H¹³CO⁺ by assuming ¹²C/¹³C = 60. A faint emission feature attributed to H¹³CO⁺ was also seen at  ≈ 10′′ from the PDR edge, corresponding to an A_V ≈ 2 with our density profile, and so lies essentially at the IR-edge.

In the immediate neighborhood of A_V = 2, however, none of the models produces enough HCO⁺, but the increase in abundance with increasing extinction is steep and by A_V = 3, all but possibly the ζ = 10^-17 s^-1 model produce a comparable result to what is observed at the IR-edge. Goicoechea et al. (2009b) did much better fitting the HCO⁺ abundances at the edge by including PAH’s. They also modeled profiles for H¹³CO⁺ and DCO⁺.

For the Cloud value, all models are in reasonable agreement with observation for HCO⁺, coming within factors of 2–5 of the observed abundance. The formation of HCO⁺ by cosmic rays is very direct at high extinction; in regions where UV photons cannot penetrate, it is almost solely the product of the destruction channel for protonated molecular hydrogen with carbon monoxide. At the IR-edge, however, the UV driven formation by the reactions H₂ + CO⁺ and H₂O + C⁺ dominates. In all regions, HCO⁺ is destroyed mainly by recombination.

For neutral HCO, all model results are too low by an order of magnitude or more at both the IR-edge and Cloud regions, even with the relatively fast reaction between CH₂ and O (from Gerin et al. 2009). Our results disagree with the model results from Gerin et al. (2009) and Goicoechea et al. (2009b) partly because the Meudon reaction network includes a formation mechanism absent in the OSU network, the photodissociation reaction (15)where hν represents an external UV photon. This reaction is also discussed in Gerin et al. (2009). Including this reaction enhances the HCO abundance by a factor of 5 in the PDR, bringing the HCO abundance within an order of magnitude of the observed value.

The ionization fraction, f(e ⁻ ), is a measure of elemental abundances, ionization rate, density, and chemistry, as well as a constraint on the coupling of the magnetic field to the matter in the cloud. The ionization fraction from our models, as shown in Fig. 5, ranges from  ~10^-4 in the PDR to  ~10^-8 in the Cloud region. This range of fractions agrees generally with the profile in Goicoechea et al. (2009b, their Fig. 4). Their inferred profile for the ionization fraction would favor the mid-range ζ(N_H) from the cosmic ray flux-spectrum of Hayakawa et al. (1961).

For the steeper density profile with β = 4 and the high ζ(N_H), our results are somewhat different. The HCO⁺ abundances are not significantly changed, and the modeled HCO abundances increase by a factor of two in the IR-Edge and IR-Peak regions (at 10⁵ yr). Significantly, our calculated abundance of HC₃N comes into good agreement with the Cloud region observation; it is a factor of 3 higher than the observed abundance at t = 10⁵ yr.

3.4. Tabulated abundances

Calculated fractional abundances (with respect to n_H) obtained with the standard elemental abundances are listed for more than twenty species in Table 4, including both observed and undetected molecules. The calculated results are for a time of 10⁵ yr and pertain to the center points of the IR-edge, IR-peak, and Cloud regions (Pety et al. 2005), for which observational results are also shown, when available. Some of the tabulated abundances, HOC⁺ especially, seem to be possible tracers for the cosmic ray ionization, because their fractional abundance becomes more dependent on the extinction when ζ depends on column density, than when ζ is a constant value.

In this table, we consider only the model with the high-range ζ(N_H), because it is evident that, at least for carbon-chain species, use of this column-dependent ζ leads generally to better agreement with observations than models with lower ionization, and it is more physical than the constant high-ionization model. Also, we do not include the case of the steeper density profile in this table. Predictions are discussed below in Sect. 3.6.

3.5. The sulfur-rich case

We considered sulfur-bearing species, both with the standard initial abundances, and also for a sulfur-rich environment. We found that the higher the elemental sulfur (up to a relative abundance of 10^-5), the closer the model matches observations for sulfur-bearing molecules. Our results and those of Goicoechea et al. (2006) for the chemistry and radiative transfer agree very well.

The results for the observed sulfur-bearing species CS and HCS⁺ vs. A_V at 10⁵ yr can be found in Fig. 6 as a function of the sulfur elemental abundance. There are two sets of curves, depending upon the rate coefficient for the charge-exchange reaction which can affect the abundances of CS and HCS⁺ at low sulfur abundances. This reaction has a listed rate coefficient of 1.3 × 10^-9 cm³ s^-1 (Prasad & Huntress 1980) but a more likely value of 1 × 10^-14 cm³ s^-1 has been calculated².

The agreement attained by increasing the elemental abundance,  [S] , to 10^-5 comes at a cost: at 10⁵ yr, all the carbon-bearing species in this scenario are reduced by up to a factor of 10 except at the IR-edge. This effect is most severe in the Cloud region. This depletion occurs in part because the high sulfur abundance destroys hydrocarbons by reactions with S⁺ and also with S at higher extinctions and because of the increased fractional ionization. The depletion of carbon-bearing species worsens agreement for all observed species except HCO⁺, which is brought to within a factor of 2 of observation in the Cloud region.

This problem may suggest that a more realistic gas-phase sulfur elemental abundance for the Horsehead nebula should lie somewhere around 10^-6, in agreement with Goicoechea et al. (2006). The abundances of observed and predicted molecules with  [S]  = 10^-6 are in Table 5 for the same species as listed in Table 4. Even with this intermediate sulfur abundance, the calculated abundances of carbon chain species in particular are lowered considerably compared with the corresponding values in Table 4, leading to worse agreement with observation.

3.6. Some predictions

A high column-dependent ζ brings with it implications for chemistry in the Horsehead PDR. This column-dependent ζ varies from  ≈ 2 × 10^-16 s^-1 at the IR-edge to  ≈ 7 × 10^-17 s^-1 in the Cloud region and so leads to profiles distinctive from models with fixed ionization rates, as can be seen for carbon-chain species in Figs. 4 to 6.

Also, other molecules are predicted to be in amounts in principle observable, and these are listed among the species in Tables 4 and 5. Because our ζ(N_H) produces reasonable abundances of C₄H and HC₃N in selected regions with a low elemental abundance of sulfur, we would also expect to observe, albeit with some difficulty, the more complex carbon-chains C₆H and HC₅N, based on our predictions for these regions. In addition, the molecule HCN should definitely be present in observable quantities, especially in inner regions, and its isomer, HNC, should also be observed with a ratio HCN/HNC ≈ 1. We predict ammonia in observable quantities at A_V > 4, for the low-sulfur case.

Given the observations of high amounts of the reactive molecular ions OH⁺ and H₂O⁺ in many molecular objects (Gerin et al. 2010; Gupta et al. 2010), it would be useful to consider predicted abundances of these species. Our model predictions for OH⁺, H₂O⁺ and H₃O⁺ in the Horsehead nebula are contained in Tables 4 and 5. These predictions show low abundances for the first two ions that are rather independent of which of the three regions we consider. The basic problem is the low abundance of atomic hydrogen except at the border of the PDR (Neufeld et al. 2010). Even at the IR-Edge, H₃O⁺ is more than an order of magnitude higher than either OH⁺ or H₂O⁺, though none of these species should be sufficiently abundant to be detected. In the Cloud Region, where the electron density is at the low level of a cold dark cloud, H₃O⁺ is depleted rather slowly by reactions with electrons, and should achieve a high enough column to be detectable.

4. Discussion

We have modeled the Horsehead nebula as a PDR with time-dependent gas-phase chemistry using a column-dependent cosmic ray ionization rate ζ(N_H), as well as the temperature and density profiles of Habart et al. (2005). At a cloud age of 10⁵ yr, the incorporation of a high ζ(N_H) improves agreement between model and observation for the small carbon-bearing molecules HCO⁺, HC₃N, C₂H, c-C₃H₂, and C₄H compared with a more standard constant ionization rate. With a higher abundance of elemental sulfur than our standard value, the results for small sulfur-bearing species are improved, but at the expense of our calculated values for carbon-chain species. There are also predictions of abundances and profiles for other species, some not yet observed in the Horsehead nebula, which should be in principle observable, including HCN, HNC, NH₃, C₆H, HC₅N, and H₃O⁺. Some of these predictions are strongly affected, however, by an increase in the assumed sulfur elemental abundance.

Our results for c-C₃H₂ and C₄H (but not for C₂H) also indicate that the fracturing of PAH’s may play an important role in the production of these molecules towards the edge of the PDR, but our model does not incorporate the effects of PAH’s. Strong aromatic emission, observed by Compiègne et al. (2007), poses some problems, however, for the hypothesis that PAH fracturing is the source of small hydrocarbons. These authors claim a high concentration of neutral PAH’s in the HII region, which suggests that PAH’s may endure the radiation at the IR-edge, instead of breaking apart into the observed hydrocarbons.

The detailed form of the calculated abundance profiles in Figs. 4 through 6 cannot be observed because observations up to the present lack sufficient resolution, and because the density profile is not well-determined. With the advent of the Atacama Large Millimeter Array (ALMA), the estimated increase in angular resolution, to  ~0.1′′ (Wootten 2003), should allow us to observe the form of these abundance profiles, so as to better determine the initial flux-spectrum for cosmic rays for the Horsehead nebula.

It appears, from Indriolo et al. (2010), that there is some environmental influence on the low energy flux of cosmic rays. It would be of great interest to not only examine the Horsehead nebula at greater angular resolution, but to also observe and model other PDR’s such as the Orion Bar, IC-63, L1688-W, and portions of Sgr B2 to determine how the low energy cosmic ray flux varies in our Galaxy. Sgr B2 is of special interest given the high values of ζ inferred from H observations in this region (Oka et al. 2005). Given the strong dependence of ζ on the path cosmic rays travel, it is very likely that the low-energy cosmic ray flux at cloud edge will be object-dependent.

Acknowledgments

We would like to thank Tom Millar, Ben McCall, and Nick Indriolo for helpful discussions about both the cosmic ray ionization rate and its implementation in chemical models. We also thank Ewine van Dishoeck, John Black, Tom Cravens and Andy Skilling for helping us properly consider the magnetic field. E.H. acknowledges the support of the National Science Foundation for his astrochemistry program through grant AST-0702876, and his program in chemical kinetics through the Center for the Chemistry of the Universe. He also acknowledges support from NASA NAI for studies in the evolution of pre-planetary matter and from JPL for Herschel studies. E.R. acknowledges the support of the Observatoire de Paris and the Programme National du CNRS-INSU “PCMI”.

References

All Tables

All Figures

	Fig. 1 Three different cosmic ray flux spectra, taken from Hayakawa et al. (1961) (dashed line), Spitzer & Tomasko (1968) (dotted line), and Nath & Biermann (1994) (solid line).
In the text

Fig. 2 The results of the one-dimensional Monte Carlo model for ζ described in Sect. 2 in terms of AV. The solid red, dashed green, and dotted blue lines derive from the flux-spectra of Nath & Biermann (1994), Hayakawa et al. (1961), and Spitzer & Tomasko (1968), respectively. These lines fit the averaged result of dozens of iterations of the Monte Carlo model. The results from a single Monte Carlo run using the flux-spectrum of Hayakawa et al. (1961) are included (pink dotted) in order to show error.

In the text

	Fig. 3 The temperature (dashed line) and density (solid line) profiles as functions of visual extinction with ζH,Nath. The density profile is in the form of Habart et al. (2005), our Eq. (11), with β = 1. The temperature is from thermal balance (Le Petit et al. 2006). At AV = 10, ζ ≈ 10-16 s-1, which raises the temperature by  ≈ 4 K at the center compared to a ζ of 10-17 s-1.
In the text

	Fig. 4 Fractional abundances of C2H, c-C3H2, C4H, and HC3N as functions of AV. The boxes represent observations with error bars, and the lines are the model results for ζ = 10-17 s-1 (blue dashed), ζ = 10-15 s-1 (pink dotted), and, from Fig. 2, the mid-range ζ(NH) (green dashed) and high-range ζ(NH) (red solid).
In the text

	Fig. 5 Relative abundances of HCO+, HCO and the ionization fraction as functions of AV. The boxes represent observations with error bars, and the lines are the model results for ζ = 10-17 s-1 (green dashed), ζ = 10-15 s-1 (pink dotted), and, from Fig. 2, the mid-range ζ(NH) ( blue dotted) and high-range ζ(NH) (red solid).
In the text

Fig. 6 Relative abundances of HCS+ and CS as a function of AV. The boxes are the observations with error bars, and the lines are the model results with  [S]  = 10-5 (red),  [S]  = 10-6 (green) and  [S]  = 7.2 × 10-8 (blue), all using the high column-dependent ζ from Fig. 2. The solid lines use a rate for S + H+ → S+ + H of 1.3 × 10-9 cm3 s-1 (Prasad & Huntress 1980) and the dashed lines use a rate of 1 × 10-14 cm3 s-1.

In the text