Issue 
A&A
Volume 616, August 2018



Article Number  A59  
Number of page(s)  4  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201832704  
Published online  17 August 2018 
Analytic model of comet ionosphere chemistry
Swedish Institute of Space Physics, Uppsala, Sweden
email: erik.vigren@irfu.se
Received:
25
January
2018
Accepted:
13
May
2018
Context. We consider a weakly to moderately active comet and make the following simplifying assumptions: (i) The partial ionization frequencies are constant throughout the considered part of the coma. (ii) All species move radially outward with the same constant speed. (iii) Ionneutral reactions affect the chemical composition of the ions, but ion removal through dissociative recombination with free electrons is negligible.
Aims. We aim to derive an analytical model for the radial variation of the abundances of various cometary ions.
Methods. We present two methods for retrieving the ion composition as a function of r. The first method, which has previously been used frequently, solves a series of coupled differential equations. The new method introduced here is based on probabilistic arguments and is analytical in nature.
Results. For a pure H_{2}O coma, the resulting closedform expressions yield results that are identical to the standard method, but are computationally much less expensive.
Conclusions. In addition to the computational simplicity, the analytical model provides insight into how the various abundances depend on parameters such as comet production rate, outflow speed, and reaction rate coefficients. It can also be used to investigate limiting cases. It cannot easily be extended to account for a radially varying flow speed or dissociative recombination in the way a code based on numerical integrations can.
Key words: comets: general / molecular processes
© ESO 2018
1. Introduction and nomenclature
When a comet comes sufficiently close to the Sun, it becomes active, releasing volatiles from nearsurface layers. Dust grains leave the comet with the gas flow. A coma and a dust tail can be seen by the reflection of sunlight on this dust. The gas in the coma can be ionized, either by particle impact or solar extreme ultraviolet radiation. The ions are eventually picked up by the solar wind, giving rise to an ion tail that is sometimes visible as a result of photon emission in the optical region of the electromagnetic spectrum. In this work we focus on ions within ~100 km from a cometary nucleus. We consider a purely waterdominated coma and assume that all considered species move radially outward with the same velocity u. As the photolysis frequencies are so low for the relevant heliocentric distances (several astronomical units), we may neglect loss of H_{2}O over the cometocentric distances we consider. We therefore assume that the H_{2}O number density, n_{N}, follows a simplified version of the Haser (1957) model:(1)
where Q is the outgassing rate (the number of H_{2}O molecules released from the cometary surface per second), u is the expansion velocity (typically on the order of 0.5–1 km s^{−1}), and r is the cometocentric distance. Assuming a constant irradiation level (optically thin medium, which is reasonable at least up to Q values of a few times 10^{27} s^{−1}, see, e.g., Vigren et al. 2015), the probability that a given H_{2}O molecule will be ionized on its journey from the surface (at r_{c}) to r is given by the ionization frequency, f, times the travel time. When the ions move with the same velocity as the neutrals (which we assume throughout this work) and when they are not lost through dissociative recombination, the electron number density, n _{e} (by quasineutrality equal to the total ion number density) is given by(2)
While Eq. (2) is a handy analytical expression for predicting the total ion number density at any cometocentric distance given a local measurement of the neutral density and an estimate of the expansion velocity, it holds no information about the ion composition. The photo or electronimpact ionization of H_{2}O leads with different branching fractions to H_{2}O^{+}, OH^{+}, O^{+}, or H^{+ }(e.g., Huebner et al. 1992; Vigren et al. 2015; Itikawa & Mason 2005); we consider here branching fractions of 0.703, 0.197, 0.011, and 0.088, respectively. Reactions of these ions with H_{2}O will affect the relative abundances of the different ion species. The five reactions considered in this work (note that we here only write out the ionic product) are H^{+} + H_{2}O → H_{2}O^{+}, O^{+} + H_{2}O → H_{2}O^{+}, OH^{+} + H_{2}O → H_{2}O^{+}, OH^{+} + H_{2}O → H_{3}O^{+}, and H_{2}O^{+} + H_{2}O → H_{3}O^{+}. Rate coefficients, which are on the order of 10^{−9} cm^{3}s^{−1}, may be found at www.udfa.net (McElroy et al. 2013). We use values at 300 K.
We start in Sect. 2 by briefly describing a standard approach of cometary ionospheric modelling as used, for example, by Giguere & Huebner (1978), Vigren & Galand (2013) and Heritier et al. (2017). We note that while details differ in these works, they all work in “r bins” from the surface and outward, and numerically solve a system of coupled differential equations (continuity equations) in each bin. In Sect. 3 we present a different approach that is based on probability arguments. The advantages and disadvantages of the new method compared with the standard method are discussed in Sect. A, together with some concluding remarks. Table 1 summarizes essential parameters of the work.
Parameters relevant for the study.
2. Standard method
The method is used in Vigren & Galand (2013), for example, and is described here only briefly. For a radially expanding coma (and constant expansion velocity u), the continuity equation for ion species i at r is given by(3)
where P_{i} is the production term and L_{i} is the chemical loss rate of i. In the case of a pure water coma and where loss through dissociative recombination is negligible, we have(4)
The numerical solution of the system consists of dividing the coma into spherical shells where densities are assumed to be constant. In the innermost shell, there is no ion inflow from below. A set of initial ion number densities (close to zero) is set in the innermost shell, and the equation system is run with an appropriate (low) time step until convergence (according to a predefined criterion) is reached. Thereafter, the second shell is considered, now accounting for the ion inflow from below and using as initial densities the “final” densities in the shell below. The procedure continues through all shells.
3. New method
The arrival of an H^{+}, O^{+}, or OH^{+} ion at r requires that an H_{2}O molecule has undergone dissociative ionization at some point r′ and that the resulting ion avoid reactions with H_{2}O molecules on its way from r′ to r. The probability expression becomes (i = 0, 1, 2)(6)
with parameters as described in Table 1, but we recall that(7)
For guidance, the factor _{fi}dr′/u in Eq. (6) corresponds to the probability that a neutral molecule is ionized into ion species i during the time interval dt = dr′/u and the exponential factor represents the probability, P(r,r′), for an ion “born” at r′ to reach r without undergoing reactions during the journey. We note that dP/P = −k_{i}n_{N}dr″/u = −η_{i}dr″/r″^{2}. Using the chain rule, we arrive at the relation d(ln P)/dr″ = − η_{i}/r″^{2}. Integration then gives ln(P(r,r′)) – ln(P(r′, r′)) = η_{i}(1/r–1/r′), and since P(r′,r′) = 1 and ln 1 = 0, we obtain the desired relation P(r, r′) = exp(η_{i}(1/r–1/r′)). The integral in Eq. (6) is readily evaluated as (see Appendix A)(8)
where we have defined the exponential integral as(9)
The arrival of an H_{2}O^{+} ion at r requires either (i) that an H_{2}O molecule undergo nondissociative ionization at some point r′ and that the resulting H_{2}O^{+} ion avoid reactions with H_{2}O molecules on its way from r′ to r, or (ii) that an H^{+}, O^{+} or OH^{+ }ion reaches r′ where it reacts with H_{2}O to produce H_{2}O^{+}, which in turn avoids reactions with H_{2}O molecules on its way from r′ to r. The resulting probability expression becomes(10)
In compact notation, we may express this as(11)
where the indexing shows the pathway to H_{2}O^{+}, for instance, 2 → 3 means that it proceeds via OH^{+}. We have from Eq. (8) that(12)
For i = 0, 1, 2 one can from Eqs. (8) and (10) arrive at(13)
where A_{i}, B_{i}, C_{i}, and D_{i} represent four integrals similar to Eq. (10) thatcanbeevaluatedbystandardmethods(e.g.,Geller & Ng1969) to yield (for i = 0, 1, 2)(14) (15) (16) (17)
Equations (10)–(18) allow calculation of K_{3}(r). We then have for i = 0, 1, 2, 3(19)
and for i = 4, that is, H_{3}O^{+},(20)
where n_{e}(r) is given by Eq. (2).
4. Discussion and conclusions
We have tested the two methods against each other and obtained good agreement (see Fig. 1 for a case example) where minor deviations (seen only near the surface) are due to the errors associated with the use of discrete time and altitude steps in the standard method. The new method is extremely fast. It does not include time steps, and it can generate densities at a given r without first having to calculate densities closer to the nucleus. These are the main advantages of the new method. A further advantage with closedform expressions is that it allows us to see how different terms scale and how the densities evolve in different limits.
Fig. 1.
Ion number densities calculated from the standard model (colored lines) and from the analytical model (overplotted dashed lines). For these simulation Q = 2 × 10^{26} s^{−1} and u = 600 m s^{−1}. Discrepancies near the surface are attributed to the discrete altitude levels and time steps considered in the standard method. 
The new method also has disadvantages. It is set up based on simplifying assumptions such as a constant expansion velocity, constant ionization frequencies, constant reaction rate coefficients, and the assumption that dissociative recombination (leading to removal of ions) is negligible. The standard method can accomodate varying parameters as well as dissociative recombination. Still, the analytical approach can be used as a limiting test for numerical codes. Another disadvantage of the new method is that it requires much more effort in order to implement an extended set of neutrals and associated ions and ionneutral reactions. For example, if NH_{3} were to be added, would overtake the role of H_{3}O^{+} as terminal ion , and the derivation of an analytical expression for the density of H_{3}O^{+} would have to include pathways involving two binary reactions. With the addition of more species, the number of ion production pathways increases, and some of the pathways can include many reaction steps. While in principle analytical solutions are there, they may be tedious to find (or even to implement).
A final limitation of the new method is that when the activity is high, the exponential factors can grow larger than what many solvers can handle. Taylor expansion on multiplications of exponential and E functions may offer a way around such obstacles. However, for activity levels that show this problem, some of the other simplifying assumptions are already highly questionable, and therefore we do not expand on the issue further here. While we here considered a coma dominated by H_{2}O, it is noted that consideration of a H_{2}O/CO_{2} coma, for instance, might to first approximation be approached using expressions in similar forms as those presented here, and we therefore postpone derivations of closedform expressions for reaction pathways involving more than one binary reaction to the future. , the dominant ion produced in the photolysis of CO_{2}, is nonreactive with CO_{2}, but reactive with H_{2}O to mainly form H_{2}O^{+}. Including this chemistry would only result in additional terms in the formalism we have presented. Extended pathways via minor primary ions are not expected to contribute greatly to the buildup of H_{2}O^{+} (and H_{3}O^{+}), although for detailed comparisons with results from standard models, closedform expressions for pathways with at least two binary reactions are desirable. The influence of the electric and magnetic field environment is a greater concern; this is not explicitly accounted for in either the standard or the new method we presented here. The ionneutral decoupling distance as fomulated in Eq. 3 in Mandt et al. (2016; see also their Fig. 5), for example, could give a rough indication of where the model is applicable for a given activity level, but the actual situation is likely very complex. The closedform expressions presented here offer a fast way to compare modelderived relative ion abundances with observations from ROSINA/DFMS measurements (e.g., Fuselier et al. 2016; Heritier et al. 2017) in the coma of comet 67P/ChuryumovGerasimenko, the target comet of the Rosetta mission. Marked variations in relative ion abundances seen in observations and not captured by the model could give important clues about spatial and/or temporal variations in the field environment.
Acknowledgments
We thank the Swedish National Space Board for financial support (Dnr 166/14) and Anders Eriksson for valuable discussions. We also acknowledge the International Space Science Institute; this work was done in the framework of the ISSI International Team on the plasma environment of comet 67P after Rosetta, number 402.
References
 Fuselier, S. A., Altwegg, K., Balsiger, H., et al. 2016, MNRAS, 462, S67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Geller, M., & Ng, E. W. 1969, J. Res. Nat. Bur. Stand. B. Math. Sci., 73, 191 [CrossRef] [Google Scholar]
 Giguere, P. T., & Huebner, W. F. 1978, ApJ, 223, 638 [NASA ADS] [CrossRef] [Google Scholar]
 Haser, L. 1957, In Bulletin de la Classe de Science, Académie Royale de Belgique, 43, 740 [Google Scholar]
 Heritier, K. L., Altwegg, K., Balsiger, H., et al. 2017, MNRAS, 469, S427 [Google Scholar]
 Huebner, W. F., Keady, J. J., & Lyon, S. P. 1992, Ap&SS, 195, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Itikawa, Y., & Mason, N. 2005, J. Phys. Chem. Ref. Data, 34, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Mandt, K. E., Eriksson, A. I., Edberg, N. J. T., et al. 2016, MNRAS, 462, S9 [NASA ADS] [CrossRef] [Google Scholar]
 McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vigren, E., & Galand, M. 2013, ApJ, 772, 33 [Google Scholar]
 Vigren, E., Galand, M., Eriksson, A. I., et al. 2015, ApJ, 812, 54 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A
Evaluation of Eq. (6)
We note that Eq. (6) can be rewritten in the form(A.1)
where A = (f_{i}/u)exp(η_{i}/r). We can integrate by parts to obtain(A.2)
For the latter integral, we make the variable substitution η_{i}/r′ = t so that dr′= η_{i}dt/t^{2} and(A.3)
where we have taken care of a sign by switching the integration limits. We realize that(A.4)
From Eqs. (9) and (A.4), we see that(A.5)
Inserting this into Eq. (A.3) and applying some trivial algebra yields the desired Eq. (8).
All Tables
All Figures
Fig. 1.
Ion number densities calculated from the standard model (colored lines) and from the analytical model (overplotted dashed lines). For these simulation Q = 2 × 10^{26} s^{−1} and u = 600 m s^{−1}. Discrepancies near the surface are attributed to the discrete altitude levels and time steps considered in the standard method. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.