A 1D Model of Radial Ion Motion Interrupted by Ion–Neutral Interactionsin a Cometary Coma

In an H₂O-dominated coma, critical ion–neutral reactions to consider are the electron transfer reaction,

$$\begin{eqnarray}&&{{\rm{H}}}_{2}{{\rm{O}}}^{+}+{{\rm{H}}}_{2}{\rm{O}}\to {{\rm{H}}}_{2}{\rm{O}}+{{\rm{H}}}_{2}{{\rm{O}}}^{+},\end{eqnarray} \tag{ 9a }$$

and the proton-transfer reaction,

$$\begin{eqnarray}&&{{\rm{H}}}_{2}{{\rm{O}}}^{+}+{{\rm{H}}}_{2}{\rm{O}}\to \mathrm{OH}+{{\rm{H}}}_{3}{{\rm{O}}}^{+}.\end{eqnarray} \tag{ 9b }$$

In both of these interactions it is assumed that if the ion is the faster of the reactants, the ion will be the slower of the products. These reactions are discussed in some depth in Fleshman et al. (2012), who refer specifically to experiments by Lishawa et al. (1990), and present usable expressions of the cross sections as a function of relative kinetic energy (defined in Table 1), E_rel, namely,

$$\begin{eqnarray}&&{\sigma }_{{et}}({E}_{\mathrm{rel}})=38{E}_{\mathrm{rel}}^{-0.5}\end{eqnarray} \tag{ 10a }$$

for Equation (9(a)) and

$$\begin{eqnarray}&&{\sigma }_{\mathrm{pt}}({E}_{\mathrm{rel}})=38{E}_{\mathrm{rel}}^{-0.88}-0.39\times \exp \left\{-\tfrac{1}{2}{\left(\tfrac{{E}_{\mathrm{rel}}-57}{12}\right)}^{2}\right\}\end{eqnarray} \tag{ 10b }$$

for Equation (9(b)), where E_rel is to be inserted in eV to get cross sections in 10⁻¹⁶ cm². Note that Lishawa et al. (1990) also examined the relative importance of H atom transfer for Equation (9(b)), finding it to be at least an order of magnitude less efficient than proton transfer. This channel is neglected in our study, so when accounting for Equations (9(a)) and (9(b)), we assume that the product ion has a velocity of u. This treatment, giving the product ion the velocity of the neutrals, is justified by the dominance of the spectator stripping type of mechanism as verified for the collision energies investigated by Lishawa et al. (1990) and also by Vančura & Herman (1991). It is noted that the contribution from a complex decomposition mechanism, proceeding via a long-lived intermediate ion, was found to increase with decreasing collision energy, implying that for collision energies below ∼1 eV our treatment may not be a particularly good approximation. Limited by the 1D nature of our model, we will stick with our assumption and merely note that it is one in which the overall momentum transfer from ions to neutrals may be somewhat overestimated. Similar to Fleshman et al. (2012), we will explore the effect of using Equation (10(b)) for the proton-transfer cross section at relative kinetic energies >1.5 eV but a different function at lower kinetic energies, namely,

$$\begin{eqnarray}&&{\sigma }_{\mathrm{pt}}({E}_{\mathrm{rel}}\lt 1.5\,\mathrm{eV})=30{E}_{\mathrm{rel}}^{-0.50},\end{eqnarray} \tag{ 10c }$$

again with E_rel in eV giving cross sections in the units of 10⁻¹⁶ cm². Model runs in which the latter relation is considered will be marked by asterisks.

It is noted that the cross section for Equation (9(a)) has been expressed differently in other studies referring also to the work of Lishawa et al. (1990). For example, Cassidy & Johnson (2010) mention the expression 10^−14.8 E_rel^−0.18 with E_rel in eV giving cross section in cm², while Sakai et al. (2016) set the cross section equal to 10⁻¹⁵ cm², independent of energy. For E_rel of 0.1, 1, and 3 eV these different relations give cross-section values that differ by up to a factor of ∼12, ∼4, and 1.2, respectively. The pronounced difference at low energies is not so critical, as interactions there anyway are dominated by proton transfer at least when utilizing the cross-section relation given by Equation (10(b)) regardless of energy. In Section 5.5 we will return to the proton-transfer reaction 9(b), but with an emphasis on very low collision energies. We may reveal already here that the use of Equation (10(b)) also at very low collision energies gives a 300 K reaction rate coefficient, higher than the one experimentally derived by Huntress & Pinizzotto (1973).

The electron transfer process

$$\begin{eqnarray}&&{{\rm{H}}}_{3}{{\rm{O}}}^{+}+{{\rm{H}}}_{2}{\rm{O}}\to {{\rm{H}}}_{3}{{\rm{O}}}^{* }+{{\rm{H}}}_{2}{{\rm{O}}}^{+}\end{eqnarray} \tag{ 11 }$$

is endoergic by ∼6 eV, and we are leaving out this reaction in our model. For the estimate of the endoergicity we assumed that the H₃O* complex rapidly dissociates into H₂O+H, and so we compared the heat of formation of H₃O⁺ (∼6.2 eV) with those of H₂O⁺ (∼10.1 eV) and H (∼2.3 eV).

Gombosi (2015, pp. 169–88), referring to Cravens & Körosmezey (1986), uses an H₃O⁺/H₂O "charge transfer" rate coefficient of k_in = 1.1 × 10⁻⁹ cm³ s⁻¹. Due to the endoergicity of the electron transfer, this "charge transfer" should refer to the proton transfer

$$\begin{eqnarray}&&{{\rm{H}}}_{3}{{\rm{O}}}^{+}+{{\rm{H}}}_{2}{\rm{O}}\to {{\rm{H}}}_{2}{\rm{O}}+{{\rm{H}}}_{3}{{\rm{O}}}^{+},\end{eqnarray} \tag{ 12 }$$

wherein a fast reactant ion is replaced by a slow product ion. While involving no direct change in the chemistry, this is very likely a reaction of great importance for the ion–neutral coupling in cometary comae. We have used the above k_in value to crudely estimate a cross-section relation for reaction (12):

