Abstract

The transformation of the energy dependence of the cosmic ray proton flux in the keV to GeV region is investigated theoretically when penetrating inside molecular clouds (Av>5 mag). The computations suggest that energy losses of the cosmic ray particles by interaction with the matter of the molecular cloud are principally caused by the inelastic (electronic) interaction potential; the transformed energy distribution of energetic protons is determined mainly by the column density of the absorbing medium. A cutoff of the cosmic ray spectrum inside clouds by their magnetic fields is also phenomenologically taken into account. This procedure allows a determination of environment-dependent ionization rates of molecular clouds. The theoretically predicted ionization rates are in good agreement with those derived from astronomical observations of H3+ absorption lines in the spectrum of the cloud connected with the Herbig Be star LkH𝛼 101.

1. Introduction

The physical characteristics of interstellar clouds such as temperature, pressure, and radiation balance as well as energy deposition from cosmic ray (CR) particles and ultraviolet photons (UV) depend strongly on the state of hydrogen (atomic versus molecular, neutral versus ionized) [4]. Interstellar clouds can be categorized thoroughly into three classes based on their size, gas densities, temperature, and ionization fraction [47]: (i) diffuse, (ii) translucent, and (iii) molecular clouds. Diffuse cold clouds are characterized by typical number densities of less than 102 cm−3 and temperatures of about 100 K. Those diffuse regions in which hydrogen is present mainly in its neutral, atomic form are called HI clouds [4]. Interstellar regions in which the gas exists predominantly in molecular form, that is primarily molecular hydrogen (H2), are referred to as molecular clouds. These areas hold average number densities of 103-104 cm−3 and temperatures of about 10 K [8]. As intermediate objects between the low-extinction diffuse clouds (Av < 1) and the optically thick molecular clouds (Av > 5), van Dishoeck and Black [9] introduced the class of translucent clouds with visual extinction spanning the range Av~1–5 mag. The interior of these clouds are entirely blocked from the external UV field at 13.6 eV (the ionization energy of atomic hydrogen) and beyond; recall that the column density along the line of sight, 𝑁 (in cm−2), is related to AV (in mag) by 𝑁=21021Av [10]. As a matter of fact, each molecular cloud has a translucent layer. Here, the edges of each molecular cloud with average concentration of 𝑛~103 cm−3 and, in the ultimate sense, the edge of every clump with 𝑛105 cm−3 must be considered as a translucent region [11]. However, although UV does not penetrate into the deeper layers, energetic CR particles can induce an ionization and destruction of molecular hydrogen in the interior of molecular clouds (1); here, p(𝐸) presents a proton of the CR with a kinetic energy before (𝐸) and after (𝐸) the ionization process. Further, the singly ionized hydrogen molecule (𝐻2+) can react with H2 to form 𝐻3+ via(2)p(𝐸)+H2H2++e𝐸+pH(1)2++H2H3++H.(2) Therefore, the measurement of the H3+ column density can be used as a direct probe of the CR ionization rate (𝜁) in dense clouds [12]. This ionization rate presents a critical and fundamental quantity to physical and chemical models of molecular clouds, for which no direct measurement exists, and usually it is chosen as a free parameter [13]. Furthermore, the H3+ cation is expected to be responsible for initiating ion-molecule reactions that lead, for instance, to the formation of simple hydrides (water, ammonia, methane) and possibly to higher order, saturated hydrocarbons [14]; unsaturated hydrocarbons and their radicals, on the other hand, can be synthesized via rapid neutral-neutral reactions [6]. Prior to the direct identification of H3+ in infrared absorption against continuum radiation emanating from background protostars (W33A and GL2136) [15], molecular ions like HCO+ and N2H+ have been utilized to estimate the ionization rate in molecular clouds [16]. On the other hand, the number density of H3+ in the steady-state, 𝑛(H2), is directly related to the CR ionization rate, 𝜁 (in s−1), via (3), with the number densities of molecular hydrogen, 𝑛(H2), and of carbon monoxide, 𝑛(CO) [17]H𝜁𝑛2𝑘COH𝑛3+𝑛(CO).(3) Here, it is assumed that H3+ is produced solely via (1-2) followed by reaction with any gas-phase constituent of the cloud, X, with the exceptions of He, O2, and 𝑁H3++XHX++H2.(4) Holding fractional abundances between 10−4 and 10−5, CO is the dominant reactant in the molecular cloud; the ion-molecule reaction rate of CO with 𝐻3+ has no entrance barrier and is very fast with gas-kinetics limits (𝑘CO=1.8109 cm3s−1 [17]. Considering that the ratio of the number densities, 𝑛(CO)/𝑛(H2), is approximately constant in dark clouds with 1.5104 [17, 18], Equation (3) reduces to (5). Assuming a constant value of 𝜁=31017 s−1, Geballe [17] estimated number densities of 𝑛(H3+) 104 cm−3𝑛H3+3.71012𝜁.(5) However, it is clear that 𝜁 can depend on the cloud’s column density, 𝑁(H2). According to observations, the ionization rate varies between 10−18 and 10−16 s−1 in dense and diffuse clouds, respectively [13]. Comparing chemical models of molecular clouds with actual astronomical observations, Caselli [13] derived an empirical relationship between the number density dependent ionization rate 𝜁(𝑛) and the number density of H2 with an empirical parameter 𝑞0.5(6)𝐻𝜁(𝑛)𝑛2𝑞.(6) Because the cross-section, 𝜎 (cm2), of the ionization of H2 via reaction (1) is well known [1], then it is possible to relate the CR flux, 𝐹s (cm−2s−1eV−1), in an energy interval Δ𝐸 (eV) to the ionization fraction and the cross-section via (7) [19]𝐹𝑠Δ𝐸=𝜁/𝜎.(7) For instance, considering 𝜁=1-51017s1 [20], 𝜎(5MeV)=11017 cm2 (Figure 1, [1]), this yields 𝜁/𝜎=1-5cm2s1. Hence, a detailed, theoretical determination of a function 𝜁[𝑛(H2)] demands a thorough incorporation of self-consistent cloud models [13] to include depth-dependent ionization rates.

In this paper, I present a numerical approach to compute 𝜁(𝑟), in a molecular cloud provided that the energetic spectral flux of the CR particles outside of the given cloud and a radial distribution of the number density of H2 inside the cloud are known. It was shown long ago that molecular clouds can magnetically reflect low-energy CR particles by a screening mechanism involving the generation of Alfven waves [21, 22]. This means that the CR fluxes may be cut off at low energies and one should be able to choose a corresponding limiting value of energy (Section 2.1). Further, the column density of H3+, 𝑁(H3+), can be calculated as the obvious consequence of 5, via (8) and can then be compared with observational column densities of H3+. This enables one to derive the depth-dependent ionization rate where 𝐷 is the cloud size𝑁H+3=3.7×1012𝐷0𝜍(𝑟)𝑑𝑟.(8)

To compute the environment-dependent ionization rates, it is crucial to extract first the energy-dependent fluxes of the galactic CR particles outside the molecular clouds taking into account a cutting off of the fluxes by magnetic fields (Section 2.1). We will investigate then the energy loss processes of these charged particles upon their penetration inside the molecular cloud structure (Section 2.2) and calculate the change in flux of the CR particles inside the molecular cloud (Section 2.3). There were many earlier works devoted to such calculations (see [2326] and references therein) where only diffuse HI clouds and protostars were considered. These authors used galactic CR fluxes that gave overestimated ionization rates at clouds edges, and clouds magnetic fields have not been taken into account. The aim of this work is to show that the ionization rate inside contracting molecular clouds is strongly dependent on environment and the limiting value of CR energy. The ionization rates in appropriate astrochemical models should be chosen taking into account both CR particles energy loss processes and the cutting off their spectral fluxes by magnetic fields of the clouds.

2. Penetration of Cosmic Ray Particles Inside the Clouds

2.1. The Non-Attenuated Cosmic Ray Spectrum

CR particles—predominantly protons (90%) and α-particles (10%)—are accelerated by magnetic fields in hydrodynamically active regions in which the kinetic energy of ejected stellar mass motions can be transferred magnetically into the energy of the particle, thus reaching up to the TeV region [10]. In the diffuse interstellar medium (ISM), these energetic particles can propagate almost freely without losing a significant fraction of their kinetic energy to the matter in the diffuse cloud [4]. Only protons will be considered here. Very recently, Cooper et al. [2] compiled observational and model data of the CR proton flux in the very local ISM (Figure 2), which is the power-law extension of the CR spectrum to an intersection with the high-energy side of the thermal plasma distribution. It should be noted that the non attenuated particle flux in [2] is different from the usually used distribution [3] (Figure 2). The total proton flux in the energy range 1 eV−10 GeV calculated from [2] is 7.4104 particlescm−2s−1 with an energy deposition of 4.3109 eV cm−2s−1, while on the basis of data from [3] one got 9.6 particles cm−2s−1 and 3.1109 eV cm−2s−1, respectively. In our model, we assume that the CR particle flux in the solar neighborhood is identical to that of in the ISM.

Most importantly, we have to consider that even the non attenuated spectrum of the CR particles can be influenced by Alfven waves (see [21, 22] and references therein). Only particles with energies greater than a minimum energy 𝐸min of about 100 MeV will enter the cloud. A more refined investigation depicted that the value of 𝐸min (10 MeV) is given by (9) where 𝑁 (in cm−3), 𝑅 (in pc), and 𝐵 (in μG) are normalization parameters of the cloud [21, 22]. Here 𝑛(H2), 𝑅c, and 𝐵c are a number density, the radius, and the magnetic field strength, respectively. Observations reported measurements of molecular cloud magnetic field strengths ranging from 0.1 μG to about 6 μG [27]. The typical value of the field strength in the ISM is 4-5 μG [28]; thus, 𝐸min varies between about 10 to 100 MeV for typical molecular clouds𝐸min50𝑁𝑅𝐵1/2,𝑛𝐻2=2103𝑅𝑁,c𝐵=8.2𝑅,c=50𝐵.(9) It is obvious that in order to estimate the environment-dependent ionization rates in molecular clouds systematically, a parameter study is crucial. Here, we compute the ionization rates for two flux distributions of the CR ([2] versus [3]); the influence of each of the distributions with the minimum energy 𝐸min induced by the magnetic shielding is investigated for 𝐸min=10 MeV, 50 MeV, and 100 MeV. Finally, these models are applied at a typical cloud scenario: a molecular cloud with H2 number densities of 𝑛=103 cm−3 and an average size of 10 pc.

2.2. Interaction of Energetic Particles with Matter

The interaction and hence energy loss processes of each charged particle from the cosmic radiation field with the matter of the molecular cloud depend primarily on the kinetic energy of the particles [6, 28] and the physical characteristics of the ISM (density, spatial size, magnetic field). The elementary composition of CR in the keV to GeV energy range for the most important elements is similar to the cosmic elementary abundance. Here, protons (p) and α-particles dominate with an α/p ratio of less than 0.1; the abundance of heavier ions accounts for less than a few percent [29].

When an energetic particle (projectile) progresses through a molecular cloud, it can interact by two principal energy loss mechanisms with the surrounding matter [30]. These are: (i) an interaction of the projectile with the electronic potential (electrons) of the atoms/molecules and (ii) with the nuclei of the atoms/molecules in the interstellar cloud (nuclear interaction/collision; elastic interaction). Inelastic energy loss processes change the kinetic energy of the projectile without modifying noticeably the propagation direction of CR particle. On the other hand, each elastic collision with an atomic nucleus changes the direction of the trajectory of the projectile [30]. These energy loss processes—the dominance of inelastic versus elastic energy loss processes—depend strongly on the kinetic energy of the projectile. A common feature for each CR particle is that the nuclear energy loss, also called nuclear stopping power (𝑆n), with a maximum at a relatively low energy of the order of 1 keVamu−1; as the kinetic energy increases, the nuclear stopping of the ion decreases. This means that the nuclear stopping power is important only for low ion holding very low velocities. For projectiles with kinetic energies larger than 1 keVamu−1, the energy loss processes are dictated by electronic interaction (electronic stopping power, 𝑆e). The latter leads predominantly to electronic excitation and ionization of the atoms/molecules in the molecular cloud (molecules can also be excited rotationally and vibrationally). It is important to stress that at the beginning of the energy transfer processes from the implant to the surrounding matter at high energies, the ion decelerates mainly by electronic stopping at high energies, therefore almost propagating in a straight path. When the ion has slowed down sufficiently to the keVamu−1 region, the nuclear interaction potential becomes increasingly important. Therefore, to quantify the energy loss of an ion while crossing a molecular cloud, one has to investigate the energy-dependence of the nuclear, (𝑆n(𝐸)), and electronic, (𝑆e(𝐸)), stopping power. The total stopping power 𝑆(𝐸) is simply the sum of both contributions (10). Here, we utilize the SRIM (Stopping and Range of Ions in Matter) code which is based on the modified Bethe-Bloch approximation [31] and is more convenient to use. Note that at kinetic energies close to 1 GeV, the energy loss of each charged particle is also determined by Bremsstrahlung, Cherenkov radiation, and nuclear reactions [30]. Here, we neglect these processes since the differential flux of the CR particles at these energies is at least 4 orders of magnitude lower than, for instance, in the MeV region [2]𝑆(𝐸)=𝑆n(𝐸)+𝑆e(𝐸)=𝑑𝐸𝑑𝑟.(10)

2.3. Energy Loss and Ranges of Charged Particles in Molecular Clouds

It is possible to quantify now the energy loss and the penetration (range) of charged particles inside molecular clouds. To describe a particle penetration depth into a cloud, the projected mean range, 𝑅, is often used in collision dynamics. Compared to the actual path length, 𝐿, of a projectile, the projected range, that is, the straight-line distance from implantation origin to the end point of the trajectory projected in the initial velocity vector, is always smaller. Here, significant angular deflections of the projectile via nuclear interaction happen to be eminent at energies of less than 1 keVamu−1. For example, considering 1 MeV protons, we can make use of (e.g., [19]) a crucial relationship between the mean free path, 𝜆, that is, the average distance between two collisions, of a proton from the CR radiation field, the ionization cross-section of molecular hydrogen 𝜎(1 MeV), and the number density of H2 (11). For 𝜎(1MeV)=31017 cm2 [1] and 𝑛(H2)=103 cm−3, a 1 MeV proton has a free mean path of about 31013 cm H𝜆𝜎(1MeV)𝑛21.(11) We computed the projected ranges and electronic, nuclear, and total stopping powers of protons in molecular hydrogen gas with the SRIM code for a density of 8.9105 gcm−3 [31]. SRIM does not allow a computation of the stopping powers at the low density of interstellar clouds directly. Since the range is inversely proportional to the density of a medium, we can utilize (11) to scale the range 𝑅1 of a proton of an energy 𝐸 in a medium of a density𝜌1to compute the range 𝑅2 of the same proton with an identical energy 𝐸 in a medium of a density𝜌2 via (12). Applying this procedure for molecular clouds with number density of 103 cm−3, we find that 𝑅22.71016𝑅1. Since H2 is the dominating component of these environments, our scaling procedure is well justified. For instance, applying this procedure to 1 MeV protons gives a computed range in molecular clouds with number densities of 103 cm−3 of 2.61017 cm. Thus, one can conclude from such computations that (i) CR protons with energies 𝐸<1MeV cannot penetrate inside molecular clouds with 𝑛(H2)~103 cm−3 deeper than 0.1 pc, (ii) for energies 𝐸1MeV, the CR particles travel essentially straight paths through the cloud since deflections from its trajectory via nuclear interactions are negligible, and (iii) the non attenuated flux of the CR particles is changed due to the energy loss; this in turn slows the particles down and shifts the flux distribution from high values to lower𝑅1𝜌1=𝑅2𝜌2.(12) In a one-dimensional steady-state approximation, neglecting diffusion and internal sources, the proton flux inside a cloud after shielding by a homogeneous layer, 𝐹internal(𝐸) with a depth 𝑑 and a density 𝜌2 can be numerically computed from the stopping power for H2𝑆(𝐸) (Figure 3) and the initial, non attenuated proton flux, 𝐹(𝐸) (Figure 2). Accounting for (10), the energy-dependent energy loss of a particle with a kinetic energy 𝐸, Δ𝐸(𝐸), can be calculated via(13)𝜌Δ𝐸(𝐸)=𝑆(𝐸)𝑑2𝜌1.(13) Considering energy conservation, we have to bear in mind that the external energy flux 𝑄external(𝐸)—in units of eV cm−2 s−1—of the protons outside the molecular cloud must be equal to the energy flux from the decelerated protons inside the cloud, 𝑄internal(𝐸), plus the energy flux given to the target atoms, 𝑄target(𝐸)𝑄external(𝐸)=𝑄internal(𝐸)+𝑄target(𝐸).(14) Here, the external energy flux 𝑄external(𝐸) can be computed as the product of the energy 𝐸 and the non attenuated flux 𝐹(𝐸) (Figure 4) to𝑄external(𝐸)=𝐹(𝐸)𝐸.(15) Likewise, the internal energy flux presents the product of the energy 𝐸 and the internal flux 𝐹internal(𝐸)—the quantity we want to compute—to𝑄internal(𝐸)=𝐹internal(𝐸)𝐸.(16) To calculate the energy flux transferred from particles of an energy 𝐸 given to the target molecules (the H2 gas inside the cloud), we have to account for the following scenario. For particles at a given energy 𝐸, the energy flux contributing from these particles is being reduced (they are actually slowed down from an energy 𝐸 to lower energies) via interaction with the target gas by 𝐹(𝐸)Δ𝐸(𝐸) (c.f. (13) for the definition of Δ𝐸(𝐸)); likewise, particles of an energy 𝐸, which is higher than 𝐸, can be decelerated to the energy 𝐸; this process actually increases the energy flux for a given energy 𝐸 by 𝐹(𝐸)Δ𝐸(𝐸).Thus, a particle with an energy 𝐸 transverses a column density𝑑𝜌2 and exits with an energy 𝐸𝐸=𝐸𝐸Δ𝐸.(17) Combining these relationships with (14)–(16) and inserting these into (13) gives the method to calculate 𝐹internal(𝐸) via (18) from which we can extract then 𝐹internal(𝐸)𝐹(𝐸)𝐸=𝐹internal𝐸(𝐸)𝐸+𝐹(𝐸)Δ𝐸(𝐸)𝐹𝐸Δ𝐸.(18) One can now utilizes these attenuated fluxes 𝐹internal(𝐸) to investigate quantitatively to what extent the CR protons ionize molecular hydrogen via (1). Accounting for the ionization cross-section (Figure 1) and the attenuated CR proton fluxes (Figure 4), the depth-dependent ionization rate, 𝜁(𝑟) (s−1), can be derived via (19). Here, the integration limits range from 𝐸min (the lowest kinetic energy of the ions penetrating to the depth, 𝑟) to the maximum energy of 1010 eV (Figures 1 and 2).𝜍(𝑟)=𝐸max𝐸min𝜎(𝐸)𝐹𝑠(𝐸)𝑑𝐸.(19)

3. Discussion

Table 2 shows the results of our computations for six molecular cloud models; the parameters of these models are compiled in Table 1. The corresponding internal fluxes and depth-dependent ionization rates are depicted in Figure 4. Based on these data, we can derive the following principal conclusions. (1)The energy absorbed by the molecular hydrogen inside the molecular cloud at a distinct depth is very small (8.3×107 eV cm−2s−1 at the most) compared to the non attenuated energy flux of 4.3×109 eV cm−2s−1 [2] and 3.1×109 eV cm−2s−1 [3]. This translates to an energy loss of the CR protons of less than 2% irrespective on the model; therefore, the charged particles transverse the cloud almost in straight trajectories nearly unperturbed. Further, as expected from the computed stopping power (Figure 3), electronic energy loss presents the dominating energy transfer mode and is three orders of magnitude more important than the nuclear energy transfer. Summarized, the energy loss of the charged particles is irrespective on the choice of the non attenuated flux, on the number density of molecular hydrogen (at least in the 103 to 104 cm−3 regime), and on the selection of 𝐸min (10 MeV, 50 MeV, and 100 MeV). (2)Even though in [2], total non attenuated proton flux of 7.4×104 cm−2 s−1 is four orders of magnitude larger than the value of [3] (9.3 cm−2 s−1), the majority of flux distribution in [2] are low-energy protons. These species are magnetically deflected and hence do not penetrate inside the cloud. This yields to a very narrow range of proton fluxes of—depending on the depth-6 protons cm−2s−1 to 21 protons cm−2s−1 [2] and 6 protons cm−2s−1 to 9 protons cm−2s−1 [3]. (3)Our computed ionization rates range between 2.51017 s−1 and 4.11016 s−1. It is interesting comparing these data with ionization rates claimed by observations, that is, 10−18 s−1–10−16 s−1. The present computations predict that—depending on the parameters utilized in models M1–M6—the ionization rate can be in certain areas of the molecular cloud higher than previously postulated. Likewise, our calculations suggest minimal ionization rates of 2.51017 s−1; this value is in strong contrast to estimated data of about 10−18 s−1. Also, the resulting ionization rates decreases slightly toward the center of the molecular cloud as predicted by [13].

It should be stressed here that these results depend essentially on the cut-off procedure of the CR spectra at low energies (10 to 100 MeV); otherwise the contradiction between calculated and observed ionization rates may be even stronger. For example, the CR spectra when protons with energies 1 MeV ≤ 𝐸 ≤ 10 GeV penetrate inside clouds cause ionization rates close to 10−15 s−1 under typical clouds conditions (𝑛~103 cm−3) (Figure 5) or to 10−16 s−1 (𝑛106 cm−3) Figure 6).