$$\begin{eqnarray}&&{\sigma }_{\mathrm{pt}}^{* }({E}_{\mathrm{rel}})=24\,{E}_{\mathrm{rel}}^{-0.50},\end{eqnarray} \tag{ 13 }$$

where, as before, insertion of E_rel in eV gives a cross section in 10⁻¹⁶ cm². The derivation of Equation (13) involved making the product of the relative velocity (proportional to the square root of E_rel) and the cross section match the value of k_in.

The cross section used for reaction (12) may be questioned, as k_in was treated as temperature independent. Should k_in have a negative temperature dependence, as may be expected from the polar nature of the water molecule, the cross section could have a steeper (negative) energy dependence than suggested by Equation (13) at low energies. Moreover, through the selected ion flow tube (SIFT) technique, Smith et al. (1980) investigated reaction (12), but with either the reagent ion replaced by D₃O⁺ or the reagent neutral replaced by D₂O. Rate coefficients for the interactions were reported as ∼2 × 10⁻⁹ cm⁻³ s⁻¹ at room temperature, roughly a factor of 2 higher than k_in and comparable with the room temperature rate coefficient of Equation (9(b)) (see Huntress & Pinizzotto 1973). In Section 5.1.2 we test, for these reasons, the effect on model output of utilizing Equation 10(b) not only for reaction 9(b) but also for the proton transfer from H₃O⁺ to H₂O. However, an interesting finding in the SIFT study of both the D₃O⁺ + H₂O and D₃O⁺ + H₂O interactions was that the dominant channel (65%–70%) involved D/H exchange, implying a mechanism in which the ion and neutral come together as an intermediate adduct ion existing long enough for "total scrambling" of the D and H atoms to occur before separating into a product ion and a product neutral. In the comet environment it may thus be the case that the interactions between a fast H₃O⁺ ion and a slow H₂O molecule do not necessarily result in a slow ion and fast neutral, similar to the discussion in Section 3.1 for the proton transfer from H₂O⁺ to H₂O at low collision energies.

The 1D nature of the TRAIN model makes a proper treatment of elastic collisions, or for that matter collisions involving transfer of energy between translational and internal degress of freedom, difficult to implement. In such ion–neutral interactions, involving no charge transfer, the ion may be significantly deflected. In head-on-collisions the momentum of the impinging ion can be transferred to the neutral, leaving the ion with a velocity close to u. However, experiments on elastic scattering often reveal a propensity for forward scattering of the ion. The total scattering cross section may be roughly similar to those for charge transfer processes (see, e.g., Cassidy & Johnson 2010, and references therein). Rather than attempting to account for elastic scattering in a detailed manner, we will consider two limiting cases. In Model A we will completely ignore elastic collisions, setting the cross section to 0. This may be regarded as a forward scattering treatment, running with the assumption that even if elastic scattering occurs, they will tend not to significantly change the trajectory of the ions or significantly reduce their velocity. In Model B, on the other hand, we assign the cross-section relation of Equation (13) and consider that whatever ion is undergoing the interaction with H₂O, it will obtain velocity u.

For transparency, we summarize in Table 2 our considered set of interactions and associated energy-dependent cross sections for models A, A*, B, and B*. Suggestions for improvements of the cross-section set are warmly welcome. As mentioned in Section 3.3, we test also one additional model (not listed in Table 2) with a slightly different cross-section set, particularly with the use of the same cross-section relation 10(b) for the proton-transfer relations 9(b) and (12). This model, giving a more extended range of a collisionally coupled plasma, is further discussed in Sections 5.1.2 and 5.5. It is also worthwhile to point out already at this stage that models A and B are associated with proton-transfer cross sections at low collision energies (sub-eV) that, when folded over Maxwellian energy distributions, yield rate coefficients in excess of experimentally determined values at room temperature (see further details in Section 5.5).

Under the assumptions of a constant ionization frequency, a constant outflow speed u (similar for neutrals and ions as E(r) = 0), and a negligible plasma neutralization through dissociative recombination, an analytical expression for the total ion number density can be derived. From charge neutrality (n_e = n_I) we will refer to this analytical expression as n_e,th (e.g., Vigren et al. 2016). It is given by

$$\begin{eqnarray}&&{n}_{e,\mathrm{th}}=\tfrac{Q\nu }{4\pi {u}^{2}{r}^{2}}(r-{r}_{c}),\end{eqnarray} \tag{ 14 }$$

where ν is the ionization frequency.

For r ≫ r_c this is well approximated by

$$\begin{eqnarray}&&{n}_{e,\mathrm{th}* }=\tfrac{Q\nu }{4\pi {u}^{2}r},\end{eqnarray} \tag{ 15 }$$

showing that a 1/r decay in the electron number density is anticipated in a field-free case and for r ≫ r_c.

The n_I,i (n_I at r_i) obtained from the TRAIN model is

$$\begin{eqnarray}&&{n}_{I,i}=\tfrac{Q\nu {dr}}{4\pi {{ur}}_{i}^{2}}\left(\sum _{j=1}^{i}\,\displaystyle \frac{{c}_{i,j}}{{v}_{i,j}}+\displaystyle \frac{{d}_{i,j}}{{v}_{i,j}^{* }}\right)\end{eqnarray} \tag{ 16 }$$

and is in excellent agreement with n_e,th when E(r) = 0 and regardless of activity level and the step size dr.

In the limit of low activity and for an electric field E inversely proportional to r, we have also verified that the resulting n_I profile from the TRAIN model reproduces the analytical expression in Appendix B of Vigren et al. (2015; see their Equation (25)).

The analytical expression for the H₃O⁺/H₂O⁺ number density ratio for a pure H₂O coma (and under similar assumptions to those listed in Section 3.1) is independent of the ionization frequency (provided that it is constant and nonzero) and given by the rather involved expression (Vigren et al. 2016)

$$\begin{eqnarray}&&{R}_{\mathrm{dens},{th}}=\displaystyle \frac{r-{r}_{c}}{r-{r}_{c}\exp \left(\eta r-\tfrac{\eta {r}^{2}}{{r}_{c}}\right)+\eta {r}^{2}\exp (\eta r)\left\{{E}_{n}\left(\tfrac{\eta {r}^{2}}{{r}_{c}}\right)-{E}_{n}(\eta r)\right\}}-1,\end{eqnarray} \tag{ 17 }$$