4. Astrophysical Implications

It is interesting to compare our calculations with those data derived from astronomical observations. We would like to stress that our computed values of the ionization rate are in good agreement with an (indirectly) observed data point. Brittain et al. [32] measured H3+ absorption lines in the spectrum of the Herbig Be star LkHα 101 suggesting a column density of 𝑁(H3+)=2.21014 cm−2 [32]. A more detailed investigation of the data revealed that the absorption lines must be correlated with a dense cloud with number densities of 𝑛(H2)~104 cm−3, the length along the line of sight of 𝐷~0.4 pc, and an averaged kinetic temperature of 𝑇~30 K. To illustrate a validity of the calculated value of the fractional ionization (Figure 4), we can apply (3)–(5) to observations [32]. Here, utilizing a ratio of the column density of CO and H2 in this object of 𝑁(CO)/𝑁(H2)=1.8104, we derive 𝑛(H3+)=3.31012𝜁 and, hence, a theoretically predicted column density of 𝑁(H3+)=2.51014 cm−2; this compares very well with the observed value of 𝑁(H3+)=2.21014cm2. Moreover, the calculated value of the ionization rate throughout the cloud, 𝜁=6.01017 s−1, while the observed value is of 6.71017 s−1. Here, it is also worth mentioning that the CR flux—assumed similar in the solar neighborhood and at the edge of the non attenuated cloud—is adequately transformed after shielding corrections inside the dense cloud toward LkHα 101. This model has been shown to correctly reproduce the astronomical observations of column densities of H3+. However, one should point out that the calculated fluxes are sensitive to the density distribution in the cloud; in approximation used one assumed a homogeneously mixed cloud; non homogeneously mixed clouds should demand on more refined procedure.