where E_n is the exponential integral function

$$\begin{eqnarray}&&{E}_{n}(x)={\int }_{x}^{\infty }\tfrac{{e}^{-t}}{t}{dt},\end{eqnarray} \tag{ 18a }$$

and the function η is given by

$$\begin{eqnarray}&&\eta (r)=\tfrac{{n}_{n}(r){k}_{\mathrm{pt}}}{u},\end{eqnarray} \tag{ 18b }$$

with k_pt being the rate coefficient for the proton-transfer reaction H₂O⁺ + H₂O H₃O⁺ + OH. In order to compare with the TRAIN model, we set k_pt/u equal to the proton-transfer cross section at the relative kinetic energy corresponding to v_rel = u_T (this does not adequately describe the proper relation between the rate coefficient and the cross section; we have, for example, not yet considered that a nonradial component of the ion motion yields extended pathways in the considered regions, a matter that we return to in Section 5.3.2).

In the TRAIN model the H₃O⁺/H₂O⁺ number density ratio R_dens,i is computed as

$$\begin{eqnarray}&&{R}_{\mathrm{dens},i}=\sum _{j=1}^{i}\,\tfrac{{d}_{i,j}}{{v}_{i,j}^{* }}/\sum _{j=1}^{i}\,\tfrac{{c}_{i,j}}{{v}_{i,j}},\end{eqnarray} \tag{ 19a }$$

while the H₃O⁺/H₂O⁺ flux ratio R_flux,i is computed as

$$\begin{eqnarray}&&{R}_{\mathrm{flux},i}=\sum _{j=1}^{i}\,{d}_{i,j}/\sum _{j=1}^{i}\,{c}_{i,j}.\end{eqnarray} \tag{ 19b }$$

For E = 0 simulations R_dens = R_flux and the level of agreement with the analytical expression (Equations (17)–(18) improves with decreasing values of dr. When using a step size of 100 m, an agreement within 5% with the analytical expression arises outside a cometocentric distance r_agree, depending on the activity level. The level of agreement remains good and approaches excellence for r > r_agree. For Q values of 10²⁵, 10²⁶, 10²⁷, and 10²⁸ s⁻¹ the values of r_agree are ∼2, 7, 23, and 70 km, respectively. If using a step size of 50 m, the last two values decrease to 16 and 48 km. Below r_agree the TRAIN model gives lower ratios than the analytical expression. This is primarily caused by the treatment of regional absorption, where slow ion species are fed to the subsequent region as either H₂O⁺ or H₃O⁺ weighted by the cross sections for electron transfer and proton transfer. It is clear that if the regional absorption probability is high, one may expect multiple interactions within the same region, enhancing the probability for H₃O⁺ delivery to the subsequent region.

In the presence of an electric field the regional absorption probability for H₂O⁺ ions is reduced, and so the neglect of multiple interactions within the same region should be less of an issue than in the field-free case. Sensitivity tests are recommended for a location of interest, assuring, for example, that computed R_dens,i or R_flux,i values are only affected in a negligible manner by utilizing a step size reduced by, say, a factor of 4.

The mean ion drift velocity at r_i in the TRAIN model is computed as

$$\begin{eqnarray}&&\langle {u}_{i}\rangle =\displaystyle \frac{{\sum }_{j=1}^{i}\,{c}_{i,j}\,{v}_{i,j}+{d}_{i,j}\,{v}_{i,j}^{* }}{{\sum }_{j=1}^{i}\,{c}_{i,j}+{d}_{i,j}},\end{eqnarray} \tag{ 20 }$$

which in the case of field-free simulations (in which case ${v}_{i,j}={{v}_{i,j}}^{* }=u$) obviously reproduces the expected value of u.

Neglecting thermal components, ion velocity and ion energy distribution functions can also be obtained from the c_i,j and v_i,j arrays (H₂O⁺) and from the d_i,j and ${{v}_{i},j}^{* }$ arrays (H₃O⁺). To realize how, an example is appropriate, so let us consider H₂O⁺ ions at r₂₀₀ (i = 200). In the case of a strictly positive electric field strength throughout, the values of v_200,1, v_200,2, ..., v_200,199 will be decreasing and pairwise distinct. These values can be converted to energies _200,j in eV (e.g., _200,73 = m₁₈v_200,732/(2q), with m₁₈ the mass of H₂O⁺ and q the elementary charge). The values of c_200,1, c_200,2, ..., c_200,199 are proportional to the fluxes of H₂O⁺ ions with corresponding energies _200,1, _200,2, ..., _200,199. By dividing the c_200,1, c_200,2, ..., c_200,199 values by the width of the energy bins, we get numbers proportional to the differential flux of H₂O⁺ ions. By interpolating to a common energy grid, the water group ion flux is obtained by summing the contributions from H₂O⁺ and H₃O⁺, and the IEDF is obtained following normalization. Special care needs to be taken in cases in which regions of zero electric field prevail, and if the electric field is zero throughout, the IEDF will consist of two lines, the relative height of which reflects R_dens. Some examples of calculated IEDFs are shown in Section 5.4.

5.1.1. Results from A, A*, B, and B* Models

Here we present results from A, A*, B, and B* model runs (see Table 2 for a description on how the cross-section sets of the models differ) with Q = 2 × 10²⁷ s⁻¹ (Figure 1) and 2 × 10²⁸ s⁻¹ (Figure 2) and with constant electric fields of magnitudes 0 (black lines), 0.03 (red lines), and 0.1 mV m⁻¹ (blue lines). Other parameters utilized in these specific simulations are listed in Table 3. Naturally, for a given nonzero field strength and model type the higher Q gives higher n_I, lower $\langle {u}_{I}\rangle$, and higher R_dens and R_flux. Likewise, for a given Q and model type the higher electric field magnitude gives lower n_I, higher $\langle {u}_{I}\rangle$, and lower R_dens and R_flux. Also as expected, for given Q and E, models with higher proton-transfer cross section at low relative kinetic energies give higher R_dens and R_flux, and models including "head-on" elastic collisions give lower $\langle {u}_{I}\rangle$ and therefore higher n_I.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  Summary of Parameters Utilized in the Simulations Associated with the Results Displayed in Figures 1–4 of the Paper

Par.	Description	Values Applied for Runs Associated with Figures
u	Neutral outflow velocity	1 km s−1 (Figures 1–4)
uT	Minimal considered relative velocity	0.5 km s−1 (Figures 1–4)
Q	Gas production rate	2 × 1027 s−1 (Figures 1 and 3(a)), 2 × 1028 s−1 (Figures 2, 3(b), and 4)
ν	Ionization frequency	4.5 × 10−7 s−1 (Figures 1–4), affects only nI.
E	Electric field	Constant (Figures 1–3), 1/r dependent for r > ren (Figure 4)
ren	Limit below which E is set to zero	100 km (Figure 4), NaN in other figures
FXH+	Fraction of ions of type XH+ at ren	0.5 (Figure 4), NaN in other figures
kDR	Rate coeff. for dissociative recombination	1 × 10−6 cm3 s−1 (Figure 4), NaN in other figures

Note. NaN means that the parameter did not need to be considered in some simulations.

Download table as:  ASCIITypeset image

Note that even for the simulations with Q = 2 × 10²⁸ s⁻¹ all calculated mean ion flow speeds at ∼200 km are markedly higher than the neutral flow speed u, set to 1 km s⁻¹. The thick dashed lines in Figure 2 are results from model B runs with Q = 2 × 10²⁹ s⁻¹. As seen in panel (b), at this high activity level the mean ion flow remains close to u for the cometocentric distances considered. It is noted that we utilized C_i,1 = 1 also in this high-activity case. Attenuation of the impinging EUV light leads in reality to decreasing ionization frequencies toward the surface. However, if we were to account for this, we would get even lower mean ion speed predictions, due to reduced fluxes of ions that have experienced acceleration over longer distances.

5.1.2. Results from a Model Assuming Efficient Proton Transfer from H₃O⁺ at Low Energies

As discussed in Section 3.2, the cross section for the proton transfer from H₃O⁺ to H₂O may have a steeper energy dependence than in relation (13). We test here the effect of using for this process a cross section given by Equation 10(b). As a reminder, Equation 10(b) is an extrapolation by Fleshman et al. (2012) of the experimental data by Lishawa et al. (1990) for the proton transfer from H₂O⁺ to H₂O, and now we use it also for the proton transfer from H₃O⁺ to H₂O, with the exception that for this change we use the cross sections of model A (neglecting elastic collisions and using Equation 10(b) for reaction 9(b) regardless of energy). Figure 3 shows mean ion flow speeds calculated from model runs with Q = 2 × 10²⁷ s⁻¹ (panel (a)) and Q = 2 × 10²⁸ s⁻¹ (panel (b)) with considerations of constant electric fields of 0.03 and 0.1 mV m⁻¹. It is seen that at least for the higher activity considered, the modification of the cross-section set has a significant influence on the coupling, i.e., on the ability of the neutrals to keep the mean ion flow speed close to that of the neutrals. However, as hinted at in Section 3.1 and further discussed in Section 5.5, the A (and B) model is most likely overestimating the H₂O⁺-to-H₂O proton-transfer cross section at very low collision energies, so that when applying the same cross-section relation for the H₃O⁺-to-H₂O proton transfer, one would expect an overestimated coupling efficiency. There is a need for investigations, experimental or theoretical, of ion–neutral reactions at low relative kinetic energies.

Zoom In Zoom Out Reset image size

A modification of the TRAIN model is to consider a limit where newly produced electrons switch from rapidly being thermalized to remaining hot. We set the border at r_en, the electron–neutral decoupling distance. This limit, certainly more gradual in reality, is anticipated well within the zone where ions and neutrals remain coupled (Mandt et al. 2016), and so we may make some simplistic assumptions. We will assume that ions from below enter the shell starting at r_en with velocity u. A fraction f_XH+ of these are of the form XH⁺ (with X being some molecule, e.g., NH₃, with higher proton affinity than H₂O), and the complementary fraction is H₃O⁺ ions. The reason for considering such molecules is because even if they prevail in low abundances, at the level of a percent or so, they can greatly influence the ion chemistry, as revealed by both modeling (e.g., Vigren & Galand 2013) and observations (Beth et al. 2017). For r > r_en we assume that there is no additional production of XH⁺ ions, and that these species will not be lost chemically but will undergo acceleration along an ambipolar electric field. We make an initial assumption of a 1/r decay in the electron number density, which, combined with the constant electron temperature, sets up an ambipolar electric field of magnitude k_BT_e/qr accelerating ions radially outward (note that we neglect magnetic forces). In principle, we could also assume a more general 1/r^γ decay of the electron number density, in which case the magnitude of the utilized ambipolar electric field (k_BT_e/qr) would need to be multiplied by γ. The outward direction of the ambipolar field may certainly be questioned in the near vicinity of r_en, where we envision a sharp increase in T_e. In the following, for the sake of simplicity we will assume that H₂O molecules interrupt the acceleration of XH⁺ ions equally efficient as they do for H₃O⁺ ions. This assumption affects also the formalism below. We do wish to stress, however, that experimental and/or theoretical studies into the mobility of XH⁺ ions in H₂O gas when subjected to low electric fields are desirable.

In order to obtain estimates of key parameters, we may still use the TRAIN model but set r₁ = r_en instead of r₁ = r_C (an action that automatically adjusts, e.g., the N_i vector). In addition, we need to construct a new matrix G_i,j to account for ions originating from below r_en, i.e., from what we refer to as the cold zone. We compute G_i,j from

$$\begin{eqnarray}&&{G}_{i,j}={G}_{i,j-1}+{G}_{j-1,j-1}{\alpha }_{i-1,j-1}^{* }\displaystyle \prod _{k=j-1}^{i-2}{\beta }_{k,j-1}^{* },\end{eqnarray} \tag{ 21 }$$

with G_i,1 = 0 for all i > 1, with the product sign replaced by 1 in cases when i = j, and with G_1,1 given by

$$\begin{eqnarray}&&{G}_{\mathrm{1,1}}=\left[\left\{{\left(1+\tfrac{{{Qk}}_{\mathrm{DR}}\nu }{\pi {u}^{3}}\right)}^{1/2}-1\right\}\tfrac{2\pi {u}^{3}}{{{Qk}}_{\mathrm{DR}}\nu }\right]\tfrac{{r}_{{en}}-{r}_{c}}{{dr}}.\end{eqnarray} \tag{ 22 }$$

The part within brackets is the ratio of Equation (12) of Gombosi (2015, pp. 169–88) and Equation (15) of the present paper. The parameter k_DR is the effective recombination coefficient, which we take as constant within the cold zone (we will stick with a conservative value of 1 × 10⁻⁶ cm³ s⁻¹). The bracketed part of Equation (22) is a measure of how much, relatively speaking, the electron density at (and within) r_en is reduced when considering dissociative recombination in addition to radial transport (in the limit of ${k}_{\mathrm{DR}}\to 0$ the bracketed expression equals 1). From G_j,j we compute g_i,j according to

$$\begin{eqnarray}&&{g}_{i,j}={G}_{j,j}\displaystyle \prod _{k=j-1}^{i-2}{\beta }_{k,j-1}^{* }.\end{eqnarray} \tag{ 23 }$$

The key parameters discussed in Section 4 are now, within this modified approach, instead given by

$$\begin{eqnarray}&&{n}_{I,i}=\tfrac{Q\nu {dr}}{4\pi {{ur}}_{i}^{2}}\left(\sum _{j=1}^{i}\,\tfrac{{c}_{i,j}}{{v}_{i,j}}+\tfrac{{d}_{i,j}+{g}_{i,j}}{{v}_{i,j}^{* }}\right),\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&\langle {u}_{I,i}\rangle =\displaystyle \frac{{\sum }_{j=1}^{i}\,{c}_{i,j}\,{v}_{i,j}+({d}_{i,j}\,+{g}_{i,j}){v}_{i,j}^{* }}{{\sum }_{j=1}^{i}{c}_{i,j}+{d}_{i,j}+{g}_{i,j}},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{R}_{\mathrm{dens},i}={\sum }_{j=1}^{i}\tfrac{{d}_{i,j}+(1-{f}_{\mathrm{XH}+}){g}_{i,j}}{{v}_{i,j}^{* }}/{\sum }_{j=1}^{i}\tfrac{{c}_{i,j}}{{v}_{i,j}},\end{eqnarray} \tag{ 26a }$$

and

$$\begin{eqnarray}&&{R}_{\mathrm{flux},i}=\tfrac{\sum _{j=1}^{i}{d}_{i,j}+(1-{f}_{\mathrm{XH}+}){g}_{i,j}}{{\sum }_{j=1}^{i}{c}_{i,j}}.\end{eqnarray} \tag{ 26b }$$

We have used this modified approach for perihelion conditions with input parameters as specified in Table 3 and with resulting key parameters shown in Figure 4 and discussed in the next paragraph. The outgassing rate Q was set to match neutral number densities measured by COPS on 2015 August 1 near r = 200 km with the assumption of radial expansion and a constant outflow speed of u = 1 km s⁻¹. We assume T_e = 10 eV for r > r_en = 100 km. The value set for r_en is taken from Figure 5 of Mandt et al. (2016). Finally, we assume f_XH+ = 0.5.

Zoom In Zoom Out Reset image size

The outputs shown in Figure 4 are from the four model types A, A*, B, and B*, as well as from a modified version of model A with the cross section for Equation (12) given by relation 10(b) instead of Equation (13). To no surprise, A* predicts the highest mean ion flow speeds and consequently the lowest ion number densities. The ion density profiles are rather well described by 1/r^γ relations for r > 110 km, and rough consistency may be achieved with small tuning of the assumed T_e and the assumed electron density profile. We find that while the models give n_I and $\langle {u}_{I}\rangle$ values at 200 km that are in reasonable agreement with limits constrained by LAP and MIP observations (typically ∼1000 cm⁻³ and a few kilometers per second, respectively), the H₃O⁺/H₂O⁺ number density and flux ratios are >5 and 7, respectively. With the extreme assumption of f_XH+ = 1, the H₃O⁺/H₂O⁺ ratios at 200 km become >3.7 (number density) and >5.1 (flux). As a reminder, DFMS observations near perihelion reveal varying H₃O⁺/H₂O⁺ ratios, not seldomly in the vicinity of unity and even below (see Figure 7 of Fuselier et al. 2016). Speculatively this may reflect nonradial ion trajectories influenced by electromagnetic forces (see Section 6 for a brief discussion) and/or plasma instabilities, with energization of the ion population by plasma waves (see Sections 5.3.2 and 6).

5.3.1. Temperature Effects

The parameter u_T has for all model runs presented hitherto been set to 500 m s⁻¹. As mentioned, u_T may be said to represent a thermal component of the mean relative speed in collisions. In a gas of molecules characterized by a temperature T the mean relative speed in collisions between molecules of types 1 and 2 is (see, e.g., Equation (24.10) in Atkins & de Paula 2002)

$$\begin{eqnarray}&&{u}_{T}=\sqrt{\tfrac{8{k}_{B}T}{\pi {\mu }_{\mathrm{1,2}}}},\end{eqnarray} \tag{ 27 }$$

with k_B being Boltzmann's constant and μ_1,2 the reduced mass of species 1 and 2 [m₁m₂/(m₁ + m₂)]. With H₂O as the molecular type 1 and either H₂O⁺ or H₃O⁺ as the type 2, a mean relative velocity of 500 m s⁻¹ corresponds to a temperature of ∼100 K. Considering a temperature as low as 20 K gives u_T in the vicinity of ∼200 m s⁻¹, while a higher temperature of 300 K gives u_T ≈ 840 m s⁻¹ (for theoretical estimates of kinetic temperatures in cometary comae, see, e.g., Tenishev et al. 2008). We have tested different values of u_T for runs with cross sections according to model A but with the cross section for Equation (12) given by relation 10(b) instead of Equation (13). In these tests, performed for both Q = 2 × 10²⁷ s⁻¹ and Q = 2 × 10²⁸ s⁻¹, we used a constant electric field of either 0.03 mV m⁻¹ or 0.1 mV m⁻¹. The values of n_I, $\langle {u}_{I}\rangle$, and R_dens at r = 200 km for different utilized values of u_T are summarized in Table 4. The parameter most influenced is the H₃O⁺/H₂O⁺ number density ratio, which drops markedly with increasing u_T. The reason for the decreased number density ratio with increased u_T is that the cross section for the proton-transfer reaction H₂O⁺ + H₂O H₃O⁺ + OH decreases with increased relative velocity between the reactants. This effect is most notable when considering the A and B models, as in these models a rather steep energy dependence is applied for the proton-transfer reaction.

Table 4.  Model Output at 200 km from Different Simulations as Described in Section 5.3.1

Input	nI (cm−3)	$\langle {u}_{I}\rangle$(m s−1)	Rdens	Input	nI (cm−3)	$\langle {u}_{I}\rangle$(m s−1)	Rdens
Q1, E1, 200	111	3807	8.2	Q2, E1, 200	3400	1042	588
Q1, E1, 500	104	4100	3.1	Q2, E1, 500	2861	1303	118
Q1, E1, 840	100	4250	2.0	Q2, E1, 840	2396	1640	45
Q1, E2, 200	52	8548	1.9	Q2, E2, 200	1727	3033	204
Q1, E2, 500	51	8651	1.1	Q2, E2, 500	1122	4443	32
Q1, E2, 840	51	8699	0.8	Q2, E2, 840	975	4976	15

Note. The columns labeled "Input" describe the difference between the models: Q1 is Q = 2 × 10²⁷ s⁻¹, Q2 is Q = 2 × 10²⁸ s⁻¹, E1 is E = 0.03 mV m⁻¹, and E2 is E = 0.1 mV m⁻¹. The number (200, 500, or 840) indicates the utilized value of u_T in m s⁻¹.

Download table as:  ASCIITypeset image

5.3.2. Ion Motion due to Waves

We have so far worked with the assumption that the ions and the neutrals are characterized by similar temperatures when calculating u_T. The motion of the ions may differ from neutral molecule motion for other reasons than the presence of electric fields or temperature considerations. In particular, low-frequency plasma waves and instabilities can provide ions with energy not directly accessible to the neutral gas. This need not be true for heating of the plasma by the wave, but the kinetic energy of the plasma motion inherent to the wave itself will for low-frequency waves reside mainly in the ions. The velocity average (with direction) in this motion is zero, but the mean energy and ion speed (magnitude of velocity) are not. We can let u_T represent this mean speed in the wave motion of the ions just as well as random thermal ion motion, as an average over a wave period.

While much remains to be done on plasma waves at 67P, some results have been published. Richter et al. (2015) reported large-amplitude (ΔB/B ∼ 1) magnetic field fluctuations during pre-perihelion conditions, also seen in hybrid simulations by Koenders et al. (2015, 2016). Though the exact nature of these waves remains to be detailed, their very low frequency means that ion velocities associated with them are likely small, as is also the case for the mirror mode waves reported by Volwerk et al. (2016). Gunell et al. (2017) used data from LAP and the Ion Composition Analyzer (ICA; Nilsson et al. 2007) to suggest that ion acoustic waves may be unstable to the ion distributions seen at the comet. Karlsson et al. (2017) identified clear spectral peaks at the lower hybrid frequency in the LAP electric field measurements. As strong plasma density gradients are observed such waves could be generated through the drift lower hybrid wave instability at very fast growth rates, on the order of a wave period. In the terrestrial magnetotail, Norgren et al. (2012) observed such waves to acquire about 10% of the energy density available in the electron gas. This energy mostly resides in the ion motion. It may be noted that while Edberg et al. (2015) found the average plasma density during a Rosetta close flyby of the 67P nucleus in 2015 February to be consistent with a 1/r profile (as mentioned in Section 1), the data showed large variations (Δn/n ∼ 1) around the mean profile. There clearly are reasons to consider a situation in which a substantial fraction of the electron thermal energy is transferred to ion kinetic energy owing to ion motion in the wave.

In Table 5, we therefore display the sensitivity of u_T to results from model runs similar to the ones described in Section 5.2 and Figure 4 (and its caption), but this time with an ambipolar electric field of magnitude E(r) = U₀/r with U₀ = 5 V, for r > r_en = 100 km. Results are displayed for u_T of 500 m s⁻¹ (our default value) and 3000 m s⁻¹, corresponding to ∼1 eV mean ion kinetic energy over a wave period. Again, among the listed parameters the H₃O⁺/H₂O⁺ number density ratio, R_dens, is the one most influenced, and the reduction becomes, of course, even more pronounced if considering even larger u_T. An interesting by-product from these simulations is that the more energized the ion population one considers for r > r_en, the more the H₃O⁺/XH⁺ flux ratio will resemble the one at r_en. This suggests that the effect of plasma waves on observables like the H₃O⁺/H₂O⁺ number density ratio and other ion number density ratios can be quite dramatic. However, if treating u_T as a velocity vector perpendicular to the radial direction, we need to modify the regional absorption and survival probabilities (see Section 2 and Table 1) since when an ion travels a distance dr in the radial direction it will also travel a distance u_T/u_I,rad dr in the perpendicular direction, where u_I,rad is the velocity of the ion in the radial direction (given by the v_i,j or ${v}_{i,{j}^{* }}$ arrays). The total distance traveled through the region will then amount to ds = (1 + (u_T/u_I,rad)²)^1/2dr, and the exponents in the expressions for the regional absorption and survival probabilities (see Table 1) will need to be multiplied by ds/dr (note that ds is a variable). The influence of this modification is revealed by comparisons of data in Table 4 with data in Table 6, as well as by comparisons of data in Table 5 with data in Table 7. For u_T = 500 m s⁻¹ (the default value used when generating the data displayed in Figures 1–5) the effect is very small; the R_dens values, which are the most affected, increase only by ∼10% at Q = 2 × 10²⁸ s⁻¹, and less so for lower activity levels. However, as expected, the influence becomes more pronounced when considering higher u_T. As seen in Table 7, the R_dens values at 200 km are still markedly decreased when increasing u_T to 3 km s⁻¹, but they fall below 3, and only marginally so, only in the A* model and when utilizing F_XH+ = 1.

Zoom In Zoom Out Reset image size

Table 5.  Results from Model Runs at 200 km Similar to the Ones Described in the Caption of Figure 4 but with U₀ Set to 5 V Instead of 10 V and with u_T Set as Specified

	uT = 0.5 km s−1			uT = 3 km s−1
Model Type	nI (103 cm−3)	$\langle {u}_{I}\rangle$ (km s−1)	Rdens	nI (103 cm−3)	$\langle {u}_{I}\rangle$ (km s−1)	Rdens
A	1.43	2.53	49.5 (38.3)	1.19	3.07	2.96 (2.17)
A*	1.45	2.50	7.67 (5.79)	1.19	3.05	1.91 (1.33)
B	1.90	1.77	65.2 (49.3)	1.43	2.46	3.42 (2.44)
B*	1.91	1.76	10.2 (7.55)	1.43	2.45	2.15 (1.45)
Special	2.49	1.24	85.0 (63.7)	1.30	2.77	3.29 (2.38)

Note. Values within parentheses are for simulation runs wherein F_XH+ was set to 1 instead of 0.5 (with our applied assumptions this affects only R_dens among the tabulated parameters).

Download table as:  ASCIITypeset image

Table 6.  Same as in Table 4 but with Consideration of Extended Pathways through Regions Following a Treatment of the Velocity u_T as Perpendicular to the Radial Direction

Input	nI (cm−3)	$\langle {u}_{I}\rangle$(m s−1)	Rdens	Input	nI (cm−3)	$\langle {u}_{I}\rangle$(m s−1)	Rdens
Q1, E1, 200	112	3798	8.4	Q2, E1, 200	3404	1041	598
Q1, E1, 500	104	4077	3.3	Q2, E1, 500	2910	1274	132
Q1, E1, 840	101	4210	2.2	Q2, E1, 840	2502	1552	57
Q1, E2, 200	52	8545	2.0	Q2, E2, 200	1749	2995	211
Q1, E2, 500	51	8642	1.1	Q2, E2, 500	1148	4359	36
Q1, E2, 840	51	8684	0.9	Q2, E2, 840	1007	4847	18

Download table as:  ASCIITypeset image

Table 7.  Same as in Table 5 but with Consideration of Extended Pathways through Regions Following a Treatment of the Velocity u_T as Perpendicular to the Radial Direction

	uT = 0.5 km s−1			uT = 3 km s−1
Model Type	nI (103 cm−3)	$\langle {u}_{I}\rangle$ (km s−1)	Rdens	nI (103 cm−3)	$\langle {u}_{I}\rangle$ (km s−1)	Rdens
A	1.45	2.51	55.0 (42.6)	1.30	2.79	6.36 (4.82)
A*	1.47	2.47	8.31 (6.28)	1.31	2.77	3.82 (2.82)
B	1.93	1.74	72.7 (55.0)	1.64	2.10	7.98 (5.90)
B*	1.94	1.73	11.2 (8.28)	1.65	2.09	4.74 (3.41)
Special	2.53	1.21	94.5 (71.1)	1.50	2.37	7.37 (5.50)

Download table as:  ASCIITypeset image

In Figure 5 we show IEDFs (dot-dashed lines) and the relative contribution to the fluxes of H₂O⁺ (solid lines) and H₃O⁺ (dashed lines) at r = 200 km calculated from model A with Q values of 2 × 10²⁶ s⁻¹ (black lines), 2 × 10²⁷ s⁻¹ (red lines), and 2 × 10²⁸ s⁻¹ (blue lines). For these runs we utilized constant electric fields of 0.03 mV m⁻¹. As expected, the functions associated with lower Q are associated with higher ion energies and a higher relative contribution by H₂O⁺ ions. Ion energy distributions are assessable from measurements by the ICA (Nilsson et al. 2007, 2015) and the Ion Electron Sensor (IES; Burch et al. 2007; Mandt et al. 2016). At low energies the measurements are affected by the spacecraft potential (typically ∼10–20 V negative; e.g., Odelstad et al. 2015; Eriksson et al. 2017), and so model–observation comparisons require detailed treatment, postponed for future studies.

We show here that if assuming Maxwellian energy distributions for the neutrals and the ions, the use of the cross-section relations 10(b) and 10(c) for reaction 9(b) at E_rel < 1.5 eV gives rate coefficients at 300 K that differ somewhat from experimental values. For binary reactions, energy-dependent collision cross sections relate to temperature-dependent rate coefficients, k_1,2(T), according to (e.g., Chapter 27.1 in Atkins & de Paula 2002)

$$\begin{eqnarray}&&{k}_{\mathrm{1,2}}(T)=\sqrt{\tfrac{8}{\pi {\mu }_{\mathrm{1,2}}}}\tfrac{1}{{k}_{B}T}{\int }_{0}^{\infty }{E}_{\mathrm{rel}}\sigma ({E}_{\mathrm{rel}})\exp \left(-\tfrac{{E}_{\mathrm{rel}}}{{k}_{B}T}\right){{dE}}_{\mathrm{rel}},\end{eqnarray} \tag{ 28 }$$

where k_B is Boltzmann's constant and μ_1,2 is the reduced mass of the reactants. The experimentally determined value for Equation 9(b) at room temperature is 2.1 × 10⁻⁹ cm³ s⁻¹ (Huntress & Pinizzotto 1973), while insertion of Equation 10(b) into Equation (28) gives a value of 7.4 × 10⁻⁹ cm³ s⁻¹ at 300 K, higher than the experimental one by a factor of ∼3.5. Likewise, doing the same for Equation 10(c) results in 1.4 × 10⁻⁹ cm³ s⁻¹, ∼34% lower than the experimental value. Obtaining the experimentally determined rate coefficient at 300 K and a temperature dependence (T/300)^−0.5 as suggested in the UMIST database for astrochemistry (see Woodall et al. 2007), requires at low collision energies an energy-dependent cross section of approximately

$$\begin{eqnarray}&&{\sigma }_{\mathrm{pt}}({E}_{\mathrm{rel}})=6{E}_{\mathrm{rel}}^{-1},\end{eqnarray} \tag{ 29 }$$

where insertion of E_rel in eV gives the cross section in 10⁻¹⁶ cm². Such a function intersects with Equation 10(c) at E_rel = 0.04 eV. Table 8 reveals the effect on predicted key parameters for model A* runs modified with the introduction of the cross-section relation (29) for E_rel < 0.04 eV. It is noted that using Equation (28) for E_rel < 0.04 and Equation 10(c) for higher E_rel (up to 1.5 eV) as the cross section for the proton-transfer reaction 9(b) gives rate coefficients within 33% agreement with the formula 2.1 × 10⁻⁹ × (T/300)^−0.5 up to temperatures of 1000 K. From Table 8 it is seen that five out of six considered models give roughly similar values of n_I and $\langle {u}_{I}\rangle$, but with marked variations in R_dens. The model labeled A_special is the one used for the data presented in Table 5 and yields, as seen in Table 6, the lowest values of $\langle {u}_{I}\rangle$. It is noted that in this model the reaction (12), i.e., the proton transfer from H₃O⁺ to H₂O, is given the cross-section relation 10(b) regardless of energy. Similar to what was noted above, applying such a cross-section relation at very low energies yields a 300 K rate coefficient significantly higher than derived experimentally by Smith et al. (1980). In summary, Equation 10(b) gives overestimated cross sections at very low energies, implying that models A and B will tend to yield overestimated H₃O⁺/H₂O⁺ number density ratios. Moreover, for a given field strength, models A and B, which in addition are modified to make H₃O⁺ proton transfer to H₂O equally efficient as the H₂O⁺ proton transfer to H₂O ("special models of A and B type"), will tend to overestimate the collisional coupling and therefore underestimate the mean ion drift speed.

Table 8.  Model Outputs for r = 200 km for Simulations with Constant Electric Field of Magnitude 0.03 mV m⁻¹ and for Q Values as Specified

Q (s−1)	A*	${{{\rm{A}}}^{* }}_{\mathrm{new}}$	A
2 × 1027	932; 4.57; 0.87	929; 4.59; 1.00	933; 4.60; 2.72
2 × 1028	1634; 2.60; 11.0	1630; 2.61; 16.1	1617; 2.62; 66.5
2 × 1029	31634; 1.13; 181	31632; 1.13; 315	31628; 1.13; 898
Q (s−1)	${{{\rm{A}}}^{* }}_{\mathrm{special}}$	${{{\rm{A}}}^{* }}_{\mathrm{new},\mathrm{special}}$	Aspecial
2 × 1027	931; 4.58; 0.87	929; 4.59; 1.00	1036; 4.10; 3.12
2 × 1028	1630; 2.61; 10.9	1672; 2.55; 16.5	2861; 1.30; 118
2 × 1029	31588; 1.13; 181	33158; 1.07; 331	34865; 1.01; 990

Note. The u and u_T parameters were set to 1 and 0.5 km s⁻¹, respectively. An array such as "932; 4.57; 0.87" refers to the calculated values of n_I (cm⁻³), $\langle {u}_{I}\rangle$ (km s⁻¹), and R_dens, respectively. Model descriptions: models A and A* are as described in Table 2. Model ${{{\rm{A}}}^{* }}_{\mathrm{new}}$ is similar to A* but uses for reaction 9(b) the cross-section relation (29) for E_rel < 0.04 eV. In models indexed by "special," reaction (12) is given the same cross-section relation as reaction 9(b) for the considered model case.

Download table as:  ASCIITypeset image

Vigren et al. (2016) reported an on average good agreement between modeled and observed electron number densities when Rosetta was at r ∼ 28 km and a heliocentric distance of ∼2.5 au in 2015 January 9–11. Their model was essentially Equation (14) with Q/u constrained from COPS measurements; so with u taken as ∼500 m s⁻¹, Q ≈ 2 × 10²⁶ s⁻¹. The good agreement (see Figure 1 in Vigren et al. 2016) was argued to be indicative of the validity of central assumptions in the model, including the one of ions being coupled to the neutrals all the way to the spacecraft. The latter is not consistent with results from the TRAIN model, as discussed below.

When applying the TRAIN model to Q = 2 × 10²⁶ s⁻¹, a neutral outflow velocity of 500 m s⁻¹, and weak electric fields, the neutrals are unable to interrupt acceleration of ions even well within cometocentric distances of 30 km. For example—with model A but with the use of an enhanced cross section for reaction (12) (see Section 5.1.2)—the use of a constant electric field of 0.03 mV m⁻¹ (0.1 mV m⁻¹) gives a mean ion flow speed at r = 28 km of ∼0.9 km s⁻¹ (2.2 km s⁻¹). For an ambipolar electric field of magnitude U₀/r (with U₀ = 5 V) for r > 20 km (for r > 10 km), the mean ion flow speed at r = 28 km becomes ∼3.2 km s⁻¹ (4.2 km s⁻¹). As electron temperatures on the order of 5–10 eV were measured by LAP in 2015 January, a strong ambipolar field could be expected, but from the observations presented in Vigren et al. (2016) it appears as if u_I ≈ u.

The results from the TRAIN model indicate that such a low mean ion velocity cannot be due to frequent collisions with the neutrals, so a more probable explanation would be that the plasma was strongly influenced by the magnetic field environment and/or plasma waves, as discussed in Section 5.3.2 above. It should be noted that our treatment of the ambipolar electric field so far has been with an implicit assumption of either a very weak or absent magnetic field or one in the radial direction (i.e., parallel to the envisioned electric field). Given a moderately strong magnetic field at moderately high angle to the radial direction, electrons would be bound closer to the nucleus, effectively reducing the ambipolar electric field. This can to some extent be explored by hybrid simulations (e.g., Koenders et al. 2015), though due to their limited possibilities for modeling electrons, full particle-in-cell (PIC) simulations, not yet available, may be needed. A detailed investigation of magnetic field effects is beyond the scope of the present study.