Further studies have the goal to expand this model to investigate how CR particles heavier than protons influence the ionization rate 𝜁. In addition, one will examine the effect of the chemical composition of the molecular cloud. Recall that in the present model it is assumed the cloud consists entirely of molecular hydrogen; therefore, the incorporation of atomic hydrogen, helium, and carbon monoxide is worth investigating. Finally, the sensitivity of the model to dust particles inside molecular clouds—although having very low number densities of about 1012𝑛 grains cm−3—will be probed. These extensions will also yield valuable information on the energy deposition inside grain particles which in turn can help to understand the energy dependent formation of molecules on icy grains deep inside molecular clouds [7]. This would allow an objective, unbiased comparison of the energy deposition inside ice-coated grains by CR versus UV photons and would expand Shen et al.’s [33] preliminary study significantly.

In summary, it is shown that the distribution of the CR particle flux inside a molecular cloud can be calculated if the energy loss of energetic protons—the principal constituent of the cosmic radiation field—in molecular hydrogen is known. The latter can be computed utilizing the SRIM program package and rescaling the computed data to the actual number densities of the cold molecular clouds. Integrating over the energy-dependent ionization cross-section of molecular hydrogen, a simple numerical procedure finally yields the ionization rates of molecular clouds. A phenomenological cutting-off procedure of the CR spectra at low energies 1MeV𝐸100MeV permits to bring together the calculated and observed ionization rates. Based on astronomically observed H3+ absorption lines, calculated values of the ionization rates are well within those derived from astronomical observations. These environment-dependent ionization rates are in turn critical parameters to establish updated physical and chemical models of molecular clouds.

Acknowledgments

Partially this work was supported by the Fulbright Scholar Program (USA). The author thanks Ralf Kaiser for helpful discussions. Useful comments from two anonymous referees helped to improve the presentation of this paper.