# Constraining ion transport in the diamagnetic cavity of comet 67P

Z. M. Lewis, A. Beth, M. Galand, P. Henri, M. Rubin, P. Stephenson

► To cite this version:

Z. M. Lewis, A. Beth, M. Galand, P. Henri, M. Rubin, et al. Constraining ion transport in the diamagnetic cavity of comet 67P. Monthly Notices of the Royal Astronomical Society, 2024, 530, pp.66-81. 10.1093/mnras/stae856. insu-04551237

HAL Id: insu-04551237

https://insu.hal.science/insu-04551237

Submitted on 18 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Constraining ion transport in the diamagnetic cavity of comet 67P

Zoe Lewis, Arnaud Beth, Marina Galand, Pierre Henri, Martin Rubin, Peter

Stephenson

► To cite this version:

Zoe Lewis, Arnaud Beth, Marina Galand, Pierre Henri, Martin Rubin, et al. Constraining ion transport in the diamagnetic cavity of comet 67P. Monthly Notices of the Royal Astronomical Society, 2024, 530 (1), pp.66-81. 10.1093/mnras/stae856. insu-04551237

HAL Id: insu-04551237

https://insu.hal.science/insu-04551237v1

Submitted on 18 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Constraining ion transport in the diamagnetic cavity of comet 67P

Z. M. Lewis<sup>1,1*</sup> A. Beth<sup>1</sup> M. Galand<sup>1</sup> P. Henri<sup>2,3</sup> M. Rubin<sup>4</sup> and P. Stephenson<sup>5</sup>

Accepted 2024 March 22. Received 2024 March 21; in original form 2024 February 2

## ABSTRACT

The European Space Agency Rosetta mission escorted comet 67P for a 2-yr section of its six and a half-year orbit around the Sun. By perihelion in 2015 August, the neutral and plasma data obtained by the spacecraft instruments showed the comet had transitioned to a dynamic object with large-scale plasma structures and a rich ion environment. One such plasma structure is the diamagnetic cavity: a magnetic field-free region formed by interaction between the unmagnetized cometary plasma and the impinging solar wind. Within this region, unexpectedly high ion bulk velocities have been observed, thought to have been accelerated by an ambipolar electric field. We have developed a 1D numerical model of the cometary ionosphere to constrain the impact of various electric field profiles on the ionospheric density profile and ion composition. In the model, we include three ion species: H<sub>2</sub>O<sup>+</sup>, H<sub>3</sub>O<sup>+</sup>, and NH<sub>4</sub><sup>+</sup>. The latter, not previously considered in ionospheric models including acceleration, is produced through the protonation of NH<sub>3</sub> and only lost through ion-electron dissociative recombination, and thus particularly sensitive to the time-scale of plasma loss through transport. We also assess the importance of including momentum transfer when assessing ion composition and densities in the presence of an electric field. By comparing simulated electron densities to Rosetta Plasma Consortium data sets, we find that to recreate the plasma densities measured inside the diamagnetic cavity near perihelion, the model requires an electric field proportional to r<sup>-1</sup> of around 0.5–2 mV m<sup>-1</sup> surface strength, leading to bulk ion speeds at Rosetta of 1.2–3.0 km s<sup>-1</sup>.

**Key words:** plasmas – comets: individual: comet 67P/CG.

## 1 INTRODUCTION

The cometary ionosphere is formed when the neutral gas coma becomes partially ionized by solar extreme ultraviolet (EUV) photons, as well as through electron-impact ionization and charge exchange with the solar wind (Cravens 1987). Unbounded by gravity, both neutral and ion populations escape to space, and have a constant source from the sublimation of the ices from the (sub-) surface of the nucleus. During its 2-yr escort of comet 67P, the Rosetta spacecraft witnessed the evolution of this cometary environment through a range of heliocentric distances (from 3.6 to 1.25 au pre-perihelion and up to 3.8 au post-perihelion). The plasma environment was probed by the Rosetta Plasma Consortium (RPC) instruments (Carr et al. 2007) and the Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA; Balsiger et al. 2007).

Near perihelion in 2015 August, the comet outgassing was at its maximum (>10<sup>28</sup> s<sup>-1</sup>; Hansen et al. 2016), leading to a more complex ion composition (Heritier et al. 2017b; Beth et al. 2022). The denser coma allows H<sub>2</sub>O<sup>+</sup> (directly produced through ionization) to readily transfer a proton to the neutral H<sub>2</sub>O to produce H<sub>3</sub>O<sup>+</sup>, which is

often the dominant ion species in the coma (Altwegg et al. 1993; Vigren & Galand 2013). H<sub>3</sub>O<sup>+</sup> can then be lost to reactions with high proton affinity (PA) neutrals (Heritier et al. 2017b) or, often to a lesser extent, by dissociative ion-electron recombination. The former process can happen repeatedly for neutrals with increasing PA until the terminal ion is reached. This is NH<sub>4</sub><sup>+</sup>, formed through the protonation of NH<sub>3</sub> that has the highest PA of the cometary neutrals (Altwegg et al. 1993; Vigren & Galand 2013). Unambiguous detection of NH<sub>4</sub><sup>+</sup> was first reported by Beth et al. (2016), using the ROSINA/ Double Focussing Mass Spectrometer instrument (DFMS; Balsiger et al. 2007). Near perihelion, detections of NH<sub>4</sub><sup>+</sup> increase (Lewis et al. 2023), sometimes overtaking H<sub>3</sub>O<sup>+</sup> at times of high outgassing.

From late 2015 April at 1.8 au to 2016 February at 2.4 au, a diamagnetic cavity was sporadically detected as regions of negligible magnetic field (Goetz et al. 2016a). The diamagnetic cavity represents the region where the mass-loaded solar wind plasma is diverted away from the unmagnetized nucleus, preventing the interplanetary magnetic field from penetrating to the surface, and leaving an inner magnetic field free region around the nucleus. This region, previously predicted by theory (Biermann, Brosowski & Schmidt 1967), was first detected by the European Space Agency (ESA) Giotto spacecraft during the 1986 fly-by of comet 1P/Halley (hereafter 1P; Neubauer

© 2024 The Author(s).

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

<footnote n="1"><sup>1</sup>Department of Physics, Imperial College London, London, SW7 2AZ, UK</footnote>

<footnote n="2"><sup>2</sup>Laboratoire de Physique et Chimie de l'Environnement et de l'Espace (LPC2E), CNRS, Université d'Orléans, Orléans, 45071, France</footnote>

<footnote n="3"><sup>3</sup>Laboratoire Lagrange, OCA, UCA, CNRS, Nice, 34229, France</footnote>

<footnote n="4"><sup>4</sup>Space Research & Planetary Sciences, Physics Institute, University of Bern, Bern, CH-3012, Switzerland</footnote>

<footnote n="5"><sup>5</sup>Lunar and Planetary Laboratory, University of Arizona, AZ 85721, USA</footnote>

<footnote>*E-mail: z.lewis21@imperial.ac.uk</footnote>
# Constraining ion transport in the diamagnetic cavity of comet 67P

Z. M. Lewis<sup>1,1*</sup> A. Beth<sup>1</sup> M. Galand<sup>1</sup> P. Henri<sup>2,3</sup> M. Rubin<sup>4</sup> and P. Stephenson<sup>5</sup>

Accepted 2024 March 22. Received 2024 March 21; in original form 2024 February 2

## ABSTRACT

The European Space Agency Rosetta mission escorted comet 67P for a 2-yr section of its six and a half-year orbit around the Sun. By perihelion in 2015 August, the neutral and plasma data obtained by the spacecraft instruments showed the comet had transitioned to a dynamic object with large-scale plasma structures and a rich ion environment. One such plasma structure is the diamagnetic cavity: a magnetic field-free region formed by interaction between the unmagnetized cometary plasma and the impinging solar wind. Within this region, unexpectedly high ion bulk velocities have been observed, thought to have been accelerated by an ambipolar electric field. We have developed a 1D numerical model of the cometary ionosphere to constrain the impact of various electric field profiles on the ionospheric density profile and ion composition. In the model, we include three ion species: $H_2O^+$, $H_3O^+$, and $NH_4^+$. The latter, not previously considered in ionospheric models including acceleration, is produced through the protonation of $NH_3$ and only lost through ion-electron dissociative recombination, and thus particularly sensitive to the time-scale of plasma loss through transport. We also assess the importance of including momentum transfer when assessing ion composition and densities in the presence of an electric field. By comparing simulated electron densities to Rosetta Plasma Consortium data sets, we find that to recreate the plasma densities measured inside the diamagnetic cavity near perihelion, the model requires an electric field proportional to $r^{-1}$ of around 0.5–2 mV m$^{-1}$ surface strength, leading to bulk ion speeds at Rosetta of 1.2–3.0 km s$^{-1}$.

**Key words:** plasmas – comets: individual: comet 67P/CG.

## 1 INTRODUCTION

The cometary ionosphere is formed when the neutral gas coma becomes partially ionized by solar extreme ultraviolet (EUV) photons, as well as through electron-impact ionization and charge exchange with the solar wind (Cravens 1987). Unbounded by gravity, both neutral and ion populations escape to space, and have a constant source from the sublimation of the ices from the (sub-) surface of the nucleus. During its 2-yr escort of comet 67P, the Rosetta spacecraft witnessed the evolution of this cometary environment through a range of heliocentric distances (from 3.6 to 1.25 au pre-perihelion and up to 3.8 au post-perihelion). The plasma environment was probed by the Rosetta Plasma Consortium (RPC) instruments (Carr et al. 2007) and the Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA; Balsiger et al. 2007).

Near perihelion in 2015 August, the comet outgassing was at its maximum (>10<sup>28</sup> s<sup>-1</sup>; Hansen et al. 2016), leading to a more complex ion composition (Heritier et al. 2017b; Beth et al. 2022). The denser coma allows $H_2O^+$ (directly produced through ionization) to readily transfer a proton to the neutral $H_2O$ to produce $H_3O^+$, which is

often the dominant ion species in the coma (Altwegg et al. 1993; Vigren & Galand 2013). $H_3O^+$ can then be lost to reactions with high proton affinity (PA) neutrals (Heritier et al. 2017b) or, often to a lesser extent, by dissociative ion-electron recombination. The former process can happen repeatedly for neutrals with increasing PA until the terminal ion is reached. This is $NH_4^+$, formed through the protonation of $NH_3$ that has the highest PA of the cometary neutrals (Altwegg et al. 1993; Vigren & Galand 2013). Unambiguous detection of $NH_4^+$ was first reported by Beth et al. (2016), using the ROSINA/ Double Focussing Mass Spectrometer instrument (DFMS; Balsiger et al. 2007). Near perihelion, detections of $NH_4^+$ increase (Lewis et al. 2023), sometimes overtaking $H_3O^+$ at times of high outgassing.

From late 2015 April at 1.8 au to 2016 February at 2.4 au, a diamagnetic cavity was sporadically detected as regions of negligible magnetic field (Goetz et al. 2016a). The diamagnetic cavity represents the region where the mass-loaded solar wind plasma is diverted away from the unmagnetized nucleus, preventing the interplanetary magnetic field from penetrating to the surface, and leaving an inner magnetic field free region around the nucleus. This region, previously predicted by theory (Biermann, Brosowski & Schmidt 1967), was first detected by the European Space Agency (ESA) Giotto spacecraft during the 1986 fly-by of comet 1P/Halley (hereafter 1P; Neubauer

© 2024 The Author(s).

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

<footnote n="1"><sup>1</sup>Department of Physics, Imperial College London, London, SW7 2AZ, UK</footnote>

<footnote n="2"><sup>2</sup>Laboratoire de Physique et Chimie de l'Environnement et de l'Espace (LPC2E), CNRS, Université d'Orléans, Orléans, 45071, France</footnote>

<footnote n="3"><sup>3</sup>Laboratoire Lagrange, OCA, UCA, CNRS, Nice, 34229, France</footnote>

<footnote n="4"><sup>4</sup>Space Research & Planetary Sciences, Physics Institute, University of Bern, Bern, CH-3012, Switzerland</footnote>

<footnote n="5"><sup>5</sup>Lunar and Planetary Laboratory, University of Arizona, AZ 85721, USA</footnote>

<footnote>*E-mail: z.lewis21@imperial.ac.uk</footnote>
et al. 1986). However, even at perihelion, 67P had a ≳10 times lower outgassing rate than 1P did during the Giotto fly-by. Consequently, the plasma environment of 67P was found to be substantially different; the diamagnetic cavity boundary location appeared to be highly variable and was inconsistent with a well-established global structure (Goetz et al. 2016a; Henri et al. 2017). This was previously predicted by magnetohydrodynamic modelling by Rubin et al. (2012), which showed that asymmetrical outgassing could lead to unstable plasma flow and Kelvin–Helmholtz instabilities along the cavity boundary.

Since cometary ions are produced from parent neutrals that are expanding from the surface at $u_n$, a first assumption would be that the ions travel with the same bulk velocity. Such an assumption has been shown to accurately reproduce the total plasma density measured by the RPC Mutual Impedance Probe (RPC-MIP) and Langmuir Probe (RPC-LAP) later in the escort phase, when the comet activity was low and Rosetta was within a few tens of kilometres of the nucleus (Vigren & Eriksson 2017). However, the same model considerably overestimates the plasma density closer to perihelion, suggesting that the ions may be undergoing significant acceleration that does not effect the neutrals.

The ion bulk speed has been derived near perihelion from Rosetta data in several ways. One method (Vigren et al. 2017) involves the use of a combination of the current-voltage characteristics from the Langmuir probe (RPC-LAP) instrument, and electron number densities from the RPC-MIP (Trotignon et al. 2007). Applied to a 3-day range in 2015 August by Vigren et al. (2017) and in 2015 November by Odelstad et al. (2018), effective ion velocities in the range 2–8 km s⁻¹ at 200 and 135 km were derived, respectively. This range is higher than the 0.5–1 km s⁻¹ neutral speed (Hansen et al. 2016), and led to the conclusion that the ions and neutrals are indeed decoupled. However, this method presents an upper estimate of the ion bulk speed, since it does not allow separation from the thermal speed (Mott-Smith & Langmuir 1926). For an ion drift speed of the order of the neutral speed, 1 km s⁻¹, the thermal speed term is of comparable size for ion temperatures ≳0.1 eV. The ion temperature is not well constrained, but values in the range 0.7–1.6 eV have been derived by Bergman et al. (2021b) using measurements from the Ion Composition Analyser (RPC-ICA; Nilsson et al. 2007). Another method applied by Vigren et al. (2017) is based on a simple flux conservation model, with the assumption of radial outflow, to estimate the ion speed. They found that the two methods produce similar values of the ion bulk speed, but note that EUV attenuation was neglected, as well as dissociative recombination. Neglecting the decrease in the photoionization frequency results in overestimated values of the ion velocity supported by Johansson et al. (2017) and their investigation regarding the attenuation due to dust. A third method uses RPC-ICA energy spectra fitted to drifting Maxwell–Boltzmann distributions for data from 88 diamagnetic cavity crossings (Bergman et al. 2021b). This fitting process allows the thermal and drift velocity components to be separated if the instrumental energy resolution allows, and a bulk speed of 5–10 km s⁻¹ was found, with a peak probability at 7 km s⁻¹ though the ability of RPC-ICA to measure the low-energy ion population is still debated.

A key challenge when interpreting these calculated ion bulk speed measurements is the influence of the negative spacecraft potential (Odelstad et al. 2017) on the detected ions. RPC-LAP was positioned on a boom, and is therefore less strongly affected by this problem than RPC-ICA, which was positioned on the main body of the spacecraft. Low-energy, positively charged ions were accelerated towards the detector, resulting in the distortion of their energy spectra

as measured by the instrument. Studies have been performed (e.g. Bergman et al. 2020; Johansson et al. 2020) using the Spacecraft-Plasma Interaction System (SPIS) to model this effect with the aim of correcting for it in the data. Bergman et al. (2021a) compared their derived ion temperature and velocity with results from the SPIS model to conclude that the effect on the measured speed would be fairly minimal (shifting the peak of the probability density function for the ion drift speed ($u_{i,D}$) to 6.9 km s⁻¹). However, other uncertainties remain, for example the variation in the spacecraft potential and the flow direction of the ions. The spacecraft potential distorts the measured direction of low-energy ions from ICA, but higher energy cometary ions have been mostly observed in the radial direction (Nilsson et al. 2017). Particle-in-cell simulations have been used to correct for the spacecraft potential in ICA data by Bergman et al. (2021b), showing a surprising ‘backstream’ of ions towards the nucleus.

Despite these uncertainties, it is likely that the cometary ions within the diamagnetic cavity are indeed travelling faster than the neutral population, and are therefore collisionally uncoupled to some extent. Since the ionizing photons and electrons carry a negligible momentum compared to the neutrals, cometary ions are produced at the neutral speed, but are then susceptible to electromagnetic fields. Collisional decoupling from the neutral flow arises when the ions are accelerated beyond the neutral speed. Determination of the extent of this decoupling is important to understand how the cometary plasma behaves (including formation of the diamagnetic cavity boundary), but it is currently not well understood.

Outside the diamagnetic cavity and solar wind ion cavity, cometary ions are accelerated by the convective electric field of the solar wind, which leads to the mass loading and ion pick up process (Szegő et al. 2000; Behar et al. 2016). Inside the unmagnetized diamagnetic cavity region, the $\vec{J} \times \vec{B}$ term in Ohm's law vanishes, and the ambipolar electric field becomes dominant. Unlike at planets, the electron pressure gradient does not arise from the different gravitational forces acting on ions and electrons (e.g. Schunk & Nagy 2009), but from the higher energy and pressure of the cometary electrons compared to the ions due to their smaller mass. An electric field arises to oppose the creation of a charge imbalance, in order to satisfy quasi-neutrality on scales longer than the Debye length. At comet 1P during the Giotto fly-by, Gan & Cravens (1990) demonstrated that the ambipolar field (and therefore the electron pressure) was negligible, since the dense neutral coma led to efficient collisional cooling of electrons. At 67P, for low outgassing activity, the ambipolar electric field has been shown to have a tangible impact on the cometary plasma environment. One such impact is keeping cometary electrons ‘trapped’ in a region close to the nucleus (e.g. Deca et al. 2019; Stephenson et al. 2022), thereby increasing the efficiency of electron-neutral collisions and leading to a cold (though minor) electron population (<0.1 eV) consistent with observations by RPC-LAP and RPC-MIP (e.g. Eriksson et al. 2017; Engelhardt et al. 2018; Gilet et al. 2020; Wattiaux et al. 2020). In addition, solar wind electrons were found to have been accelerated towards the nucleus by such an electric field, detected as a suprathermal electron population (Deca et al. 2017; Madanian et al. 2017) and responsible for generating aurora (Galand et al. 2020; Stephenson et al. 2021) and ionization (Stephenson et al. 2023).

The magnitude and functional form of the ambipolar electric field are difficult to constrain. Cravens (2004) derived an expression for the radial electric field ($\vec{E} = E(r) \vec{r}$) by first enforcing $n_e \approx n_i \propto 1/r$ while $T_e$ is assumed constant and isotropic throughout the coma. This arises from a field-free and chemistry-free model of the coma,
et al. 1986). However, even at perihelion, 67P had a ≳10 times lower outgassing rate than 1P did during the Giotto fly-by. Consequently, the plasma environment of 67P was found to be substantially different; the diamagnetic cavity boundary location appeared to be highly variable and was inconsistent with a well-established global structure (Goetz et al. 2016a; Henri et al. 2017). This was previously predicted by magnetohydrodynamic modelling by Rubin et al. (2012), which showed that asymmetrical outgassing could lead to unstable plasma flow and Kelvin–Helmholtz instabilities along the cavity boundary.

Since cometary ions are produced from parent neutrals that are expanding from the surface at $u_n$, a first assumption would be that the ions travel with the same bulk velocity. Such an assumption has been shown to accurately reproduce the total plasma density measured by the RPC Mutual Impedance Probe (RPC-MIP) and Langmuir Probe (RPC-LAP) later in the escort phase, when the comet activity was low and Rosetta was within a few tens of kilometres of the nucleus (Vigren & Eriksson 2017). However, the same model considerably overestimates the plasma density closer to perihelion, suggesting that the ions may be undergoing significant acceleration that does not effect the neutrals.

The ion bulk speed has been derived near perihelion from Rosetta data in several ways. One method (Vigren et al. 2017) involves the use of a combination of the current-voltage characteristics from the Langmuir probe (RPC-LAP) instrument, and electron number densities from the RPC-MIP (Trotignon et al. 2007). Applied to a 3-day range in 2015 August by Vigren et al. (2017) and in 2015 November by Odelstad et al. (2018), effective ion velocities in the range 2–8 km s⁻¹ at 200 and 135 km were derived, respectively. This range is higher than the 0.5–1 km s⁻¹ neutral speed (Hansen et al. 2016), and led to the conclusion that the ions and neutrals are indeed decoupled. However, this method presents an upper estimate of the ion bulk speed, since it does not allow separation from the thermal speed (Mott-Smith & Langmuir 1926). For an ion drift speed of the order of the neutral speed, 1 km s⁻¹, the thermal speed term is of comparable size for ion temperatures ≳0.1 eV. The ion temperature is not well constrained, but values in the range 0.7–1.6 eV have been derived by Bergman et al. (2021b) using measurements from the Ion Composition Analyser (RPC-ICA; Nilsson et al. 2007). Another method applied by Vigren et al. (2017) is based on a simple flux conservation model, with the assumption of radial outflow, to estimate the ion speed. They found that the two methods produce similar values of the ion bulk speed, but note that EUV attenuation was neglected, as well as dissociative recombination. Neglecting the decrease in the photoionization frequency results in overestimated values of the ion velocity supported by Johansson et al. (2017) and their investigation regarding the attenuation due to dust. A third method uses RPC-ICA energy spectra fitted to drifting Maxwell–Boltzmann distributions for data from 88 diamagnetic cavity crossings (Bergman et al. 2021b). This fitting process allows the thermal and drift velocity components to be separated if the instrumental energy resolution allows, and a bulk speed of 5–10 km s⁻¹ was found, with a peak probability at 7 km s⁻¹ though the ability of RPC-ICA to measure the low-energy ion population is still debated.

A key challenge when interpreting these calculated ion bulk speed measurements is the influence of the negative spacecraft potential (Odelstad et al. 2017) on the detected ions. RPC-LAP was positioned on a boom, and is therefore less strongly affected by this problem than RPC-ICA, which was positioned on the main body of the spacecraft. Low-energy, positively charged ions were accelerated towards the detector, resulting in the distortion of their energy spectra

as measured by the instrument. Studies have been performed (e.g. Bergman et al. 2020; Johansson et al. 2020) using the Spacecraft-Plasma Interaction System (SPIS) to model this effect with the aim of correcting for it in the data. Bergman et al. (2021a) compared their derived ion temperature and velocity with results from the SPIS model to conclude that the effect on the measured speed would be fairly minimal (shifting the peak of the probability density function for the ion drift speed ($u_{i,D}$) to 6.9 km s⁻¹). However, other uncertainties remain, for example the variation in the spacecraft potential and the flow direction of the ions. The spacecraft potential distorts the measured direction of low-energy ions from ICA, but higher energy cometary ions have been mostly observed in the radial direction (Nilsson et al. 2017). Particle-in-cell simulations have been used to correct for the spacecraft potential in ICA data by Bergman et al. (2021b), showing a surprising ‘backstream’ of ions towards the nucleus.

Despite these uncertainties, it is likely that the cometary ions within the diamagnetic cavity are indeed travelling faster than the neutral population, and are therefore collisionally uncoupled to some extent. Since the ionizing photons and electrons carry a negligible momentum compared to the neutrals, cometary ions are produced at the neutral speed, but are then susceptible to electromagnetic fields. Collisional decoupling from the neutral flow arises when the ions are accelerated beyond the neutral speed. Determination of the extent of this decoupling is important to understand how the cometary plasma behaves (including formation of the diamagnetic cavity boundary), but it is currently not well understood.

Outside the diamagnetic cavity and solar wind ion cavity, cometary ions are accelerated by the convective electric field of the solar wind, which leads to the mass loading and ion pick up process (Szegö et al. 2000; Behar et al. 2016). Inside the unmagnetized diamagnetic cavity region, the $\vec{J} \times \vec{B}$ term in Ohm’s law vanishes, and the ambipolar electric field becomes dominant. Unlike at planets, the electron pressure gradient does not arise from the different gravitational forces acting on ions and electrons (e.g. Schunk & Nagy 2009), but from the higher energy and pressure of the cometary electrons compared to the ions due to their smaller mass. An electric field arises to oppose the creation of a charge imbalance, in order to satisfy quasi-neutrality on scales longer than the Debye length. At comet 1P during the Giotto fly-by, Gan & Cravens (1990) demonstrated that the ambipolar field (and therefore the electron pressure) was negligible, since the dense neutral coma led to efficient collisional cooling of electrons. At 67P, for low outgassing activity, the ambipolar electric field has been shown to have a tangible impact on the cometary plasma environment. One such impact is keeping cometary electrons ‘trapped’ in a region close to the nucleus (e.g. Deca et al. 2019; Stephenson et al. 2022), thereby increasing the efficiency of electron-neutral collisions and leading to a cold (though minor) electron population (<0.1 eV) consistent with observations by RPC-LAP and RPC-MIP (e.g. Eriksson et al. 2017; Engelhardt et al. 2018; Gilet et al. 2020; Wattieaux et al. 2020). In addition, solar wind electrons were found to have been accelerated towards the nucleus by such an electric field, detected as a suprathermal electron population (Deca et al. 2017; Madanian et al. 2017) and responsible for generating aurora (Galand et al. 2020; Stephenson et al. 2021) and ionization (Stephenson et al. 2023).

The magnitude and functional form of the ambipolar electric field are difficult to constrain. Cravens (2004) derived an expression for the radial electric field ($\vec{E} = E(r) \vec{r}$) by first enforcing $n_e \approx n_i \propto 1/r$ while $T_e$ is assumed constant and isotropic throughout the coma. This arises from a field-free and chemistry-free model of the coma,
with a pure H₂O coma and the comet as a point source. Under these assumptions, the electric field becomes

$$E(r) = -\frac{1}{qn_e(r)} \frac{dp_e(r)}{dr} = -\frac{k_B T_e}{qn_e(r)} \frac{d(n_e(r))}{dr} = \frac{k_B T_e}{qr}, \quad (1)$$

suggesting a $\propto 1/r$ dependence. However, that may not be applicable between the surface and the radial location of the ionospheric peak located around $2r_c$ (or lower) as $n_i$ increases (Heritier et al. 2017a).

An argument can also be made for a $1/r^2$ dependence of the electric field by enforcing neutrality of the plasma. Gauss's law under spherical symmetry gives

$$\frac{1}{r^2} \frac{d(E(r)r^2)}{dr} = \frac{q(n_i - n_e)}{\epsilon_0}, \quad (2)$$

which, assuming strict neutrality ($n_i = n_e$), leads to

$$E(r) = E(r_c) \frac{r_c^2}{r^2}, \quad (3)$$

where $E(r_c)$ is the electric field strength at the comet surface. This is equivalent to the electric field generated by a charged sphere.

This study aims to further constrain the electric field with a 1D ionospheric model within the diamagnetic cavity including acceleration of the cometary ions by an electric field. The model is first described in Section 2, including validation against different approaches (see Section 2.4). In Section 3, the sensitivity of the ion composition and total plasma density to the electric field and momentum transfer is assessed. Finally in Section 4, we compare the total plasma density from the model to the measured electron density from RPC-MIP and instruments, which we use to constrain ion acceleration within the diamagnetic cavity.

## 2 MODEL DESCRIPTION

The 1D ionospheric model used in this paper is based on the one described in Heritier et al. (2017b), restricted to the ion species H₂O⁺, H₃O⁺, and NH₄⁺. We have, however, updated it substantially to include acceleration of the ions above the neutral speed $u_n$ by an ambipolar electric field, as well as momentum transfer collisions.<fnref n="1" /> The coupled continuity equation applied to each ion species produced at a given cometocentric distance $r$ are solved using a finite difference method for the three ion species for logarithmically spaced spherical shells from the comet nucleus at $r = 2$ km up to $r = 1 \times 10^3$ km (300 bins), until a steady state is reached.

The simplified neutral gas coma model is comprised of a small fraction of NH₃ ($f_{NH_3}$) and the rest is water ($f_{H_2O} = 1 - f_{NH_3}$). For input values of $f_{NH_3}$, outgassing Q, and neutral speed $u_n$, the neutral density $n_n(r)$ of each species n is calculated using the simplified Haser model (Haser 1957):

$$n_n(r) = \frac{Q f_n}{4\pi u_n r^2}. \quad (4)$$

### 2.1 Ion continuity equation

Consider an ion population (j, s), of species j born at a cometocentric distance $r_s$. In 1D and spherical symmetry, the number density $n_{j,s}(r)$ of this species, along a given radial line, is governed by the continuity

equation

$$\frac{\partial n_{j,s}}{\partial t} + \frac{1}{r^2} \frac{\partial}{\partial r} (r^2 n_{j,s} u_{j,s}) = P_{j,s} - R_{j,s} n_{j,s}, \quad (5)$$

where $\vec{u}_{j,s}(r)$ is the bulk velocity of the species (j, s) at r (see Section 2.3). All ions are produced initially at the neutral radial velocity, so at $r = r_s$, $\vec{u}_{j,s}(r) = \vec{u}_n$.

$P_{j,s}(r, t)$ is the production rate of the ion population (j, s) (in cm⁻³ s⁻<fnref n="1" />), which is comprised of the contributions from ionization, ion-neutral chemistry (Section 2.2), and momentum transfer (see Section 2.3):

$$P_{j,s} = P_{j}^{\text{ioni}} + P_{j}^{\text{chem}} + P_{j,s}^{\text{MT}}. \quad (6)$$

Equivalently, $R_{j,s}(r, t)$ is the loss frequency of the ion population (j, s) (in s⁻<fnref n="1" />), and comprises losses due to chemistry (see Section 2.2) and momentum transfer (see Section 2.3):

$$R_{j,s} = R_{j}^{\text{chem}} + R_{j,s}^{\text{MT}}. \quad (7)$$

### 2.2 Ionization and ion chemistry

The two main ionization sources for cometary ions are photoionization and electron impact. However, the importance of electron impact has been shown to be low at the location of Rosetta compared to photoionization near perihelion when the outgassing is high (Stephenson et al. 2023). It is neglected in the model over the full range of cometocentric distances. Large electron impact ionization frequencies can be driven by acceleration of solar wind electrons (such as 67P away from perihelion; Madanian et al. 2017) or by absorption of X-rays at very large outgassing rates (such as 1P). Neither of these mechanisms are significant in the diamagnetic cavity of 67P near perihelion, so it is reasonable to neglect ionization by electron impact. Photoionization of the most dominant neutral species, H₂O, is the most important ionization reaction and results in the production of H₂O⁺ at the same speed as the neutral coma ($u_n$).

The ionization production rate is given by

$$P_{j,s}^{\text{ioni}} = v_{n \to j}^{\text{ioni}} n_n, \quad (8)$$

where $v_{n \to j}^{\text{ioni}}$ is the ionization frequency: the number of neutral molecules n (here H₂O) per second that are photoionized to form ion species j (here H₂O⁺). This frequency is calculated using the photoionization cross-sections for each neutral species $\sigma_{n \to j}^{\text{ioni}}(\lambda)$, combined with the attenuated solar EUV spectrum, calculated from the Lamber–Beer law using the unattenuated solar flux from TIMED/SEE (Woods et al. 2005) adjusted to the heliocentric distance of the comet. A time shift is applied to correct for the difference in solar phase angle between the comet and the Earth (Galand et al. 2016). The profiles of the photoionization frequency with cometocentric distance, for the two case studies discussed in Section 4, are shown in Fig. 1, attesting greater attenuation by the neutral coma when the outgassing is higher in 2015 July compared to 2015 November.

Once ions have been created through ionization, they can then undergo chemical reactions with other neutral species. A list of all the reactions considered in the model is given in Appendix A. The rates depend on the neutral temperature, for which we use the profile from the adiabatic model in Heritier et al. (2017b).

The production rate of ion population (j, s) at the neutral speed $u_n$ due to ion-neutral chemistry is given by

$$P_{j}^{\text{chem}} = \sum_{n, j' \neq j} k_{j', n \to j}^{\text{IN}} n_n n_{j'}, \quad (9)$$

<footnote n="1">¹It should be noted that there is a typo in the caption of fig. 4 of Heritier et al. (2017b); the neutral number density used to calculate the photoionization rates was $3 \times 10^7$ cm⁻³, instead of $2 \times 10^6$ cm⁻³.</footnote>
with a pure H₂O coma and the comet as a point source. Under these assumptions, the electric field becomes

$$E(r) = -\frac{1}{qn_e(r)} \frac{dp_e(r)}{dr} = -\frac{k_B T_e}{qn_e(r)} \frac{d(n_e(r))}{dr} = \frac{k_B T_e}{qr}, \quad (1)$$

suggesting a $\propto 1/r$ dependence. However, that may not be applicable between the surface and the radial location of the ionospheric peak located around $2r_c$ (or lower) as $n_i$ increases (Heritier et al. 2017a).

An argument can also be made for a $1/r^2$ dependence of the electric field by enforcing neutrality of the plasma. Gauss's law under spherical symmetry gives

$$\frac{1}{r^2} \frac{d(E(r)r^2)}{dr} = \frac{q(n_i - n_e)}{\epsilon_0}, \quad (2)$$

which, assuming strict neutrality ($n_i = n_e$), leads to

$$E(r) = E(r_c) \frac{r_c^2}{r^2}, \quad (3)$$

where $E(r_c)$ is the electric field strength at the comet surface. This is equivalent to the electric field generated by a charged sphere.

This study aims to further constrain the electric field with a 1D ionospheric model within the diamagnetic cavity including acceleration of the cometary ions by an electric field. The model is first described in Section 2, including validation against different approaches (see Section 2.4). In Section 3, the sensitivity of the ion composition and total plasma density to the electric field and momentum transfer is assessed. Finally in Section 4, we compare the total plasma density from the model to the measured electron density from RPC-MIP and instruments, which we use to constrain ion acceleration within the diamagnetic cavity.

## 2 MODEL DESCRIPTION

The 1D ionospheric model used in this paper is based on the one described in Heritier et al. (2017b), restricted to the ion species H₂O⁺, H₃O⁺, and NH₄⁺. We have, however, updated it substantially to include acceleration of the ions above the neutral speed $u_n$ by an ambipolar electric field, as well as momentum transfer collisions.<fnref n="1" /> The coupled continuity equation applied to each ion species produced at a given cometocentric distance $r$ are solved using a finite difference method for the three ion species for logarithmically spaced spherical shells from the comet nucleus at $r = 2$ km up to $r = 1 \times 10^3$ km (300 bins), until a steady state is reached.

The simplified neutral gas coma model is comprised of a small fraction of NH₃ ($f_{NH_3}$) and the rest is water ($f_{H_2O} = 1 - f_{NH_3}$). For input values of $f_{NH_3}$, outgassing $Q$, and neutral speed $u_n$, the neutral density $n_n(r)$ of each species $n$ is calculated using the simplified Haser model (Haser 1957):

$$n_n(r) = \frac{Q f_n}{4 \pi u_n r^2}. \quad (4)$$

### 2.1 Ion continuity equation

Consider an ion population ($j, s$), of species $j$ born at a cometocentric distance $r_s$. In 1D and spherical symmetry, the number density $n_{j,s}(r)$ of this species, along a given radial line, is governed by the continuity

equation

$$\frac{\partial n_{j,s}}{\partial t} + \frac{1}{r^2} \frac{\partial}{\partial r} (r^2 n_{j,s} u_{j,s}) = P_{j,s} - R_{j,s} n_{j,s}, \quad (5)$$

where $\vec{u}_{j,s}(r)$ is the bulk velocity of the species ($j, s$) at $r$ (see Section 2.3). All ions are produced initially at the neutral radial velocity, so at $r = r_s$, $\vec{u}_{j,s}(r) = \vec{u}_n$.

$P_{j,s}(r, t)$ is the production rate of the ion population ($j, s$) (in cm⁻³ s⁻<fnref n="1" />), which is comprised of the contributions from ionization, ion-neutral chemistry (Section 2.2), and momentum transfer (see Section 2.3):

$$P_{j,s} = P_{j}^{\text{ioni}} + P_{j}^{\text{chem}} + P_{j,s}^{\text{MT}}. \quad (6)$$

Equivalently, $R_{j,s}(r, t)$ is the loss frequency of the ion population ($j, s$) (in s⁻<fnref n="1" />), and comprises losses due to chemistry (see Section 2.2) and momentum transfer (see Section 2.3):

$$R_{j,s} = R_{j}^{\text{chem}} + R_{j,s}^{\text{MT}}. \quad (7)$$

### 2.2 Ionization and ion chemistry

The two main ionization sources for cometary ions are photoionization and electron impact. However, the importance of electron impact has been shown to be low at the location of Rosetta compared to photoionization near perihelion when the outgassing is high (Stephenson et al. 2023). It is neglected in the model over the full range of cometocentric distances. Large electron impact ionization frequencies can be driven by acceleration of solar wind electrons (such as 67P away from perihelion; Madanian et al. 2017) or by absorption of X-rays at very large outgassing rates (such as 1P). Neither of these mechanisms are significant in the diamagnetic cavity of 67P near perihelion, so it is reasonable to neglect ionization by electron impact. Photoionization of the most dominant neutral species, H₂O, is the most important ionization reaction and results in the production of H₂O⁺ at the same speed as the neutral coma ($u_n$).

The ionization production rate is given by

$$P_{j,s}^{\text{ioni}} = v_{n \to j}^{\text{ioni}} n_n, \quad (8)$$

where $v_{n \to j}^{\text{ioni}}$ is the ionization frequency: the number of neutral molecules $n$ (here H₂O) per second that are photoionized to form ion species $j$ (here H₂O⁺). This frequency is calculated using the photoionization cross-sections for each neutral species $\sigma_{n \to j}^{\text{ioni}}(\lambda)$, combined with the attenuated solar EUV spectrum, calculated from the Lamber–Beer law using the unattenuated solar flux from TIMED/SEE (Woods et al. 2005) adjusted to the heliocentric distance of the comet. A time shift is applied to correct for the difference in solar phase angle between the comet and the Earth (Galand et al. 2016). The profiles of the photoionization frequency with cometocentric distance, for the two case studies discussed in Section 4, are shown in Fig. 1, attesting greater attenuation by the neutral coma when the outgassing is higher in 2015 July compared to 2015 November.

Once ions have been created through ionization, they can then undergo chemical reactions with other neutral species. A list of all the reactions considered in the model is given in Appendix A. The rates depend on the neutral temperature, for which we use the profile from the adiabatic model in Heritier et al. (2017b).

The production rate of ion population ($j, s$) at the neutral speed $u_n$ due to ion-neutral chemistry is given by

$$P_{j}^{\text{chem}} = \sum_{n, j' \neq j} k_{j', n \to j}^{\text{IN}} n_n n_{j'}, \quad (9)$$

<footnote n="1">¹It should be noted that there is a typo in the caption of fig. 4 of Heritier et al. (2017b); the neutral number density used to calculate the photoionization rates was $3 \times 10^7$ cm⁻³, instead of $2 \times 10^6$ cm⁻³.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Photoionization rate profiles for the two cases given in Section 4 (see inputs in Table 1).</figcaption></figure>

where $k_{j',n \to j}^{\text{IN}}$ is the reaction rate coefficient for ion-neutral chemical reactions between the ion species $j'$ and the neutral species $n$, leading to the production of ion species $j$ (here $\text{H}_2\text{O}^+$ or $\text{NH}_4^+$). We exclude symmetric charge exchange. Since the production of one ion species through ion-neutral chemistry results in the loss of another, there is an equivalent term in the loss rate ($R_{j,s}$).

The other way that an ion of a particular species can be chemically lost is through ion-electron dissociative recombination (often abbreviated to DR herein). A list of all the DR reactions considered in the model is given in Appendix A. The recombination rate coefficient $\alpha_j^{\text{DR}}$ [cm$^3$ s$^{-1}$] is dependent on the electron temperature, assumed to be $10^4$ K (using average combined temperature in the cavity from the RPC-MIP data set; Wattiaux et al. 2020). This assumption of constant electron temperature through the coma is unrealistic, but the sensitivity of our modelled electron density to the DR rate is limited (Beth, Galand & Heritier 2019).

Putting together the loss frequency for ion species $j$ through both ion-neutral chemistry and DR, it follows that:

$$R_j^{\text{chem}} = \sum_{\substack{n, j' \neq j \\ \text{ion-neutral chemistry}}} k_{j,n \to j}^{\text{IN}} n_n + \alpha_j^{\text{DR}} n_e . \quad (10)$$

## 2.3 Ion acceleration by an electric field in presence of momentum transfer

The ion populations in the model are distinguished by their source location $r_s$ as well as their species $j$. The modelled ions are born at the neutral speed and subsequently accelerated by an electric field as they travel radially outwards from the nucleus. This means that at a given shell $r_q$, the ion density is the sum of the populations from each source below and equal to it ($r_s \le r_q$), which are all travelling with different bulk velocities.

The velocity profiles $u_{j,s}(r_q)$ are calculated before solving the continuity equations. We first define the electric field as $E(r_q) = E_c(r_c/r_q)^m$, where coefficient $m \in \{0, 1, 2\}$ is an integer and $E_c$ [V m$^{-1}$] is the electric field at the comet surface ($r_c$). Conservation

of energy then dictates that

$$\frac{1}{2} m_j u_{j,s}(r_q)^2 = \frac{1}{2} m_j u_n^2 + \int_{r_s}^{r_q} E_{\text{amb}}(r_c) \left(\frac{r_c}{r}\right)^m dr, \quad (11)$$

where $u_{j,s}(r_q)$ is the radial speed of the ions at $r_q$ that were produced at $r_s$, either through ionization, ion-neutral chemistry, or momentum transfer, and subsequently accelerated from the neutral speed $u_n$ by the electric field. For example, for an electric field $E \propto r^{-1}$, this then leads to an ion velocity profile:

$$u_{j,s}(r) = \sqrt{u_n^2 + \frac{2E_c r_c}{m_j} \ln \frac{r_q}{r_s}} \quad (12)$$

The velocity profile is calculated for the ion population of species $j$ produced at each shell $r_s$. The continuity equation is then solved for each population ($j, s$), and the production term $P_{j,s}$ is only non-zero in the source shells. The contributions from each source are then summed at each cometocentric distance shell to find the total density of each species:

$$n_j(r_q) = \sum_s n_{j,s}(r_q), \quad (13)$$

where $s$ varies to convert cometocentric distances from $r_s = r_c$ to $r_s = r_q$.

The acceleration of ions can then be interrupted by ion-neutral elastic collisions, where there is no change in ion species involved, but the momentum of the fast ion species is transferred to the slow neutral. We assume that this collision is completely elastic, leaving the previously fast ion now at the neutral speed and the neutral at the accelerated ion speed. The case we consider is then an upper limit for the contribution of momentum transfer that said the lower limit of the ion bulk velocity. We explore the sensitivity in Section 3.

The contribution of momentum transfer to the production rate of the ion population ($j, s$) is

$$P_{j,s}^{\text{MT}}(r_s) = \sum_{s' \le s} F_{j,s'}^{u_{j,s'} \to u_n}(r_q) n_{j,s'}(r_s), \quad (14)$$

where $F_{j,s'}^{u_{j,s'} \to u_n}$ [s$^{-1}$] is the collision frequency for ions in population ($j, s'$), where $s' < s$. Put simply, this is the rate at which ions slowed down have been 'left behind' by their original population produced at $r_{s'}$ and will then contribute to the population ($j, s$).

Conversely, the ions originally in population ($j, s$) and 'left behind' by the accelerating ions after undergoing an ion-neutral collision at radial distance $r_q$ will then make up the momentum transfer loss frequency. They are produced at neutral velocity $u_n$. In this model, we use simple constant collision rate coefficients $k_{j,n}^{\text{MT}}$ [cm$^3$s$^{-1}$] (given in Appendix A) such that

$$R_{j,s}^{\text{MT}}(r_q) = F_{j,s}^{u_{j,s} \to u_n}(r_q) = \sum_n k_{j,n}^{\text{MT}} n_n(r_q). \quad (15)$$

The bulk velocity of each species can then be calculated by summing over the sources

$$u_j(r_q) = \frac{\sum_s n_{j,s}(r_q) u_{j,s}(r_q)}{\sum_s n_{j,s}(r_q)} \quad (16)$$

and the total ion bulk velocity $u_i$ is the weighted mean of all species $u_i = (\sum_j n_j u_j) / \sum_j n_j$.

## 2.4 Validation of the model

In this section, the model described in Sections 2.1–2.3 is validated against analytical solutions for each of the three core components
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Photoionization rate profiles for the two cases given in Section 4 (see inputs in Table 1).</figcaption></figure>

where $k_{j',n \to j}^{\text{IN}}$ is the reaction rate coefficient for ion-neutral chemical reactions between the ion species $j'$ and the neutral species $n$, leading to the production of ion species $j$ (here $\text{H}_2\text{O}^+$ or $\text{NH}_4^+$). We exclude symmetric charge exchange. Since the production of one ion species through ion-neutral chemistry results in the loss of another, there is an equivalent term in the loss rate ($R_{j,s}$).

The other way that an ion of a particular species can be chemically lost is through ion-electron dissociative recombination (often abbreviated to DR herein). A list of all the DR reactions considered in the model is given in Appendix A. The recombination rate coefficient $\alpha_j^{\text{DR}}$ [cm$^3$ s$^{-1}$] is dependent on the electron temperature, assumed to be $10^4$ K (using average combined temperature in the cavity from the RPC-MIP data set; Wattiaux et al. 2020). This assumption of constant electron temperature through the coma is unrealistic, but the sensitivity of our modelled electron density to the DR rate is limited (Beth, Galand & Heritier 2019).

Putting together the loss frequency for ion species $j$ through both ion-neutral chemistry and DR, it follows that:

$$R_j^{\text{chem}} = \sum_{\substack{n, j' \neq j \\ \text{ion-neutral chemistry}}} k_{j,n \to j}^{\text{IN}} n_n + \alpha_j^{\text{DR}} n_e . \quad (10)$$

## 2.3 Ion acceleration by an electric field in presence of momentum transfer

The ion populations in the model are distinguished by their source location $r_s$ as well as their species $j$. The modelled ions are born at the neutral speed and subsequently accelerated by an electric field as they travel radially outwards from the nucleus. This means that at a given shell $r_q$, the ion density is the sum of the populations from each source below and equal to it ($r_s \le r_q$), which are all travelling with different bulk velocities.

The velocity profiles $u_{j,s}(r_q)$ are calculated before solving the continuity equations. We first define the electric field as $E(r_q) = E_c(r_c/r_q)^m$, where coefficient $m \in \{0, 1, 2\}$ is an integer and $E_c$ [V m$^{-1}$] is the electric field at the comet surface ($r_c$). Conservation

of energy then dictates that

$$\frac{1}{2} m_j u_{j,s}(r_q)^2 = \frac{1}{2} m_j u_n^2 + \int_{r_s}^{r_q} E_{\text{amb}}(r_c) \left(\frac{r_c}{r}\right)^m dr, \quad (11)$$

where $u_{j,s}(r_q)$ is the radial speed of the ions at $r_q$ that were produced at $r_s$, either through ionization, ion-neutral chemistry, or momentum transfer, and subsequently accelerated from the neutral speed $u_n$ by the electric field. For example, for an electric field $E \propto r^{-1}$, this then leads to an ion velocity profile:

$$u_{j,s}(r) = \sqrt{u_n^2 + \frac{2E_c r_c}{m_j} \ln \frac{r_q}{r_s}} \quad (12)$$

The velocity profile is calculated for the ion population of species $j$ produced at each shell $r_s$. The continuity equation is then solved for each population ($j, s$), and the production term $P_{j,s}$ is only non-zero in the source shells. The contributions from each source are then summed at each cometocentric distance shell to find the total density of each species:

$$n_j(r_q) = \sum_s n_{j,s}(r_q), \quad (13)$$

where $s$ varies to convert cometocentric distances from $r_s = r_c$ to $r_s = r_q$.

The acceleration of ions can then be interrupted by ion-neutral elastic collisions, where there is no change in ion species involved, but the momentum of the fast ion species is transferred to the slow neutral. We assume that this collision is completely elastic, leaving the previously fast ion now at the neutral speed and the neutral at the accelerated ion speed. The case we consider is then an upper limit for the contribution of momentum transfer that said the lower limit of the ion bulk velocity. We explore the sensitivity in Section 3.

The contribution of momentum transfer to the production rate of the ion population ($j, s$) is

$$P_{j,s}^{\text{MT}}(r_s) = \sum_{s' \le s} F_{j,s'}^{u_{j,s'} \to u_n}(r_q) n_{j,s'}(r_s), \quad (14)$$

where $F_{j,s'}^{u_{j,s'} \to u_n}$ [s$^{-1}$] is the collision frequency for ions in population ($j, s'$), where $s' < s$. Put simply, this is the rate at which ions slowed down have been 'left behind' by their original population produced at $r_{s'}$ and will then contribute to the population ($j, s$).

Conversely, the ions originally in population ($j, s$) and 'left behind' by the accelerating ions after undergoing an ion-neutral collision at radial distance $r_q$ will then make up the momentum transfer loss frequency. They are produced at neutral velocity $u_n$. In this model, we use simple constant collision rate coefficients $k_{j,n}^{\text{MT}}$ [cm$^3$s$^{-1}$] (given in Appendix A) such that

$$R_{j,s}^{\text{MT}}(r_q) = F_{j,s}^{u_{j,s} \to u_n}(r_q) = \sum_n k_{j,n}^{\text{MT}} n_n(r_q). \quad (15)$$

The bulk velocity of each species can then be calculated by summing over the sources

$$u_j(r_q) = \frac{\sum_s n_{j,s}(r_q) u_{j,s}(r_q)}{\sum_s n_{j,s}(r_q)} \quad (16)$$

and the total ion bulk velocity $u_i$ is the weighted mean of all species $u_i = (\sum_j n_j u_j) / \sum_j n_j$.

## 2.4 Validation of the model

In this section, the model described in Sections 2.1–2.3 is validated against analytical solutions for each of the three core components
of the numerical scheme: momentum transfer (Model I), ion-neutral chemistry (Model II), and electric field acceleration (Model III). In all cases, dissociative recombination is neglected, the neutral speed $u_n = 1$ km s⁻¹, the photoionization frequency $1 \times 10^7$ s⁻¹, and comet radius $r_c = 2$ km.

### 2.4.1 Analytical model I: momentum transfer

First, we check the momentum transfer numerical scheme (see Section 2.3) against a simple analytical model where we consider a population of ions (single species) produced from the surface at a speed $u_i > u_n$. The ions are then transported radially outwards, undergoing momentum transfer collisions with the neutrals, but with no ion-neutral chemistry or dissociative recombination loss. Ions that have collided with a neutral have their speed reduced to the neutral speed. We then have two populations of ions, ‘hot’ ($n_H(r)$) and ‘cold’ ($n_C(r)$) at constant speeds $u_i$ and $u_n$, respectively. The two populations follow the coupled continuity equations:

$$ \frac{1}{r^2} \frac{d}{dr} (n_H u_i r^2) = -k^{\text{MT}} n_H n_n \quad (17) $$

$$ \frac{1}{r^2} \frac{d}{dr} (n_C u_n r^2) = k^{\text{MT}} n_H n_n, \quad (18) $$

which can be solved to find

$$ n_C = \frac{n_H(r_c) r_c^2 u_i}{r^2 u_n} \left( 1 - \exp \left( \frac{C_0}{r} - \frac{C_0}{r_c} \right) \right) \quad (19) $$

$$ n_H = \frac{n_H(r_c) r_c^2}{r^2} \exp \left( \frac{C_0}{r} - \frac{C_0}{r_c} \right), \quad (20) $$

where $C_0 = \frac{k^{\text{MT}} Q}{4\pi u_n u_i}$. The results of comparison between equations (19) and (20), and the numerical model with the same assumptions are given in Fig. 2a, for $Q = 10^{27}$ s⁻¹ and $k^{\text{MT}} = 5 \times 10^{-11}$ cm³ s⁻¹. The numerical and analytical models show a very good agreement, with the maximum relative difference in the total ion density of 0.5 per cent for the logarithmically spaced bins described in Section 2, and for $u_n = 1$ km s⁻¹ and $u_i = 3$ km s⁻¹. The difference results from the finite cell size in the numerical approach.

### 2.4.2 Analytical model II: ion-neutral chemistry

Next, we validate the implementation of ion-neutral chemistry in the model. With acceleration, dissociative recombination, and momentum transfer neglected, and considering a monoenergetic solar radiation at noon, it is possible to derive an analytical solution to the coupled continuity equations for three ion species ($H_2O^+$, $H_3O^+$, and $NH_4^+$):

$$ \frac{1}{r^2} \frac{dn_{H_2O^+} u_n r^2}{dr} = \frac{\nu Q}{4\pi u_n r^2} \exp(-\tau r_c/r) - (k_1 f_{H_2O} + k_2 f_{NH_3}) \frac{Q}{4\pi u_n r^2} n_{H_2O^+}, \quad (21) $$

$$ \frac{1}{r^2} \frac{dn_{H_3O^+} u_n r^2}{dr} = (k_1 f_{H_2O}) \frac{Q}{4\pi u_n r^2} n_{H_2O^+} - (k_3 f_{NH_3}) \frac{Q}{4\pi u_n r^2} n_{H_3O^+}, \quad (22) $$

$$ \frac{1}{r^2} \frac{dn_{NH_4^+} u_n r^2}{dr} = (k_2 f_{NH_3}) \frac{Q}{4\pi u_n r^2} n_{H_2O^+} + (k_3 f_{NH_3}) \frac{Q}{4\pi u_n r^2} n_{H_3O^+}, \quad (23) $$

where $k_1, k_2$, and $k_3$ are the relevant chemical rate coefficients (see Section 2.2 and Appendix A) and $\tau$ is the optical depth. The solution (derived in Appendix B) is then given by

$$ n_i = \frac{\nu Q}{4\pi u_n^2 r_c} \frac{r_c^2}{r^2} \left[ \frac{r}{r_c} E_2 \left( \tau \frac{r_c}{r} \right) - E_2(\tau) \right] \quad (24) $$

$$ n_{H_2O^+} = \frac{\nu Q}{4\pi u_n^2 r_c} \frac{r_c^2}{r^2} \left[ \frac{r}{r_c} E_2 \left( (K_1 + K_2 + \tau) \frac{r_c}{r} \right) - E_2(K_1 + K_2 + \tau) \right] \quad (25) $$

$$ \times \exp \left[ (K_1 + K_2) \frac{r_c}{r} \right] \quad (26) $$

$$ n_{H_3O^+} = \frac{\nu Q}{4\pi u_n^2 r_c} \frac{r_c^2}{r^2} \frac{K_1}{K_1 + (K_2 - K_3)} \times \left( \left[ \frac{r}{r_c} E_2 \left( (K_1 + K_2 + \tau) \frac{r_c}{r} \right) - E_2(K_1 + K_2 + \tau) \right] \exp \left( (K_1 + K_2) \frac{r_c}{r} \right) - \left[ \frac{r}{r_c} E_2 \left( (K_3 + \tau) \frac{r_c}{r} \right) - E_2(K_3 + \tau) \right] \exp \left( (K_3 \frac{r_c}{r}) \right) \right) \quad (27) $$

$$ n_{NH_4^+} = n_i - n_{H_2O^+} - n_{H_3O^+}, \quad (28) $$

where $E_2(x) = \int_1^\infty \frac{e^{-xt}}{t^2} dt$ is the exponential integral function and

$$ K_1 = \frac{k_1 f_{H_2O} Q}{4\pi u_n^2 r_c} \quad (29) $$

$$ K_2 = \frac{k_2 f_{NH_3} Q}{4\pi u_n^2 r_c} \quad (30) $$

$$ K_3 = \frac{k_3 f_{NH_3} Q}{4\pi u_n^2 r_c} \quad (31) $$

$$ \tau = \frac{\sigma Q}{4\pi u_n r_c}. \quad (32) $$

The solutions are overplotted in Fig. 2b, with the result of the numerical model with the same assumptions and for $Q = 10^{27}$ s⁻¹, $f_{NH_3} = 0.01$, $f_{H_2O} = 0.99$, and $\tau = 0$ (optically thin coma at the surface). The two models show again very good agreement, with maximum variation 5 per cent, for again the logarithmically spaced bins described in Section 2.

### 2.4.3 Analytical model III: electric field acceleration

Finally, the effect on the number density of the electric field can be validated analytically by the consideration of a single species model with no momentum transfer or ion-neutral chemistry. The model is described in Appendix C and is compared to the present numerical model in Fig. 3, for electric fields $\propto \text{const.}$ (a), $\propto r^{-1}$ (b), and $\propto r^{-2}$ (c). Again, the maximum difference between the analytical model and our numerical scheme is around 5 per cent.

## 3 ION COMPOSITION AND IONOSPHERIC DENSITIES

In this section, we explore the sensitivity of the model output to momentum transfer and to the electric field.

First, we examine the effect of the electric field on each ion species. Fig. 4 (left) shows the ion density profiles for input electric field profiles $\propto 1/r$ with three different surface field strengths: 0, 1, and 10 mV m⁻¹. The time-scales for each ion loss process are also shown in Fig. 4, right panel. The density profile of $H_2O^+$ varies very little with the electric field strength (Fig. 4a), owing to the proton transfer time-scale that is always shorter than both the transport and dissociative recombination time-scales (Fig. 4b). As
of the numerical scheme: momentum transfer (Model I), ion-neutral chemistry (Model II), and electric field acceleration (Model III). In all cases, dissociative recombination is neglected, the neutral speed $u_n = 1$ km s⁻¹, the photoionization frequency $1 \times 10^7$ s⁻¹, and comet radius $r_c = 2$ km.

### 2.4.1 Analytical model I: momentum transfer

First, we check the momentum transfer numerical scheme (see Section 2.3) against a simple analytical model where we consider a population of ions (single species) produced from the surface at a speed $u_i > u_n$. The ions are then transported radially outwards, undergoing momentum transfer collisions with the neutrals, but with no ion-neutral chemistry or dissociative recombination loss. Ions that have collided with a neutral have their speed reduced to the neutral speed. We then have two populations of ions, ‘hot’ ($n_H(r)$) and ‘cold’ ($n_C(r)$) at constant speeds $u_i$ and $u_n$, respectively. The two populations follow the coupled continuity equations:

$$ \frac{1}{r^2} \frac{d}{dr} (n_H u_i r^2) = -k^{\text{MT}} n_H n_n \quad (17) $$

$$ \frac{1}{r^2} \frac{d}{dr} (n_C u_n r^2) = k^{\text{MT}} n_H n_n, \quad (18) $$

which can be solved to find

$$ n_C = \frac{n_H(r_c) r_c^2 u_i}{r^2 u_n} \left( 1 - \exp \left( \frac{C_0}{r} - \frac{C_0}{r_c} \right) \right) \quad (19) $$

$$ n_H = \frac{n_H(r_c) r_c^2}{r^2} \exp \left( \frac{C_0}{r} - \frac{C_0}{r_c} \right), \quad (20) $$

where $C_0 = \frac{k^{\text{MT}} Q}{4\pi u_n u_i}$. The results of comparison between equations (19) and (20), and the numerical model with the same assumptions are given in Fig. 2a, for $Q = 10^{27}$ s⁻¹ and $k^{\text{MT}} = 5 \times 10^{-11}$ cm³ s⁻¹. The numerical and analytical models show a very good agreement, with the maximum relative difference in the total ion density of 0.5 per cent for the logarithmically spaced bins described in Section 2, and for $u_n = 1$ km s⁻¹ and $u_i = 3$ km s⁻¹. The difference results from the finite cell size in the numerical approach.

### 2.4.2 Analytical model II: ion-neutral chemistry

Next, we validate the implementation of ion-neutral chemistry in the model. With acceleration, dissociative recombination, and momentum transfer neglected, and considering a monoenergetic solar radiation at noon, it is possible to derive an analytical solution to the coupled continuity equations for three ion species ($H_2O^+$, $H_3O^+$, and $NH_4^+$):

$$ \frac{1}{r^2} \frac{dn_{H_2O^+} u_n r^2}{dr} = \frac{\nu Q}{4\pi u_n r^2} \exp(-\tau r_c/r) - (k_1 f_{H_2O} + k_2 f_{NH_3}) \frac{Q}{4\pi u_n r^2} n_{H_2O^+}, \quad (21) $$

$$ \frac{1}{r^2} \frac{dn_{H_3O^+} u_n r^2}{dr} = (k_1 f_{H_2O}) \frac{Q}{4\pi u_n r^2} n_{H_2O^+} - (k_3 f_{NH_3}) \frac{Q}{4\pi u_n r^2} n_{H_3O^+}, \quad (22) $$

$$ \frac{1}{r^2} \frac{dn_{NH_4^+} u_n r^2}{dr} = (k_2 f_{NH_3}) \frac{Q}{4\pi u_n r^2} n_{H_2O^+} + (k_3 f_{NH_3}) \frac{Q}{4\pi u_n r^2} n_{H_3O^+}, \quad (23) $$

where $k_1, k_2$, and $k_3$ are the relevant chemical rate coefficients (see Section 2.2 and Appendix A) and $\tau$ is the optical depth. The solution (derived in Appendix B) is then given by

$$ n_i = \frac{\nu Q}{4\pi u_n^2 r_c} \frac{r_c^2}{r^2} \left[ \frac{r}{r_c} E_2 \left( \tau \frac{r_c}{r} \right) - E_2(\tau) \right] \quad (24) $$

$$ n_{H_2O^+} = \frac{\nu Q}{4\pi u_n^2 r_c} \frac{r_c^2}{r^2} \left[ \frac{r}{r_c} E_2 \left( (K_1 + K_2 + \tau) \frac{r_c}{r} \right) - E_2(K_1 + K_2 + \tau) \right] \quad (25) $$

$$ \times \exp \left[ (K_1 + K_2) \frac{r_c}{r} \right] \quad (26) $$

$$ n_{H_3O^+} = \frac{\nu Q}{4\pi u_n^2 r_c} \frac{r_c^2}{r^2} \frac{K_1}{K_1 + (K_2 - K_3)} \times \left( \left[ \frac{r}{r_c} E_2 \left( (K_1 + K_2 + \tau) \frac{r_c}{r} \right) - E_2(K_1 + K_2 + \tau) \right] \exp \left( (K_1 + K_2) \frac{r_c}{r} \right) \right. \\ \left. - \left[ \frac{r}{r_c} E_2 \left( (K_3 + \tau) \frac{r_c}{r} \right) - E_2(K_3 + \tau) \right] \exp \left( K_3 \frac{r_c}{r} \right) \right) \quad (27) $$

$$ n_{NH_4^+} = n_i - n_{H_2O^+} - n_{H_3O^+}, \quad (28) $$

where $E_2(x) = \int_1^\infty \frac{e^{-xt}}{t^2} dt$ is the exponential integral function and

$$ K_1 = \frac{k_1 f_{H_2O} Q}{4\pi u_n^2 r_c} \quad (29) $$

$$ K_2 = \frac{k_2 f_{NH_3} Q}{4\pi u_n^2 r_c} \quad (30) $$

$$ K_3 = \frac{k_3 f_{NH_3} Q}{4\pi u_n^2 r_c} \quad (31) $$

$$ \tau = \frac{\sigma Q}{4\pi u_n r_c}. \quad (32) $$

The solutions are overplotted in Fig. 2b, with the result of the numerical model with the same assumptions and for $Q = 10^{27}$ s⁻¹, $f_{NH_3} = 0.01$, $f_{H_2O} = 0.99$, and $\tau = 0$ (optically thin coma at the surface). The two models show again very good agreement, with maximum variation 5 per cent, for again the logarithmically spaced bins described in Section 2.

### 2.4.3 Analytical model III: electric field acceleration

Finally, the effect on the number density of the electric field can be validated analytically by the consideration of a single species model with no momentum transfer or ion-neutral chemistry. The model is described in Appendix C and is compared to the present numerical model in Fig. 3, for electric fields $\propto$const. (a), $\propto r^{-1}$ (b), and $\propto r^{-2}$ (c). Again, the maximum difference between the analytical model and our numerical scheme is around 5 per cent.

## 3 ION COMPOSITION AND IONOSPHERIC DENSITIES

In this section, we explore the sensitivity of the model output to momentum transfer and to the electric field.

First, we examine the effect of the electric field on each ion species. Fig. 4 (left) shows the ion density profiles for input electric field profiles $\propto 1/r$ with three different surface field strengths: 0, 1, and 10 mV m⁻¹. The time-scales for each ion loss process are also shown in Fig. 4, right panel. The density profile of $H_2O^+$ varies very little with the electric field strength (Fig. 4a), owing to the proton transfer time-scale that is always shorter than both the transport and dissociative recombination time-scales (Fig. 4b). As
<figure><img src="image_3.png" /><figcaption>**Figure 2.** [a] Density profile of ions at $u_i$ and $u_n$ and the total density for the present numerical model in comparison with analytical model I described in equations (19) and (20). [b] Density profile of $H_2O^+$, $H_3O^+$, and $NH_4^+$, as well as the total density for the present numerical model in comparison with analytical model II given in B. In both [a] and [b], the coloured lines show the numerical model, and the corresponding analytical model is shown by black dashed lines.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Total ion density profiles from the numerical model (coloured lines) compared to analytical model III (black dashed lines) given in Appendix C for $Q = 10^{28} s^{-1}$. Electric fields are of varying strengths and $\propto$const. (a), $\propto r^{-1}$ (b), and $\propto r^{-2}$ (c).</figcaption></figure>

previously documented (e.g. Vigren & Galand 2013), $H_2O^+$ is close to photochemical equilibrium at this outgassing rate ($10^{28} s^{-1}$) and low cometocentric distances, such that $n_{H_2O^+} \sim \nu/k$ and mostly governed by the change in ionization frequency. At high cometocentric distance, the $10 \text{ mV m}^{-1}$ case begins to exhibit a small departure from photochemical equilibrium, as the transport time-scale is decreased such that it is more similar in magnitude to the proton transfer time-scale.

$H_3O^+$ is only produced through the proton transfer from $H_2O^+$ and is therefore more sensitive to the increasing electric field than $H_2O^+$ (see Figs 4 c and 4d). With the strongest electric field, the proton transfer has the shortest time-scale below 10 km, so loss of $H_3O^+$ to $NH_4^+$ dominates close to the nucleus. With no electric field, this continues up to 100 km. Transport then takes over (shorter time-scale), decreasing the ion density further from the surface. In the electric field case, the $H_3O^+$ density profile reaches a slope in $1/r$ above ~100 km. When the electric field is applied, the dissociative recombination time-scale is always at least an order of magnitude

larger than the transport time-scale, so this is not a significant loss process for this ion species.

$NH_4^+$ is the most sensitive of the three ion species to the increasing electric field, with the density over an order of magnitude lower at $100 \text{ km}$ with a $10 \text{ mV m}^{-1}$ electric field than when there is none. This means that the presence and strength of the ambipolar electric field may have important consequences for the detectability of $NH_4^+$ within the diamagnetic cavity. The transport time-scale is always shorter than the dissociative recombination time-scale. The difference is particularly marked when the strong ($10 \text{ mV m}^{-1}$) electric field is applied: transport dominated. This, combined with the lack of $NH_4^+$ production at high cometocentric distances means that it is quickly transported away, and the ion density slope is in $1/r^2$ above ~300 km. This is the same as the pure transport solution (no production or chemical loss) to the continuity equation (equation 5). With lower electric field, DR becomes increasingly important at high cometocentric distances and the slope of the $NH_4^+$ profile is between $r^{-1}$ and $r^{-2}$.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** [a] Density profile of ions at $u_i$ and $u_n$ and the total density for the present numerical model in comparison with analytical model I described in equations (19) and (20). [b] Density profile of $H_2O^+$, $H_3O^+$, and $NH_4^+$, as well as the total density for the present numerical model in comparison with analytical model II given in B. In both [a] and [b], the coloured lines show the numerical model, and the corresponding analytical model is shown by black dashed lines.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Total ion density profiles from the numerical model (coloured lines) compared to analytical model III (black dashed lines) given in Appendix C for $Q = 10^{28} \text{ s}^{-1}$. Electric fields are of varying strengths and $\propto \text{const.}$ (a), $\propto r^{-1}$ (b), and $\propto r^{-2}$ (c).</figcaption></figure>

previously documented (e.g. Vigren & Galand 2013), $H_2O^+$ is close to photochemical equilibrium at this outgassing rate ($10^{28} \text{ s}^{-1}$) and low cometocentric distances, such that $n_{H_2O^+} \sim \nu/k$ and mostly governed by the change in ionization frequency. At high cometocentric distance, the $10 \text{ mV m}^{-1}$ case begins to exhibit a small departure from photochemical equilibrium, as the transport time-scale is decreased such that it is more similar in magnitude to the proton transfer time-scale.

$H_3O^+$ is only produced through the proton transfer from $H_2O^+$ and is therefore more sensitive to the increasing electric field than $H_2O^+$ (see Figs 4 c and 4d). With the strongest electric field, the proton transfer has the shortest time-scale below 10 km, so loss of $H_3O^+$ to $NH_4^+$ dominates close to the nucleus. With no electric field, this continues up to 100 km. Transport then takes over (shorter time-scale), decreasing the ion density further from the surface. In the electric field case, the $H_3O^+$ density profile reaches a slope in $1/r$ above ~100 km. When the electric field is applied, the dissociative recombination time-scale is always at least an order of magnitude

larger than the transport time-scale, so this is not a significant loss process for this ion species.

$NH_4^+$ is the most sensitive of the three ion species to the increasing electric field, with the density over an order of magnitude lower at $100 \text{ km}$ with a $10 \text{ mV m}^{-1}$ electric field than when there is none. This means that the presence and strength of the ambipolar electric field may have important consequences for the detectability of $NH_4^+$ within the diamagnetic cavity. The transport time-scale is always shorter than the dissociative recombination time-scale. The difference is particularly marked when the strong ($10 \text{ mV m}^{-1}$) electric field is applied: transport dominated. This, combined with the lack of $NH_4^+$ production at high cometocentric distances means that it is quickly transported away, and the ion density slope is in $1/r^2$ above ~300 km. This is the same as the pure transport solution (no production or chemical loss) to the continuity equation (equation 5). With lower electric field, DR becomes increasingly important at high cometocentric distances and the slope of the $NH_4^+$ profile is between $r^{-1}$ and $r^{-2}$.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Ion density (left column) and loss time-scales (right column) for H₂O⁺ (top), H₃O⁺ (middle), and NH₄⁺ (bottom). Three different electric field conditions are considered: no electric field, *E*(*r*<sub>c</sub>) = 1 mV m⁻¹, and *E*(*r*<sub>c</sub>) = 10 mV m⁻¹. The three loss time-scales are transport (solid line), dissociative recombination (DR, dashed), and proton transfer (PT, dot-dashed). The electric field is assumed to be radial and proportional to 1/r, the outgassing *Q* = 10²⁸ s⁻¹, and the neutral composition to be 1 per cent NH₃ and 99 per cent H₂O.</figcaption></figure>

Fig. 5 shows the effect of including momentum transfer on the density of the three ion species as well as their total. As anticipated, the addition of momentum transfer increases the ion density, since the acceleration process is interrupted, slowing the bulk ion speed and allowing the ionosphere to build up more before it is transported. H₂O⁺, however, is unaffected since it is close to photochemical equilibrium and therefore not sensitive to the transport time-scale. Momentum transfer has the greatest impact on the NH₄⁺ density, not only because of its decreased transport loss, but also because of its enhanced production. The slowing of H₃O⁺ ions through ion-neutral collisions allows protonation to occur more readily.

The influence of both the electric field and momentum transfer on the total density and ion bulk velocity is summarized in Fig. 6, for electric field proportional to 1/r (a and b) and 1/r² (c and d). In all cases, the total density is higher when momentum transfer is included, and the bulk speed is lower. Vigren et al. (2015) found for the case of an electric field decreasing in r⁻¹, in the absence of momentum transfer and for the case of water ions, the *n*ᵢ ∝ 1/r relation is recovered above the density peak for all electric field strengths. We found that this remains valid for the case of three ions species with momentum transfer included.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Ion density (left column) and loss time-scales (right column) for H₂O⁺ (top), H₃O⁺ (middle), and NH₄⁺ (bottom). Three different electric field conditions are considered: no electric field, *E*(*r*<sub>c</sub>) = 1 mV m⁻¹, and *E*(*r*<sub>c</sub>) = 10 mV m⁻¹. The three loss time-scales are transport (solid line), dissociative recombination (DR, dashed), and proton transfer (PT, dot-dashed). The electric field is assumed to be radial and proportional to 1/r, the outgassing *Q* = 10²⁸ s⁻¹, and the neutral composition to be 1 per cent NH₃ and 99 per cent H₂O.</figcaption></figure>

Fig. 5 shows the effect of including momentum transfer on the density of the three ion species as well as their total. As anticipated, the addition of momentum transfer increases the ion density, since the acceleration process is interrupted, slowing the bulk ion speed and allowing the ionosphere to build up more before it is transported. H₂O⁺, however, is unaffected since it is close to photochemical equilibrium and therefore not sensitive to the transport time-scale. Momentum transfer has the greatest impact on the NH₄⁺ density, not only because of its decreased transport loss, but also because of its enhanced production. The slowing of H₃O⁺ ions through ion-neutral collisions allows protonation to occur more readily.

The influence of both the electric field and momentum transfer on the total density and ion bulk velocity is summarized in Fig. 6, for electric field proportional to 1/r (a and b) and 1/r² (c and d). In all cases, the total density is higher when momentum transfer is included, and the bulk speed is lower. Vigren et al. (2015) found for the case of an electric field decreasing in r⁻¹, in the absence of momentum transfer and for the case of water ions, the *n*ᵢ ∝ 1/r relation is recovered above the density peak for all electric field strengths. We found that this remains valid for the case of three ions species with momentum transfer included.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** [a] Ion density for H₂O⁺, H₃O⁺, and NH₄⁺, both with momentum transfer (solid line) and without (dashed line) momentum transfer. [b] Ion bulk velocity for each ion shown in [a], calculated according to equation (16). The electric field was assumed to be *E* = 1 × 10⁻³(*r*ᵣ/*r*) V m⁻¹, the outgassing *Q* = 10²⁸ s⁻¹, and the neutral composition to be 1 per cent NH₃ and 99 per cent H₂O.</figcaption></figure>

# 4 COMPARISON WITH ROSETTA DATA

We now focus on two key time periods during the *Rosetta* escort phase, comparing the total density derived from the RPC-MIP and instruments with the total ion density from our model for various electric field profiles. In doing so, we constrain the strength of the electric field which is required to explain the measurements. We also compare the predicted versus the measured NH₄⁺ density with the same inputs and electric field strengths.

The first period we consider is 2015 November 20–21, and a summary of the key data is shown in Fig. 7. *Rosetta* spent a significant amount of time inside the diamagnetic cavity (blue shaded vertical boxes) during this window, and there is strong coverage of the combined MIP/LAP data set (black, bottom panel). On November 21, the ROSINA/DFMS high-resolution ion mode was briefly active, and captured signatures of NH₄⁺ consistently (red, vertical lines, middle panel) (see Lewis et al. 2023), hence this makes an interesting case study for comparison with our model.

We also consider 2015 July 29–30, two weeks before perihelion, which saw some of the strongest NH₄⁺ signatures and many diamagnetic cavity crossings (see Fig. 8). The MIP/LAP combined data set was not available for this period, so we use electron density data from RPC-MIP only (red, bottom panel).

The distribution of electron density measurements that were taken while *Rosetta* was inside the diamagnetic cavity for each time period is shown in the histograms in Fig. 9. The data are similarly distributed, with less variation in November when more data were available. For comparison with the model, we take the range of electron densities one standard deviation on either side of the mean (indicated by the blue shaded regions), which was around twice as high in 2015 July (955.1–1517 cm⁻³) compared to 2015 November (544.4–791.8 cm⁻³).

A summary of the model inputs and range of electron densities for each period is shown in Table 1. The neutral speed was taken from Biver et al. (2019), and the outgassing *Q* was derived from this and the COPS neutral density (see equation 4). It is worth noting that the data availability of COPS was limited for the November period, so only a best estimate of the average outgassing over the whole period can be used. The fraction of NH₃ was estimated from the DFMS neutral mode data shown in Figs 7 and 8.

For a range of electric field surface strengths *E*(*r*ᵣ), the ionospheric model was run for the inputs in Table 1, for the time periods in 2015 November and July. Fig. 10 a shows the modelled total ion density at the cometocentric distance of *Rosetta* for each electric field strength. The profiles are calculated for both *E* ∝ *r*⁻¹ (solid lines) and *E* ∝ *r*⁻² (dashed lines). The horizontal shaded regions in red and grey show the electron density data for 2015 July and November, respectively (derived in Fig. 9). The *r*⁻² profiles are clearly unable to explain the measured electron densities and ion speeds, since an unreasonably high electric field would be required to reduce the electron density to the measured range at the *Rosetta* location. Hence, we follow the *r*⁻¹ profiles to derive the surface electric field strength range that would explain the measured electron density for each time period. We find electric field strengths for July and November of around 1.1–2.7 and 0.2–0.6 mV m⁻¹, respectively. Stronger electric fields (>3 mV m⁻¹) lead to enhanced ion transport, reducing the plasma density below the observed range.

We now examine the predicted NH₄⁺ density for the same input conditions and range of electric field measurements as in the previous section. Fig. 10 b shows a decreasing trend of NH₄⁺ density with increased electric field strength, similar to the total plasma density. As discussed in Section 3, NH₄⁺ is more sensitive to the change in the electric field strength, and therefore varies over more orders of magnitude than the total ion density. For the total plasma densities measured in the two time periods, the model predicts that this would correspond to 20–40 cm⁻³ of NH₄⁺ on 2015 November 20–21 and 50–110 NH₄⁺ for 2015 July 29–30.

While we know NH₄⁺ was consistently detected in the diamagnetic cavity during the two case study periods while the mass spectrometer was in ion mode (see Figs 8 and 7), the number density of NH₄⁺ cannot be inferred from the ROSINA/DFMS spectra with identified NH₄⁺ peaks. It is likely that even the lowest modelled densities in Fig. 10 b would be detectable by the instrument in the high-resolution mode (see Lewis et al. 2023, Appendix A), and NH₄⁺ detection with DFMS is also dependent on the field of view and energy acceptance window of the instrument. Therefore, it is not possible to use the detection of NH₄⁺ to constrain the electric field in this way.

Fig. 10c shows the ion bulk velocity for the same model runs as Figs 10 a and 10b. For the lowest electric fields (≈1 × 10⁻⁵ Vm⁻¹), the ions are not significantly accelerated above the neutral speed. When the electric field is high (>1 × 10⁻² V m⁻¹), the ions are further accelerated, up to 10 km s⁻¹. Such high speeds have been derived from RPC-ICA measurements (Bergman et al. 2021a), but the present ionospheric modelling suggests that such strong acceleration would lead to electron densities around 6 times lower than what have
<figure><img src="image_6.png" /><figcaption>**Figure 5.** [a] Ion density for H₂O⁺, H₃O⁺, and NH₄⁺, both with momentum transfer (solid line) and without (dashed line) momentum transfer. [b] Ion bulk velocity for each ion shown in [a], calculated according to equation (16). The electric field was assumed to be *E* = 1 × 10⁻³(*r*<sub>c</sub>/r) V m⁻¹, the outgassing *Q* = 10²⁸ s⁻¹, and the neutral composition to be 1 per cent NH₃ and 99 per cent H₂O.</figcaption></figure>

## 4 COMPARISON WITH ROSETTA DATA

We now focus on two key time periods during the *Rosetta* escort phase, comparing the total density derived from the RPC-MIP and instruments with the total ion density from our model for various electric field profiles. In doing so, we constrain the strength of the electric field which is required to explain the measurements. We also compare the predicted versus the measured NH₄⁺ density with the same inputs and electric field strengths.

The first period we consider is 2015 November 20–21, and a summary of the key data is shown in Fig. 7. *Rosetta* spent a significant amount of time inside the diamagnetic cavity (blue shaded vertical boxes) during this window, and there is strong coverage of the combined MIP/LAP data set (black, bottom panel). On November 21, the ROSINA/DFMS high-resolution ion mode was briefly active, and captured signatures of NH₄⁺ consistently (red, vertical lines, middle panel) (see Lewis et al. 2023), hence this makes an interesting case study for comparison with our model.

We also consider 2015 July 29–30, two weeks before perihelion, which saw some of the strongest NH₄⁺ signatures and many diamagnetic cavity crossings (see Fig. 8). The MIP/LAP combined data set was not available for this period, so we use electron density data from RPC-MIP only (red, bottom panel).

The distribution of electron density measurements that were taken while *Rosetta* was inside the diamagnetic cavity for each time period is shown in the histograms in Fig. 9. The data are similarly distributed, with less variation in November when more data were available. For comparison with the model, we take the range of electron densities one standard deviation on either side of the mean (indicated by the blue shaded regions), which was around twice as high in 2015 July (955.1–1517 cm⁻³) compared to 2015 November (544.4–791.8 cm⁻³).

A summary of the model inputs and range of electron densities for each period is shown in Table 1. The neutral speed was taken from Biver et al. (2019), and the outgassing *Q* was derived from this and the COPS neutral density (see equation 4). It is worth noting that the data availability of COPS was limited for the November period, so only a best estimate of the average outgassing over the whole period can be used. The fraction of NH₃ was estimated from the DFMS neutral mode data shown in Figs 7 and 8.

For a range of electric field surface strengths *E*(*r*<sub>c</sub>), the ionospheric model was run for the inputs in Table 1, for the time periods in 2015 November and July. Fig. 10 a shows the modelled total ion density at the cometocentric distance of *Rosetta* for each electric field strength. The profiles are calculated for both *E* ∝ *r*⁻¹ (solid lines) and *E* ∝ *r*⁻² (dashed lines). The horizontal shaded regions in red and grey show the electron density data for 2015 July and November, respectively (derived in Fig. 9). The *r*⁻² profiles are clearly unable to explain the measured electron densities and ion speeds, since an unreasonably high electric field would be required to reduce the electron density to the measured range at the *Rosetta* location. Hence, we follow the *r*⁻¹ profiles to derive the surface electric field strength range that would explain the measured electron density for each time period. We find electric field strengths for July and November of around 1.1–2.7 and 0.2–0.6 mV m⁻¹, respectively. Stronger electric fields (>3 mV m⁻¹) lead to enhanced ion transport, reducing the plasma density below the observed range.

We now examine the predicted NH₄⁺ density for the same input conditions and range of electric field measurements as in the previous section. Fig. 10 b shows a decreasing trend of NH₄⁺ density with increased electric field strength, similar to the total plasma density. As discussed in Section 3, NH₄⁺ is more sensitive to the change in the electric field strength, and therefore varies over more orders of magnitude than the total ion density. For the total plasma densities measured in the two time periods, the model predicts that this would correspond to 20–40 cm⁻³ of NH₄⁺ on 2015 November 20–21 and 50–110 NH₄⁺ for 2015 July 29–30.

While we know NH₄⁺ was consistently detected in the diamagnetic cavity during the two case study periods while the mass spectrometer was in ion mode (see Figs 8 and 7), the number density of NH₄⁺ cannot be inferred from the ROSINA/DFMS spectra with identified NH₄⁺ peaks. It is likely that even the lowest modelled densities in Fig. 10 b would be detectable by the instrument in the high-resolution mode (see Lewis et al. 2023, Appendix A), and NH₄⁺ detection with DFMS is also dependent on the field of view and energy acceptance window of the instrument. Therefore, it is not possible to use the detection of NH₄⁺ to constrain the electric field in this way.

Fig. 10c shows the ion bulk velocity for the same model runs as Figs 10 a and 10b. For the lowest electric fields (≈1 × 10⁻⁵ Vm⁻¹), the ions are not significantly accelerated above the neutral speed. When the electric field is high (>1 × 10⁻² V m⁻¹), the ions are further accelerated, up to 10 km s⁻¹. Such high speeds have been derived from RPC-ICA measurements (Bergman et al. 2021a), but the present ionospheric modelling suggests that such strong acceleration would lead to electron densities around 6 times lower than what have
<figure><img src="image_7.png" /><figcaption>**Figure 6.** [a] Total ion density for different electric field strengths with $E \propto 1/r$, with (solid line) and without (dashed line) momentum transfer. [b] Ion bulk velocity as a fraction of the neutral velocity for the case in [a]. [c] Same as [a] but for $E \propto 1/r^2$. [d] Ion bulk velocities for the case in [c]. The outgassing used is $Q = 10^{28} \text{ s}^{-1}$, and the neutral composition is 1 per cent NH$_3$ and 99 per cent H$_2$O.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Time series of key data inputs into the model for 2015 November 20–21. (a) Local outgassing from ROSINA/COPS with spikes from spacecraft thruster firings removed and spacecraft-comet distance. (b) NH$_3$ mixing ratio from ROSINA/DFMS neutral mode (black crosses) and integrated counts from scans of NH$_4^+$ by ROSINA/DFMS (method described by Lewis et al. 2023). (c) Total electron density from the RPC-MIP/LAP combined data set or RPC-MIP only (red, Henri et al. 2017) where was unavailable. In all panels, shaded regions show where Rosetta was in the diamagnetic cavity, according to RPC-MAG data (Goetz et al. 2016b).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** [a] Total ion density for different electric field strengths with $E \propto 1/r$, with (solid line) and without (dashed line) momentum transfer. [b] Ion bulk velocity as a fraction of the neutral velocity for the case in [a]. [c] Same as [a] but for $E \propto 1/r^2$. [d] Ion bulk velocities for the case in [c]. The outgassing used is $Q = 10^{28} \text{ s}^{-1}$, and the neutral composition is 1 per cent NH$_3$ and 99 per cent H$_2$O.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Time series of key data inputs into the model for 2015 November 20–21. (a) Local outgassing from ROSINA/COPS with spikes from spacecraft thruster firings removed and spacecraft-comet distance. (b) NH$_3$ mixing ratio from ROSINA/DFMS neutral mode (black crosses) and integrated counts from scans of NH$_4^+$ by ROSINA/DFMS (method described by Lewis et al. 2023). (c) Total electron density from the RPC-MIP/LAP combined data set or RPC-MIP only (red, Henri et al. 2017) where was unavailable. In all panels, shaded regions show where Rosetta was in the diamagnetic cavity, according to RPC-MAG data (Goetz et al. 2016b).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Same as Fig. 7, but for 2015 July 29–30.</figcaption></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Histograms showing the spread of total electron density measurements within the diamagnetic cavity for the two time periods considered. The shaded region shows the data one standard deviation either side of the mean.</figcaption></figure>

<table><caption>**Table 1.** Inputs in the model as used to represent the two time periods discussed in Section 4. Local outgassing rate ($Q$), cometocentric distance $r$ of Rosetta, and NH$_3$ mixing ratio are estimated from the data in Figs 7 and 8. The neutral speed $u_n$ is taken from Biver et al. (2019), and the electron density range is as shown in Fig. 9.</caption><thead><tr><th>Input parameter</th><th>2015 July 29–30</th><th>2015 November 20–21</th></tr></thead><tbody><tr><td>Heliocentric dist. [au]</td><td>1.254</td><td>1.715</td></tr><tr><td><i>Q</i> [s<sup>-1</sup>]</td><td>2.3 × 10<sup>28</sup></td><td>7 × 10<sup>27</sup></td></tr><tr><td><i>r</i>[km]</td><td>180</td><td>135</td></tr><tr><td><i>u</i><sub>n</sub> [km s<sup>-1</sup>]</td><td>0.90</td><td>0.75</td></tr><tr><td>NH<sub>3</sub> %</td><td>0.3</td><td>0.2</td></tr><tr><td><i>n</i><sub>e</sub> [cm<sup>-3</sup>]</td><td>955.1–1517</td><td>544.4–791.8</td></tr></tbody></table>

been measured by RPC-MIP and RPC-LAP. Instead, the observed densities are best explained by ion bulk velocities of 1.2–1.7 km s⁻¹ on 2015 November 20–21 and 2–3 km s⁻¹ for July 29–30. The finding of higher ion velocity for 2015 July compared to 2015 November is as may be expected. Rosetta was further away from the comet during July, allowing the cometary ions more time to accelerate.

# 5 DISCUSSION

## 5.1 Ionospheric composition

One novelty of our approach is the inclusion of NH$_4^+$, the dominant ion species in the inner part of the ionosphere near perihelion. This does not have a significant impact on the total ion density, since the dissociative recombination rate coefficient is similar for all the species, but it does reduce the density of H$_3$O$^+$ since loss through proton transfer to NH$_3$ happens very readily in a collisional coma at high outgassing rates. The NH$_4^+$ density is much more sensitive to the enhanced ion speed than the water ion species. However, it is not possible to directly compare modelled ion densities with counts from ROSINA/DFMS ion mode scans, limiting the possibility for using NH$_4^+$ detections with ion mass spectrometer observations within the diamagnetic cavity to constrain the ambipolar electric field.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Same as Fig. 7, but for 2015 July 29–30.</figcaption></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Histograms showing the spread of total electron density measurements within the diamagnetic cavity for the two time periods considered. The shaded region shows the data one standard deviation either side of the mean.</figcaption></figure>

<table><caption>**Table 1.** Inputs in the model as used to represent the two time periods discussed in Section 4. Local outgassing rate ($Q$), cometocentric distance $r$ of Rosetta, and NH$_3$ mixing ratio are estimated from the data in Figs 7 and 8. The neutral speed $u_n$ is taken from Biver et al. (2019), and the electron density range is as shown in Fig. 9.</caption><thead><tr><th>Input parameter</th><th>2015 July 29–30</th><th>2015 November 20–21</th></tr></thead><tbody><tr><td>Heliocentric dist. [au]</td><td>1.254</td><td>1.715</td></tr><tr><td><i>Q</i> [s<sup>-1</sup>]</td><td>2.3 × 10<sup>28</sup></td><td>7 × 10<sup>27</sup></td></tr><tr><td><i>r</i>[km]</td><td>180</td><td>135</td></tr><tr><td><i>u</i><sub>n</sub> [km s<sup>-1</sup>]</td><td>0.90</td><td>0.75</td></tr><tr><td>NH<sub>3</sub></td><td>%</td><td>0.2</td></tr><tr><td><i>n</i><sub>e</sub> [cm<sup>-3</sup>]</td><td>955.1–1517</td><td>544.4–791.8</td></tr></tbody></table>

been measured by RPC-MIP and RPC-LAP. Instead, the observed densities are best explained by ion bulk velocities of 1.2–1.7 km s⁻¹ on 2015 November 20–21 and 2–3 km s⁻¹ for July 29–30. The finding of higher ion velocity for 2015 July compared to 2015 November is as may be expected. Rosetta was further away from the comet during July, allowing the cometary ions more time to accelerate.

# 5 DISCUSSION

## 5.1 Ionospheric composition

One novelty of our approach is the inclusion of NH$_4^+$, the dominant ion species in the inner part of the ionosphere near perihelion. This does not have a significant impact on the total ion density, since the dissociative recombination rate coefficient is similar for all the species, but it does reduce the density of H$_3$O$^+$ since loss through proton transfer to NH$_3$ happens very readily in a collisional coma at high outgassing rates. The NH$_4^+$ density is much more sensitive to the enhanced ion speed than the water ion species. However, it is not possible to directly compare modelled ion densities with counts from ROSINA/DFMS ion mode scans, limiting the possibility for using NH$_4^+$ detections with ion mass spectrometer observations within the diamagnetic cavity to constrain the ambipolar electric field.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** [a] Total plasma density, [b] NH₄⁺ density, and [c] bulk ion velocity, for model runs for input conditions representing 2015 July 29–30 (red lines) and 2015 November 20–21 (black lines), as a function of electric field surface strength. Field profiles ∝ 1/r (solid lines) and ∝ 1/r² (dashed lines) are shown. Shaded red and grey regions show how we derive the NH₄⁺ density, ion bulk velocity, and electric field strength (when E ∝ 1/r) that corresponds to the measured electron density from the RPC instruments over the two periods considered (see Table 1).</figcaption></figure>

The density of H₂O⁺ is resistant to changes in the electric field, since it is close to photochemical equilibrium (Galand et al. 2016; Heritier et al. 2017b) and therefore unaffected by the changing transport term in the continuity equation. As a result, it is also unaffected by the inclusion of momentum transfer in the model. For the other ion species, the momentum transfer has the impact of interrupting the process of ion acceleration, therefore reducing their transport loss and increasing the density. NH₄⁺ is again more sensitive to momentum transfer than H₃O⁺ in the presence of an electric field.

## 5.2 Ion bulk velocity

We compared the ionospheric simulation with total electron density data from RPC instruments (see Section 4). We find that to explain the measured plasma density in the diamagnetic cavity in 2015 November, our model requires an ambipolar field of around 0.5 mV m⁻¹, leading to a bulk ion speed of ~1.2–1.7 km s⁻¹ at the Rosetta location. On 2015 July 29–30, we derive a stronger electric field, up to 1.5 mV m⁻¹, leading to slightly faster bulk speeds of ~2–3 km s⁻¹. The ion speeds we derive are lower than those derived from RPC-ICA, but are broadly consistent with the measured values from RPC-LAP and from the flux conservation method based on observations of the electron density around perihelion (Vigren et al. 2017).

To understand Fig. 10 further, we return to the continuity equation (equation 5), but now considering the total plasma density (i.e. nᵢ(r) = Σⱼ,s nⱼ,s). By making the assumption that dissociative recombination is negligible (justified by its relatively long time-scale, see Fig. 4), then in steady state the equation then reduces to

$$ \frac{1}{r^2} \frac{d}{dr}(n_i r^2 u_i) = P^{\text{ioni}}, \qquad (33) $$

where uᵢ is the bulk ion velocity (equation 16). This leads to the simple relation between the ion density nᵢ and the ion bulk speed uᵢ,

$$ n_i u_i = \frac{\nu^{\text{ioni}} Q}{4\pi u_n r^2} (r - r_c), \qquad (34) $$

showing that at a given cometocentric distance r, the modelled nᵢ ∝ 1/uᵢ. The results of this section then are not dependent on what we assumed for the momentum transfer rate coefficients or on the included ion-neutral chemistry. The most critical assumptions are then regarding the neutral speed uₙ and the ionization frequency νᵢₐₗ.

For the neutral speed, we have taken the values given by the power-law fits given in Biver et al. (2019) (see Table 1), and assumed that this speed is constant down to the surface of the nucleus. In reality, the neutral gas would be better described with an adiabatic expansion model (Heritier et al. 2017b; Huebner & Markiewicz 2000). We can explore the sensitivity of our results to a slower neutral expansion velocity, taking 400 m s⁻¹ as a lower limit – this is the surface speed used in Heritier et al. (2017b). For the same model runs as in Fig. 10, we find this increases the bulk ion speed derived for 2015 July 29–30 to 3–5 km s⁻¹, and for 2015 November 20–21 to 2–3 km s⁻¹.

The ionization frequency (see Fig. 1) is derived from the appropriate TIMED/SEE data set, adjusted to the heliocentric distance of 67P and time-shifted to account for the change in solar phase. The largest uncertainty in this approach is in the solar flux data themselves, which are up to 20 per cent (Woods et al. 2005). Repeating the analysis of Fig. 10 to include maximum and minimum photoionization frequencies (assuming a 20 per cent uncertainty), leads to electric field estimates of 0.1–1 mV m⁻¹ for November 20–21 and 0.7–4 mV m⁻¹ for July 29–30.

We also neglect any attenuation of the solar flux due to absorption from dust grains beyond the cometocentric distance of Rosetta (Johansson et al. 2017). The effect of this attenuation on our model would be to decrease the photoionization frequency, reducing the production of ions, meaning a lower bulk ion speed is required to produce the same plasma density (see equation 34).

## 5.3 Nature of the ambipolar electric field

We find that the electric field following a r⁻¹ dependence is most plausible (compared to r⁻², see Fig. 10). This is, however, likely to not be the case at very large cometocentric distances as would lead to an unbounded potential. Since our simulation is magnetic-field free, it is only valid within the diamagnetic cavity. It is clear that the electric field becomes more complex close to the boundary and outside of this region owing to the presence of the solar wind, and therefore a simple function of cometocentric distance is unlikely.

In addition, we have assumed a constant electron temperature, but in reality it varies with cometocentric distance due to increased electron-neutral collisions in the dense coma close to the surface. A significant population of cold electrons (~0.1 eV) has been observed through much of the escort phase of Rosetta (Eriksson et al. 2017; Henri et al. 2017; Engelhardt et al. 2018; Wattiaux et al. 2020). They dominate over the warm population particularly post-perihelion, and
<figure><img src="image_12.png" /><figcaption>**Figure 10.** [a] Total plasma density, [b] NH₄⁺ density, and [c] bulk ion velocity, for model runs for input conditions representing 2015 July 29–30 (red lines) and 2015 November 20–21 (black lines), as a function of electric field surface strength. Field profiles ∝ 1/r (solid lines) and ∝ 1/r² (dashed lines) are shown. Shaded red and grey regions show how we derive the NH₄⁺ density, ion bulk velocity, and electric field strength (when E ∝ 1/r) that corresponds to the measured electron density from the RPC instruments over the two periods considered (see Table 1).</figcaption></figure>

The density of H₂O⁺ is resistant to changes in the electric field, since it is close to photochemical equilibrium (Galand et al. 2016; Heritier et al. 2017b) and therefore unaffected by the changing transport term in the continuity equation. As a result, it is also unaffected by the inclusion of momentum transfer in the model. For the other ion species, the momentum transfer has the impact of interrupting the process of ion acceleration, therefore reducing their transport loss and increasing the density. NH₄⁺ is again more sensitive to momentum transfer than H₃O⁺ in the presence of an electric field.

## 5.2 Ion bulk velocity

We compared the ionospheric simulation with total electron density data from RPC instruments (see Section 4). We find that to explain the measured plasma density in the diamagnetic cavity in 2015 November, our model requires an ambipolar field of around 0.5 mV m⁻¹, leading to a bulk ion speed of ~1.2–1.7 km s⁻¹ at the Rosetta location. On 2015 July 29–30, we derive a stronger electric field, up to 1.5 mV m⁻¹, leading to slightly faster bulk speeds of ~2–3 km s⁻¹. The ion speeds we derive are lower than those derived from RPC-ICA, but are broadly consistent with the measured values from RPC-LAP and from the flux conservation method based on observations of the electron density around perihelion (Vigren et al. 2017).

To understand Fig. 10 further, we return to the continuity equation (equation 5), but now considering the total plasma density (i.e. nᵢ(r) = Σⱼ,s nⱼ,s). By making the assumption that dissociative recombination is negligible (justified by its relatively long time-scale, see Fig. 4), then in steady state the equation then reduces to

$$ \frac{1}{r^2} \frac{d}{dr}(n_i r^2 u_i) = P^{\text{ioni}}, \qquad (33) $$

where uᵢ is the bulk ion velocity (equation 16). This leads to the simple relation between the ion density nᵢ and the ion bulk speed uᵢ,

$$ n_i u_i = \frac{\nu^{\text{ioni}} Q}{4\pi u_n r^2} (r - r_c), \qquad (34) $$

showing that at a given cometocentric distance r, the modelled nᵢ ∝ 1/uᵢ. The results of this section then are not dependent on what we assumed for the momentum transfer rate coefficients or on the included ion-neutral chemistry. The most critical assumptions are then regarding the neutral speed uₙ and the ionization frequency νᵢₐₗ.

For the neutral speed, we have taken the values given by the power-law fits given in Biver et al. (2019) (see Table 1), and assumed that this speed is constant down to the surface of the nucleus. In reality, the neutral gas would be better described with an adiabatic expansion model (Heritier et al. 2017b; Huebner & Markiewicz 2000). We can explore the sensitivity of our results to a slower neutral expansion velocity, taking 400 m s⁻¹ as a lower limit – this is the surface speed used in Heritier et al. (2017b). For the same model runs as in Fig. 10, we find this increases the bulk ion speed derived for 2015 July 29–30 to 3–5 km s⁻¹, and for 2015 November 20–21 to 2–3 km s⁻¹.

The ionization frequency (see Fig. 1) is derived from the appropriate TIMED/SEE data set, adjusted to the heliocentric distance of 67P and time-shifted to account for the change in solar phase. The largest uncertainty in this approach is in the solar flux data themselves, which are up to 20 per cent (Woods et al. 2005). Repeating the analysis of Fig. 10 to include maximum and minimum photoionization frequencies (assuming a 20 per cent uncertainty), leads to electric field estimates of 0.1–1 mV m⁻¹ for November 20–21 and 0.7–4 mV m⁻¹ for July 29–30.

We also neglect any attenuation of the solar flux due to absorption from dust grains beyond the cometocentric distance of Rosetta (Johansson et al. 2017). The effect of this attenuation on our model would be to decrease the photoionization frequency, reducing the production of ions, meaning a lower bulk ion speed is required to produce the same plasma density (see equation 34).

## 5.3 Nature of the ambipolar electric field

We find that the electric field following a r⁻¹ dependence is most plausible (compared to r⁻², see Fig. 10). This is, however, likely to not be the case at very large cometocentric distances as would lead to an unbounded potential. Since our simulation is magnetic-field free, it is only valid within the diamagnetic cavity. It is clear that the electric field becomes more complex close to the boundary and outside of this region owing to the presence of the solar wind, and therefore a simple function of cometocentric distance is unlikely.

In addition, we have assumed a constant electron temperature, but in reality it varies with cometocentric distance due to increased electron-neutral collisions in the dense coma close to the surface. A significant population of cold electrons (~0.1 eV) has been observed through much of the escort phase of Rosetta (Eriksson et al. 2017; Henri et al. 2017; Engelhardt et al. 2018; Wattiaux et al. 2020). They dominate over the warm population particularly post-perihelion, and
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Ion density profiles [a] and ion bulk velocities [b] for runs of the model where an electric field $E = 10^{-3} r_c/r$ Vm$^{-1}$ is applied, starting at the surface, from $r = 50$ km, and from $r = 100$ km.</figcaption></figure>

were observed to decrease with increasing cometocentric distance (though comet latitude also plays a key role; Gilet et al. 2020). These observations lead to the question of whether it is realistic to assume an ambipolar electric field in $r^{-1}$, or whether it would be better described with a decreasing or even zero electric field strength close to the surface. A difference in energy between the ions and electrons is necessary for an electron pressure gradient, and therefore an ambipolar electric field, to be set up. This possibility was examined in Vigren & Eriksson (2017), who implemented an electric field that abruptly ‘switches on’ at a radius of 100 km, assuming that below this threshold the electrons are too cold to produce such a field. In contrast, we have assumed throughout this work that the ions are accelerated from the nucleus surface.

To test the impact of implementing a ‘cold zone’ with no ion acceleration on our results, Fig. 11 shows ion density and bulk velocity profiles from the present model but with acceleration from 50 km (blue) and from 100 km (red). This treatment is non-physical given that a decrease in the electric field strength towards the nucleus would be gradual and not abrupt, but it is the simplest approximation to make given the electron temperature profile at 67P is not well constrained. It allows the assessment of the sensitivity of the electron density when no electric field is present at the surface. In Fig. 11, we see that the ion bulk velocity and ion density converge to the same profile in all three cases by around $2r_E$, where $r_E$ is the upper limit of the ‘cold zone’ and the start of the ion acceleration (see Fig. 11). This result is not wholly unexpected when we consider equation (12)

in the case where acceleration is only from $r_E$ to a distance $r$:

$$u_i(r) = \sqrt{u_n^2 + \frac{2qE_c r_c}{m} \ln\left(\frac{r}{r_E}\right)} \quad (35)$$

If $E_c r_c$ is held constant (i.e. the field starts with the same strength but from a different initial $r$), then in the limit $r \gg r_E$ there is no $r_E$ dependence. At the cometocentric distances considered in Section 4 (135 and 180 km, see Table 1) a cold electric field-free zone is unlikely to significantly affect the ion bulk velocities we derive, unless it extends close enough to the spacecraft location. The electron exobase has been shown to be close to or even beyond the diamagnetic cavity boundary (Henri et al. 2017), but to explain the accelerated plasma speed observed beyond the neutral speed within this region, an ambipolar electric field must exist below this limit. Hence, although the cold electron population dominates within the cavity, the warm electron presence must be enough to set up a weak ambipolar field for at least some of the coma below the spacecraft location. Both populations are indeed always observed with RPC-MIP (Wattiaux et al. 2020).

For a higher outgassing comet such as 1P, the electrons are much colder compared with 67P and an ambipolar electric field is less likely to result. Electron temperatures in the unmagnetized region at Halley during the Giotto fly-by were calculated by solving coupled continuity, momentum, and energy equations by Korosmezey et al. (1987), finding values around 100 K (similar to the 0.01 eV cold population at Rosetta). In contrast, the total electron temperature (weighted mean of the warm and cold populations) in the diamagnetic cavity at 67P was on the order of $10^4$ K [~1 eV], according to spectra from RPC-MIP (Wattiaux et al. 2020).

In the lower outgassing regime, such as 67P during the Rosetta escort at larger heliocentric distances, the coma is less dense and the cold electron population less significant than near perihelion (Gilet et al. 2020). The ambipolar electric field has been shown to play a role in the plasma environment at these times, both through the presence of suprathermal electrons leading to UV aurora (Galand et al. 2020; Stephenson et al. 2021) and through trapping of electrons allowing them to be cooled despite the thin coma (Stephenson et al. 2022; Stephenson 2024). However, the total electron density measured before 2015 January or from around 2016 February onwards is well explained even with the assumption that the ions travel at the neutral velocity (Galand et al. 2016; Heritier et al. 2018; Vigren et al. 2019), suggesting no acceleration of cometary ions by the ambipolar field. This could be explained by the low cometocentric distances explored by Rosetta during these times (≤60 km) not allowing for sufficient travel times for measurable acceleration, and by the flattening of the electric potential well as shown by test particle simulations (Stephenson et al. 2023).

# 6 CONCLUSION

We have assessed the effect of ion acceleration on ion densities in the presence of an electric field and momentum transfer, for three key ion species: $H_2O^+$, $H_3O^+$, and $NH_4^+$. Using the 1D ionospheric model we developed, we were able to assess the sensitivity of each species to various electric field profiles and strengths. Expanding from the two water ion case (Vigren & Eriksson 2017), we found that ion acceleration by the ambipolar electric field increases the transport loss and therefore reduces the total ion density. We have also shown that each ion species is influenced differently by the electric field and momentum transfer. High PA ions, such as $NH_4^+$, are particularly sensitive to both processes, and where such ions
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Ion density profiles [a] and ion bulk velocities [b] for runs of the model where an electric field $E = 10^{-3} r_c/r$ Vm$^{-1}$ is applied, starting at the surface, from $r = 50$ km, and from $r = 100$ km.</figcaption></figure>

were observed to decrease with increasing cometocentric distance (though comet latitude also plays a key role; Gilet et al. 2020). These observations lead to the question of whether it is realistic to assume an ambipolar electric field in $r^{-1}$, or whether it would be better described with a decreasing or even zero electric field strength close to the surface. A difference in energy between the ions and electrons is necessary for an electron pressure gradient, and therefore an ambipolar electric field, to be set up. This possibility was examined in Vigren & Eriksson (2017), who implemented an electric field that abruptly ‘switches on’ at a radius of 100 km, assuming that below this threshold the electrons are too cold to produce such a field. In contrast, we have assumed throughout this work that the ions are accelerated from the nucleus surface.

To test the impact of implementing a ‘cold zone’ with no ion acceleration on our results, Fig. 11 shows ion density and bulk velocity profiles from the present model but with acceleration from 50 km (blue) and from 100 km (red). This treatment is non-physical given that a decrease in the electric field strength towards the nucleus would be gradual and not abrupt, but it is the simplest approximation to make given the electron temperature profile at 67P is not well constrained. It allows the assessment of the sensitivity of the electron density when no electric field is present at the surface. In Fig. 11, we see that the ion bulk velocity and ion density converge to the same profile in all three cases by around $2r_E$, where $r_E$ is the upper limit of the ‘cold zone’ and the start of the ion acceleration (see Fig. 11). This result is not wholly unexpected when we consider equation (12)

in the case where acceleration is only from $r_E$ to a distance $r$:

$$u_i(r) = \sqrt{u_n^2 + \frac{2qE_c r_c}{m} \ln\left(\frac{r}{r_E}\right)} \quad (35)$$

If $E_c r_c$ is held constant (i.e. the field starts with the same strength but from a different initial $r$), then in the limit $r \gg r_E$ there is no $r_E$ dependence. At the cometocentric distances considered in Section 4 (135 and 180 km, see Table 1) a cold electric field-free zone is unlikely to significantly affect the ion bulk velocities we derive, unless it extends close enough to the spacecraft location. The electron exobase has been shown to be close to or even beyond the diamagnetic cavity boundary (Henri et al. 2017), but to explain the accelerated plasma speed observed beyond the neutral speed within this region, an ambipolar electric field must exist below this limit. Hence, although the cold electron population dominates within the cavity, the warm electron presence must be enough to set up a weak ambipolar field for at least some of the coma below the spacecraft location. Both populations are indeed always observed with RPC-MIP (Wattiaux et al. 2020).

For a higher outgassing comet such as 1P, the electrons are much colder compared with 67P and an ambipolar electric field is less likely to result. Electron temperatures in the unmagnetized region at Halley during the Giotto fly-by were calculated by solving coupled continuity, momentum, and energy equations by Korosmezey et al. (1987), finding values around 100 K (similar to the 0.01 eV cold population at Rosetta). In contrast, the total electron temperature (weighted mean of the warm and cold populations) in the diamagnetic cavity at 67P was on the order of $10^4$ K [~1 eV], according to spectra from RPC-MIP (Wattiaux et al. 2020).

In the lower outgassing regime, such as 67P during the Rosetta escort at larger heliocentric distances, the coma is less dense and the cold electron population less significant than near perihelion (Gilet et al. 2020). The ambipolar electric field has been shown to play a role in the plasma environment at these times, both through the presence of suprathermal electrons leading to UV aurora (Galand et al. 2020; Stephenson et al. 2021) and through trapping of electrons allowing them to be cooled despite the thin coma (Stephenson et al. 2022; Stephenson 2024). However, the total electron density measured before 2015 January or from around 2016 February onwards is well explained even with the assumption that the ions travel at the neutral velocity (Galand et al. 2016; Heritier et al. 2018; Vigren et al. 2019), suggesting no acceleration of cometary ions by the ambipolar field. This could be explained by the low cometocentric distances explored by Rosetta during these times ($\le 60$ km) not allowing for sufficient travel times for measurable acceleration, and by the flattening of the electric potential well as shown by test particle simulations (Stephenson et al. 2023).

# 6 CONCLUSION

We have assessed the effect of ion acceleration on ion densities in the presence of an electric field and momentum transfer, for three key ion species: $H_2O^+$, $H_3O^+$, and $NH_4^+$. Using the 1D ionospheric model we developed, we were able to assess the sensitivity of each species to various electric field profiles and strengths. Expanding from the two water ion case (Vigren & Eriksson 2017), we found that ion acceleration by the ambipolar electric field increases the transport loss and therefore reduces the total ion density. We have also shown that each ion species is influenced differently by the electric field and momentum transfer. High PA ions, such as $NH_4^+$, are particularly sensitive to both processes, and where such ions
are strongly present may be an indicator of limited ion transport (and vice versa). Momentum transfer reduces the effect of the electric field, which affects both ion composition and electron density.

We find that to reproduce the plasma densities measured by RPC inside the diamagnetic cavity of comet 67P near perihelion, our 1D ionospheric model requires an electric field in $r^{-1}$ of around 2 mV m$^{-1}$ at the surface, leading to bulk ion speeds at Rosetta of 1.4–3.0 km s$^{-1}$. Although likely weakened by efficient electron cooling, this electric field is strong enough to accelerate cometary ions above the neutral speed, as observed by RPC. This may have implications for the nature of the diamagnetic cavity boundary. It also contrasts 67P with the higher outgassing case of comet 1P during the Giotto fly-by, during which the overall electron population was likely to be much colder preventing an ambipolar field from forming (Gan & Cravens 1990). The Giotto and Rosetta missions provide a context for the interpretation of future plasma observations in a diamagnetic cavity, as planned for Comet Interceptor (Snodgrass & Jones 2019, Jones et al. 2024). As the target comet for this mission has not yet been identified, modelling will continue to play an important role in exploring the parameter space to build up a picture of the full plasma environment and its key drivers.

## ACKNOWLEDGEMENTS

We would like to acknowledge the invaluable work of the ROSINA team (PI K. Altwegg), the RPC team (including the PI of RPC-MAG, K.-H. Glassmeier, PI of RPC-LAP, A. Eriksson), the whole ESA Rosetta team, and the ESA Planetary Science Archive team. Rosetta is a European Space Agency (ESA) mission with contributions from its member states and NASA. Work at Imperial College London was supported by the Science and Technology Facilities Council (STFC) of the UK under studentship ST/W507519/1 and grant ST/W001071/1, and by the UK Space Agency (UKSA) under grant ST/X002349/1. Work at the University of Bern was funded by the Canton of Bern and the Swiss National Science Foundation (200020_207312). The analytical models in Appendices B and C were derived by AB.

## DATA AVAILABILITY

The data used in this article are available on the Planetary Science Archive at https://psa.esa.int

## REFERENCES

* Altwegg K. et al., 1993, A&A, 279, 260

* Anicich V., Kim J., Huntress W., 1977, Int. J. Mass Spectrom. Ion Phys., 25, 433

* Balsiger H. et al., 2007, Space Sci. Rev., 128, 745

* Behar E., Nilsson H., Wieser G. S., Nemeth Z., Broiles T. W., Richter I., 2016, Geophys. Res. Lett., 43, 1411

* Bergman S., Stenberg Wieser G., Wieser M., Johansson F. L., Eriksson A., 2020, J. Geophys. Res., 125, e27478

* Bergman S. et al., 2021a, MNRAS, 503, 2733

* Bergman S., Wieser G. S., Wieser M., Nilsson H., Vigren E., Beth A., Masunaga K., Eriksson A., 2021b, MNRAS, 507, 4900

* Beth A. et al., 2016, MNRAS, 462, S562

* Beth A., Galand M., Heritier K. L., 2019, A&A, 630, A47

* Beth A., Galand M., Wedlund C., Eriksson A., 2022, Cometary Ionospheres: An Updated Tutorial

* Biermann L., Brosowski B., Schmidt H. U., 1967, Sol. Phys., 1, 254

* Biver N. et al., 2019, A&A, 630, A19

* Carr C. et al., 2007, Space Sci. Rev., 128, 629

* Cravens T. E., 1987, Adv. Space Res, 7, 147

* Cravens T. E., Cambridge University Press UK 2004, Physics of Solar System Plasmas

* Deca J., Divin A., Henri P., Eriksson A., Markidis S., Olshevsky V., Horányi M., 2017, Phys. Rev. Lett., 118, 205101

* Deca J., Henri P., Divin A., Eriksson A., Galand M., Beth A., Ostaszewski K., Horányi M., 2019, Phys. Rev. Lett., 123, 055101

* Engelhardt I. A., Eriksson A. I., Vigren E., Vallières X., Rubin M., Gilet N., Henri P., 2018, A&A, 616, A51

* Eriksson A. I. et al., 2017, A&A, 605, A15

* Galand M. et al., 2016, MNRAS, 462, S331

* Galand M. et al., 2020, Nat. Astron., 4, 1084

* Gan L., Cravens T. E., 1990, J. Geophys. Res., 95, 6285

* Gilet N. et al., 2020, A&A, 640, A110

* Goetz C. et al., 2016a, MNRAS, 462, S459

* Goetz C. et al., 2016b, A&A, 588, A24

* Gombosi T. I., De Zeeuw D. L., Häberli R. M., Powell K. G., 1996, J. Geophys. Res., 101, 15233

* Hansen K. C. et al., 2016, MNRAS, 462, S491

* Haser L., 1957, Bull. Soc. R. Sci. Liege, 43, 740

* Henri P. et al., 2017, MNRAS, 469, S372

* Heritier K. L. et al., 2017a, MNRAS, 469, S118

* Heritier K. L. et al., 2017b, MNRAS, 469, S427

* Heritier K. L. et al., 2018, A&A, 618, A77

* Huebner W. F., Markiewicz W. J., 2000, Icarus, 148, 594

* Huntress W. T. J., Pinizzotto R. F.J., 2003, J. Chem. Phys., 59, 4742

* Johansson F. L. et al., 2017, MNRAS, 469, S626

* Johansson F. L., Eriksson A. I., Gilet N., Henri P., Wattieaux G., Taylor M. G. G. T., Imhof C., Cipriani F., 2020, A&A, 642, A43

* Jones G., Snodgrass C., Tubiana C., Kueppers M., Kawakita H., Lara L., 2024, Space Sci. Rev., 881, 6

* Korosmezey A., Cravens T. E., Gombosi T. I., Nagy A. F., Mendis D. A., Szegő6 K., 1987, J. Geophys. Res., 92, 7331

* Lewis Z. M. et al., 2023, MNRAS, 523, 6208

* McElroy D., Walsh C., Markwick A. J., Cordiner M. A., Smith K., Millar T. J., 2013, A&A, 550, A36

* Madanian H. et al., 2017, AJ, 153 30

* Mott-Smith H. M., Langmuir I., 1926, Phys. Rev., 28, 727

* Neubauer F. M. et al., 1986, Geophys. Res. Lett, 2, 356

* Nilsson H. et al., 2007, Space Sci. Rev., 128, 671

* Nilsson H. et al., 2017, MNRAS, 469, S252

* Novotný O. et al., 2010, J. Phys. Chem. A, 114, 4870

* Odelstad E., Stenberg-Wieser G., Wieser M., Eriksson A. I., Nilsson H., Johansson F. L., 2017, MNRAS, 469, S568

* Odelstad E. et al., 2018, J. Geophys. Res.: Space Phys., 123, 5870

* Öjekull J. et al., 2004, J. Chem. Phys., 120, 7391

* Rosén S. et al., 2000, Faraday Discuss., 115, 295

* Rubin M., Hansen K. C., Combi M. R., Daldorff L. K. S., Gombosi T. I., Tenishev V. M., 2012, J. Geophys. Res., 117, A06227

* Schunk R., Nagy A., 2009, Ionospheres: Physics, Plasma Physics, and Chemistry Cambridge University Press UK

* Smith D., Adams N. G., Henchman M. J., 2008, J. Chem. Phys., 72, 4951

* Snodgrass C., Jones G. H., 2019, Nat. Commun., 10, 5418

* Stephenson P., et al. 2024, MNRAS, 529 3 2854–2865

* Stephenson P. et al., 2021, A&A, 647, A119

* Stephenson P., Galand M., Deca J., Henri P., Carnielli G., 2022, MNRAS, 511, 4090

* Stephenson P. et al., 2023, 525 4 5041–5065 MNRAS

* Szegő K. et al., 2000, Space Sci. Rev., 94, 429

* Trotignon J. G. et al., 2007, Space Sci. Rev., 128, 713

* Vigren E., Eriksson A. I., 2017, AJ, 153, 150

* Vigren E., Galand M., 2013, ApJ, 772, 33

* Vigren E., Galand M., Eriksson A. I., Edberg N. J. T., Odelstad E., Schwartz S. J., 2015, ApJ, 812, 54

* Vigren E. et al., 2017, MNRAS, 469, S142

* Vigren E. et al., 2019, ApJ, 881, 6

* Wattieaux G., Henri P., Gilet N., Vallières X., Deca J., 2020, A&A, 638, A124

* Woods T. N. et al., 2005, J. Geophys. Res., 110, A01312
are strongly present may be an indicator of limited ion transport (and vice versa). Momentum transfer reduces the effect of the electric field, which affects both ion composition and electron density.

We find that to reproduce the plasma densities measured by RPC inside the diamagnetic cavity of comet 67P near perihelion, our 1D ionospheric model requires an electric field in $r^{-1}$ of around 2 mV m$^{-1}$ at the surface, leading to bulk ion speeds at Rosetta of 1.4–3.0 km s$^{-1}$. Although likely weakened by efficient electron cooling, this electric field is strong enough to accelerate cometary ions above the neutral speed, as observed by RPC. This may have implications for the nature of the diamagnetic cavity boundary. It also contrasts 67P with the higher outgassing case of comet 1P during the Giotto fly-by, during which the overall electron population was likely to be much colder preventing an ambipolar field from forming (Gan & Cravens 1990). The Giotto and Rosetta missions provide a context for the interpretation of future plasma observations in a diamagnetic cavity, as planned for Comet Interceptor (Snodgrass & Jones 2019, Jones et al. 2024). As the target comet for this mission has not yet been identified, modelling will continue to play an important role in exploring the parameter space to build up a picture of the full plasma environment and its key drivers.

## ACKNOWLEDGEMENTS

We would like to acknowledge the invaluable work of the ROSINA team (PI K. Altwegg), the RPC team (including the PI of RPC-MAG, K.-H. Glassmeier, PI of RPC-LAP, A. Eriksson), the whole ESA Rosetta team, and the ESA Planetary Science Archive team. Rosetta is a European Space Agency (ESA) mission with contributions from its member states and NASA. Work at Imperial College London was supported by the Science and Technology Facilities Council (STFC) of the UK under studentship ST/W507519/1 and grant ST/W001071/1, and by the UK Space Agency (UKSA) under grant ST/X002349/1. Work at the University of Bern was funded by the Canton of Bern and the Swiss National Science Foundation (200020_207312). The analytical models in Appendices B and C were derived by AB.

## DATA AVAILABILITY

The data used in this article are available on the Planetary Science Archive at https://psa.esa.int

## REFERENCES

* Altwegg K. et al., 1993, A&A, 279, 260

* Anicich V., Kim J., Huntress W., 1977, Int. J. Mass Spectrom. Ion Phys., 25, 433

* Balsiger H. et al., 2007, Space Sci. Rev., 128, 745

* Behar E., Nilsson H., Wieser G. S., Nemeth Z., Broiles T. W., Richter I., 2016, Geophys. Res. Lett., 43, 1411

* Bergman S., Stenberg Wieser G., Wieser M., Johansson F. L., Eriksson A., 2020, J. Geophys. Res., 125, e27478

* Bergman S. et al., 2021a, MNRAS, 503, 2733

* Bergman S., Wieser G. S., Wieser M., Nilsson H., Vigren E., Beth A., Masunaga K., Eriksson A., 2021b, MNRAS, 507, 4900

* Beth A. et al., 2016, MNRAS, 462, S562

* Beth A., Galand M., Heritier K. L., 2019, A&A, 630, A47

* Beth A., Galand M., Wedlund C., Eriksson A., 2022, Cometary Ionospheres: An Updated Tutorial

* Biermann L., Brosowski B., Schmidt H. U., 1967, Sol. Phys., 1, 254

* Biver N. et al., 2019, A&A, 630, A19

* Carr C. et al., 2007, Space Sci. Rev., 128, 629

* Cravens T. E., 1987, Adv. Space Res, 7, 147

* Cravens T. E., Cambridge University Press UK 2004, Physics of Solar System Plasmas

* Deca J., Divin A., Henri P., Eriksson A., Markidis S., Olshevsky V., Horányi M., 2017, Phys. Rev. Lett., 118, 205101

* Deca J., Henri P., Divin A., Eriksson A., Galand M., Beth A., Ostaszewski K., Horányi M., 2019, Phys. Rev. Lett., 123, 055101

* Engelhardt I. A., Eriksson A. I., Vigren E., Vallières X., Rubin M., Gilet N., Henri P., 2018, A&A, 616, A51

* Eriksson A. I. et al., 2017, A&A, 605, A15

* Galand M. et al., 2016, MNRAS, 462, S331

* Galand M. et al., 2020, Nat. Astron., 4, 1084

* Gan L., Cravens T. E., 1990, J. Geophys. Res., 95, 6285

* Gilet N. et al., 2020, A&A, 640, A110

* Goetz C. et al., 2016a, MNRAS, 462, S459

* Goetz C. et al., 2016b, A&A, 588, A24

* Gombosi T. I., De Zeeuw D. L., Häberli R. M., Powell K. G., 1996, J. Geophys. Res., 101, 15233

* Hansen K. C. et al., 2016, MNRAS, 462, S491

* Haser L., 1957, Bull. Soc. R. Sci. Liege, 43, 740

* Henri P. et al., 2017, MNRAS, 469, S372

* Heritier K. L. et al., 2017a, MNRAS, 469, S118

* Heritier K. L. et al., 2017b, MNRAS, 469, S427

* Heritier K. L. et al., 2018, A&A, 618, A77

* Huebner W. F., Markiewicz W. J., 2000, Icarus, 148, 594

* Huntress W. T. J., Pinizzotto R. F.J., 2003, J. Chem. Phys., 59, 4742

* Johansson F. L. et al., 2017, MNRAS, 469, S626

* Johansson F. L., Eriksson A. I., Gilet N., Henri P., Wattieaux G., Taylor M. G. G. T., Imhof C., Cipriani F., 2020, A&A, 642, A43

* Jones G., Snodgrass C., Tubiana C., Kueppers M., Kawakita H., Lara L., 2024, Space Sci. Rev., 881, 6

* Korosmezey A., Cravens T. E., Gombosi T. I., Nagy A. F., Mendis D. A., Szegő6 K., 1987, J. Geophys. Res., 92, 7331

* Lewis Z. M. et al., 2023, MNRAS, 523, 6208

* McElroy D., Walsh C., Markwick A. J., Cordiner M. A., Smith K., Millar T. J., 2013, A&A, 550, A36

* Madanian H. et al., 2017, AJ, 153 30

* Mott-Smith H. M., Langmuir I., 1926, Phys. Rev., 28, 727

* Neubauer F. M. et al., 1986, Geophys. Res. Lett, 2, 356

* Nilsson H. et al., 2007, Space Sci. Rev., 128, 671

* Nilsson H. et al., 2017, MNRAS, 469, S252

* Novotný O. et al., 2010, J. Phys. Chem. A, 114, 4870

* Odelstad E., Stenberg-Wieser G., Wieser M., Eriksson A. I., Nilsson H., Johansson F. L., 2017, MNRAS, 469, S568

* Odelstad E. et al., 2018, J. Geophys. Res.: Space Phys., 123, 5870

* Öjekull J. et al., 2004, J. Chem. Phys., 120, 7391

* Rosén S. et al., 2000, Faraday Discuss., 115, 295

* Rubin M., Hansen K. C., Combi M. R., Daldorff L. K. S., Gombosi T. I., Tenishev V. M., 2012, J. Geophys. Res., 117, A06227

* Schunk R., Nagy A., 2009, Ionospheres: Physics, Plasma Physics, and Chemistry Cambridge University Press UK

* Smith D., Adams N. G., Henchman M. J., 2008, J. Chem. Phys., 72, 4951

* Snodgrass C., Jones G. H., 2019, Nat. Commun., 10, 5418

* Stephenson P., et al. 2024, MNRAS, 529 3 2854–2865

* Stephenson P. et al., 2021, A&A, 647, A119

* Stephenson P., Galand M., Deca J., Henri P., Carnielli G., 2022, MNRAS, 511, 4090

* Stephenson P. et al., 2023, 525 4 5041–5065 MNRAS

* Szegő K. et al., 2000, Space Sci. Rev., 94, 429

* Trotignon J. G. et al., 2007, Space Sci. Rev., 128, 713

* Vigren E., Eriksson A. I., 2017, AJ, 153, 150

* Vigren E., Galand M., 2013, ApJ, 772, 33

* Vigren E., Galand M., Eriksson A. I., Edberg N. J. T., Odelstad E., Schwartz S. J., 2015, ApJ, 812, 54

* Vigren E. et al., 2017, MNRAS, 469, S142

* Vigren E. et al., 2019, ApJ, 881, 6

* Wattieaux G., Henri P., Gilet N., Vallières X., Deca J., 2020, A&A, 638, A124

* Woods T. N. et al., 2005, J. Geophys. Res., 110, A01312
## APPENDIX A: REACTIONS INCLUDED IN THE IONOSPHERIC MODEL

The rate coefficients for dissociative recombination, momentum transfer and ion-neutral chemistry are given in Tables A1, A2, and A3. The

<table><caption>**Table A1.** Dissociative recombination rate coefficients used in the model. $T_e$ is the electron temperature in K.</caption><thead><tr><th>Reaction</th><th>kinetic rate coefficient α<sup>DR</sup><sub>j</sub> [cm<sup>3</sup> s<sup>-1</sup>]</th><th>Temp. range [K]</th><th>Reference</th></tr></thead><tbody><tr><td>H<sub>2</sub>O<sup>+</sup> + e<sup>-</sup> → O + H<sub>2</sub></td><td>3.9 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>10–41000</td><td>Rosén et al. 2000</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + e<sup>-</sup> → O + H + H</td><td>3.05 × 10<sup>-7</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>10–1000</td><td>Rosén et al. 2000</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + e<sup>-</sup> → OH + H</td><td>8.6 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>10–1000</td><td>Rosén et al. 2000</td></tr><tr><td>H<sub>3</sub>O<sup>+</sup> + e<sup>-</sup> → H<sub>2</sub>O + H</td><td>7.09 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>10–1000</td><td>Novotný et al. 2010</td></tr><tr><td>H<sub>3</sub>O<sup>+</sup> + e<sup>-</sup> → O + H<sub>2</sub> + H</td><td>5.60 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>10–1000</td><td>Novotný et al. 2010</td></tr><tr><td>H<sub>3</sub>O<sup>+</sup> + e<sup>-</sup> → OH + H<sub>2</sub></td><td>5.37 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>10–1000</td><td>Novotný et al. 2010</td></tr><tr><td>H<sub>3</sub>O<sup>+</sup> + e<sup>-</sup> → OH + H + H</td><td>3.05 × 10<sup>-7</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>10–1000</td><td>Novotný et al. 2010</td></tr><tr><td>NH<sub>4</sub><sup>+</sup> + e<sup>-</sup> → NH<sub>2</sub> + H<sub>2</sub></td><td>4.72 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.6</sup></td><td>10–2000</td><td>Öjekull et al. 2004</td></tr><tr><td>NH<sub>4</sub><sup>+</sup> + e<sup>-</sup> → NH<sub>2</sub> + H + H</td><td>3.77 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.6</sup></td><td>10–2000</td><td>Öjekull et al. 2004</td></tr><tr><td>NH<sub>4</sub><sup>+</sup> + e<sup>-</sup> → NH<sub>3</sub> + H</td><td>8.49 × 10<sup>-7</sup>(300/T<sub>e</sub>)<sup>0.6</sup></td><td>10–2000</td><td>Öjekull et al. 2004</td></tr></tbody></table>

<table><caption>**Table A2.** Ion-neutral collision coefficients for processes included in the model, assumed to be entirely elastic. The star (★) denotes a fast ion or neutral.</caption><thead><tr><th>Reaction</th><th>Collision rate coefficient k<sup>MT</sup> [cm<sup>3</sup>s<sup>-1</sup>]</th><th>Reference</th></tr></thead><tbody><tr><td>H<sub>2</sub>O<sup>+★</sup> + H<sub>2</sub>O → H<sub>2</sub>O<sup>+</sup> + H<sub>2</sub>O<sup>*</sup></td><td>1.7 × 10<sup>-9</sup></td><td>Gombosi et al. 1996</td></tr><tr><td>H<sub>3</sub>O<sup>+★</sup> + H<sub>2</sub>O → H<sub>3</sub>O<sup>+</sup> + H<sub>2</sub>O<sup>*</sup></td><td>5.04 × 10<sup>-10</sup></td><td>Schunk & Nagy 2009</td></tr><tr><td>NH<sub>4</sub><sup>+</sup><sup>★</sup> + H<sub>2</sub>O → NH<sub>4</sub><sup>+</sup> + H<sub>2</sub>O<sup>★</sup></td><td>5.25 × 10<sup>-10</sup></td><td>Schunk & Nagy 2009</td></tr></tbody></table>

<table><caption>**Table A3.** Ion-neutral chemical reaction rates for reactions included in the model.</caption><thead><tr><th>Reaction</th><th>Kinetic rate coefficient k<sup>IN</sup><sub>j,n→j'</sub> [cm<sup>3</sup> s<sup>-1</sup>]</th><th>Temp. range [K]</th><th>Reference</th></tr></thead><tbody><tr><td>H<sub>2</sub>O<sup>+★</sup> + H<sub>2</sub>O → H<sub>3</sub>O<sup>+</sup> + OH<sup>★</sup></td><td>2.10 × 10<sup>-9</sup>(300/T<sub>n</sub>)<sup>0.5</sup></td><td>10–41000</td><td>Huntress & Pinizzotto 2003</td></tr><tr><td>H<sub>3</sub>O<sup>+★</sup> + NH<sub>3</sub> → NH<sub>4</sub><sup>+</sup> + H<sub>2</sub>O<sup>★</sup></td><td>2.20 × 10<sup>-9</sup>(300/T<sub>n</sub>)<sup>0.5</sup></td><td>10–41000</td><td>Smith, Adams & Henchman 2008</td></tr><tr><td>H<sub>2</sub>O<sup>+★</sup> + NH<sub>3</sub> → NH<sub>4</sub><sup>+</sup> + OH<sup>★</sup></td><td>9.45 × 10<sup>-10</sup>(300/T<sub>n</sub>)<sup>0.5</sup></td><td>10–41000</td><td>Anicich, Kim & Huntress 1977</td></tr></tbody></table>

dissociative recombination and ion-neutral chemical reaction rates are taken from the UMIST data base (McElroy et al. 2013) and their accuracy is within 25 per cent.

## APPENDIX B: ANALYTICAL MODEL II

To solve equations (21)–(23), we first consider the case without photoabsorption, i.e. $\tau = 0$, and set $n_j = \frac{\nu Q}{4\pi u_n^2 r_c} N_j$ with $j = 1, 2, 3$ corresponding to H<sub>2</sub>O<sup>+</sup>, H<sub>3</sub>O<sup>+</sup>, and NH<sub>4</sub><sup>+</sup> respectively.

$$ \frac{\nu Q}{4\pi u_n^2 r_c} \frac{1}{r^2} \frac{dN_1 r^2}{dr} = \frac{\nu Q}{4\pi u_n^2 r_c^2} - \frac{\nu Q}{4\pi u_n^2 r_c} (k_1 f_{H_2O} + k_2 f_{NH_3}) \frac{Q}{4\pi u_n r^2} N_1 \quad (B1) $$

$$ \frac{\nu Q}{4\pi u_n^2 r_c} \frac{1}{r^2} \frac{dN_2 r^2}{dr} = \frac{\nu Q}{4\pi u_n^2 r_c} (k_1 f_{H_2O}) \frac{Q}{4\pi u_n r^2} N_1 - \frac{\nu Q}{4\pi u_n^2 r_c} (k_3 f_{NH_3}) \frac{Q}{4\pi u_n r^2} N_2 \quad (B2) $$

$$ \frac{\nu Q}{4\pi u_n^2 r_c} \frac{1}{r^2} \frac{dN_3 r^2}{dr} = \frac{\nu Q}{4\pi u_n^2 r_c} (k_2 f_{NH_3}) \frac{Q}{4\pi u_n^2 r^2} N_1 + \frac{\nu Q}{4\pi u_n^2 r_c} (k_3 f_{NH_3}) \frac{Q}{4\pi u_n r^2} N_2 \quad (B3) $$

We then set $x = r_c/r$ and $y_i = n_i x^2$ to obtain the simplified set of differential equations:

$$ \frac{dy_1}{dx} = 1 - \frac{B y_1}{x^2} \quad (B4) $$
## APPENDIX A: REACTIONS INCLUDED IN THE IONOSPHERIC MODEL

The rate coefficients for dissociative recombination, momentum transfer and ion-neutral chemistry are given in Tables A1, A2, and A3. The

<table><caption>**Table A1.** Dissociative recombination rate coefficients used in the model. $T_e$ is the electron temperature in K.</caption><thead><tr><th>Reaction</th><th>kinetic rate coefficient α<sup>DR</sup><sub>j</sub> [cm<sup>3</sup> s<sup>-1</sup>]</th><th>Temp. range [K]</th><th>Reference</th></tr></thead><tbody><tr><td>H<sub>2</sub>O<sup>+</sup> + e<sup>-</sup> → O + H<sub>2</sub></td><td>3.9 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>10–41000</td><td>Rosén et al. 2000</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + e<sup>-</sup> → O + H + H</td><td>3.05 × 10<sup>-7</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>10–1000</td><td>Rosén et al. 2000</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + e<sup>-</sup> → OH + H</td><td>8.6 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>10–1000</td><td>Rosén et al. 2000</td></tr><tr><td>H<sub>3</sub>O<sup>+</sup> + e<sup>-</sup> → H<sub>2</sub>O + H</td><td>7.09 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>10–1000</td><td>Novotný et al. 2010</td></tr><tr><td>H<sub>3</sub>O<sup>+</sup> + e<sup>-</sup> → O + H<sub>2</sub> + H</td><td>5.60 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>10–1000</td><td>Novotný et al. 2010</td></tr><tr><td>H<sub>3</sub>O<sup>+</sup> + e<sup>-</sup> → OH + H<sub>2</sub></td><td>5.37 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>10–1000</td><td>Novotný et al. 2010</td></tr><tr><td>H<sub>3</sub>O<sup>+</sup> + e<sup>-</sup> → OH + H + H</td><td>3.05 × 10<sup>-7</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>10–1000</td><td>Novotný et al. 2010</td></tr><tr><td>NH<sub>4</sub><sup>+</sup> + e<sup>-</sup> → NH<sub>2</sub> + H<sub>2</sub></td><td>4.72 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.6</sup></td><td>10–2000</td><td>Öjekull et al. 2004</td></tr><tr><td>NH<sub>4</sub><sup>+</sup> + e<sup>-</sup> → NH<sub>2</sub> + H + H</td><td>3.77 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.6</sup></td><td>10–2000</td><td>Öjekull et al. 2004</td></tr><tr><td>NH<sub>4</sub><sup>+</sup> + e<sup>-</sup> → NH<sub>3</sub> + H</td><td>8.49 × 10<sup>-7</sup>(300/T<sub>e</sub>)<sup>0.6</sup></td><td>10–2000</td><td>Öjekull et al. 2004</td></tr></tbody></table>

<table><caption>**Table A2.** Ion-neutral collision coefficients for processes included in the model, assumed to be entirely elastic. The star (★) denotes a fast ion or neutral.</caption><thead><tr><th>Reaction</th><th>Collision rate coefficient k<sup>MT</sup> [cm<sup>3</sup>s<sup>-1</sup>]</th><th>Reference</th></tr></thead><tbody><tr><td>H<sub>2</sub>O<sup>+★</sup> + H<sub>2</sub>O → H<sub>2</sub>O<sup>+</sup> + H<sub>2</sub>O<sup>*</sup></td><td>1.7 × 10<sup>-9</sup></td><td>Gombosi et al. 1996</td></tr><tr><td>H<sub>3</sub>O<sup>+★</sup> + H<sub>2</sub>O → H<sub>3</sub>O<sup>+</sup> + H<sub>2</sub>O<sup>*</sup></td><td>5.04 × 10<sup>-10</sup></td><td>Schunk & Nagy 2009</td></tr><tr><td>NH<sub>4</sub><sup>+</sup><sup>★</sup> + H<sub>2</sub>O → NH<sub>4</sub><sup>+</sup> + H<sub>2</sub>O<sup>★</sup></td><td>5.25 × 10<sup>-10</sup></td><td>Schunk & Nagy 2009</td></tr></tbody></table>

<table><caption>**Table A3.** Ion-neutral chemical reaction rates for reactions included in the model.</caption><thead><tr><th>Reaction</th><th>Kinetic rate coefficient k<sup>IN</sup><sub>j,n→j'</sub> [cm<sup>3</sup> s<sup>-1</sup>]</th><th>Temp. range [K]</th><th>Reference</th></tr></thead><tbody><tr><td>H<sub>2</sub>O<sup>+★</sup> + H<sub>2</sub>O → H<sub>3</sub>O<sup>+</sup> + OH<sup>★</sup></td><td>2.10 × 10<sup>-9</sup>(300/T<sub>n</sub>)<sup>0.5</sup></td><td>10–41000</td><td>Huntress & Pinizzotto 2003</td></tr><tr><td>H<sub>3</sub>O<sup>+★</sup> + NH<sub>3</sub> → NH<sub>4</sub><sup>+</sup> + H<sub>2</sub>O<sup>★</sup></td><td>2.20 × 10<sup>-9</sup>(300/T<sub>n</sub>)<sup>0.5</sup></td><td>10–41000</td><td>Smith, Adams & Henchman 2008</td></tr><tr><td>H<sub>2</sub>O<sup>+★</sup> + NH<sub>3</sub> → NH<sub>4</sub><sup>+</sup> + OH<sup>★</sup></td><td>9.45 × 10<sup>-10</sup>(300/T<sub>n</sub>)<sup>0.5</sup></td><td>10–41000</td><td>Anicich, Kim & Huntress 1977</td></tr></tbody></table>

dissociative recombination and ion-neutral chemical reaction rates are taken from the UMIST data base (McElroy et al. 2013) and their accuracy is within 25 per cent.

## APPENDIX B: ANALYTICAL MODEL II

To solve equations (21)–(23), we first consider the case without photoabsorption, i.e. $\tau = 0$, and set $n_j = \frac{\nu Q}{4\pi u_n^2 r_c} N_j$ with $j = 1, 2, 3$ corresponding to $H_2O^+$, $H_3O^+$, and $NH_4^+$ respectively.

$$ \frac{\nu Q}{4\pi u_n^2 r_c} \frac{1}{r^2} \frac{dN_1 r^2}{dr} = \frac{\nu Q}{4\pi u_n^2 r_c} \frac{1}{r^2} \frac{dN_2 r^2}{dr} - \frac{\nu Q}{4\pi u_n^2 r_c} (k_1 f_{H_2O} + k_2 f_{NH_3}) \frac{Q}{4\pi u_n r^2} N_1 \quad (B1) $$

$$ \frac{\nu Q}{4\pi u_n^2 r_c} \frac{1}{r^2} \frac{dN_2 r^2}{dr} = \frac{\nu Q}{4\pi u_n^2 r_c} (k_1 f_{H_2O}) \frac{Q}{4\pi u_n r^2} N_1 - \frac{\nu Q}{4\pi u_n^2 r_c} (k_3 f_{NH_3}) \frac{Q}{4\pi u_n r^2} N_2 \quad (B2) $$

$$ \frac{\nu Q}{4\pi u_n^2 r_c} \frac{1}{r^2} \frac{dN_3 r^2}{dr} = \frac{\nu Q}{4\pi u_n^2 r_c} (k_2 f_{NH_3}) \frac{Q}{4\pi u_n^2 r^2} N_1 + \frac{\nu Q}{4\pi u_n^2 r_c} (k_3 f_{NH_3}) \frac{Q}{4\pi u_n r^2} N_2 \quad (B3) $$

We then set $x = r_c/r$ and $y_i = n_i x^2$ to obtain the simplified set of differential equations:

$$ \frac{dy_1}{dx} = 1 - \frac{B y_1}{x^2} \quad (B4) $$
$$ \frac{dy_2}{dx} = \frac{Cy_1}{x^2} - \frac{Dy_2}{x^2} \qquad (B5) $$

$$ \frac{dy_3}{dx} = \frac{Ey_1}{x^2} + \frac{Fy_2}{x^2} \qquad (B6) $$

which can be rewritten as

$$ \frac{dy_1 \exp(-B/x)}{dx} = \exp(-B/x), \qquad (B7) $$

$$ \frac{dy_2 \exp(-D/x)}{dx} = \frac{Cy_1}{x^2} \exp(-D/x), \qquad (B8) $$

$$ \frac{dy_3}{dx} = \frac{Ey_1}{x^2} + \frac{Fy_2}{x^2}. \qquad (B9) $$

By setting $Y_1 = y_1 \exp(-B/x)$ and $Y_2 = y_2 \exp(-D/x)$, and using $Ei(-x) = -E_1(x)$ and $E_2(x) = \exp(-x) - E_1(x)$, these can now be solved for $Y_1$, $Y_2$, and $Y_3$:

$$ Y_1 = xE_2(B/x) - E_2(B) \qquad (B10) $$

$$ Y_2 \exp(D/x) = - \frac{C}{B-D} [xE_2(B/x) - E_2(B)] \exp(B/x) + \frac{C}{B-D} [xE_2(D/x) - E_2(D)] \exp(D/x) \qquad (B11) $$

With photoabsorption (equations 21–23), one has to replace $B$ with $B + \tau$ and $D$ with $D + \tau$ in equations (B10) and (B11). Putting it all together, we obtain:

$$ n_i = \frac{\nu Q}{4\pi u_n^2 r_c} \frac{r_c^2}{r^2} \left[ \frac{r}{r_c} E_2 \left( \tau \frac{r_c}{r} \right) - E_2(\tau) \right] \qquad (B12) $$

$$ n_{H_2O^+} = \frac{\nu Q}{4\pi u_n^2 r_c} \frac{r_c^2}{r^2} \left[ \frac{r}{r_c} E_2 \left( (B + \tau) \frac{r_c}{r} \right) - E_2(B + \tau) \right] \exp \left( B \frac{r_c}{r} \right) \qquad (B13) $$

$$ n_{H_3O^+} = \frac{\nu Q}{4\pi u_n^2 r_c} \frac{r_c^2}{r^2} \frac{k_1 f_{H_2O}}{k_1 f_{H_2O} + (k_2 - k_3) f_{NH_3}} \left( \left[ \frac{r}{r_c} E_2 \left( (B + \tau) \frac{r_c}{r} \right) - E_2(B + \tau) \right] \exp \left( B \frac{r_c}{r} \right) - \left[ \frac{r}{r_c} E_2 \left( (D + \tau) \frac{r_c}{r} \right) - E_2(D + \tau) \right] \exp \left( D \frac{r_c}{r} \right) \right) \qquad (B14) $$

$$ n_{NH_4^+} = n_i - n_{H_2O^+} - n_{H_3O^+} \qquad (B15) $$

where the constants are given by

$$ B = \frac{(k_1 f_{H_2O} + k_2 f_{NH_3}) Q}{4\pi u_n^2 r_c} \qquad (B16) $$

$$ C = \frac{k_1 f_{H_2O} Q}{4\pi u_n^2 r_c} \qquad (B17) $$

$$ D = \frac{k_3 f_{NH_3} Q}{4\pi u_n^2 r_c} \qquad (B18) $$

$$ \tau = \frac{\sigma Q}{4\pi u_n r_c}. \qquad (B19) $$

## APPENDIX C: ANALYTICAL MODEL III

Assuming now that ions are still produced radially at $u_n$, the ions accelerate from $u_n$ at their birthplace $r_s$ to $u_j$ at the final location $r$. By the conservation of mechanical energy, $u_j$ is given by

$$ \frac{1}{2} m u_j^2(r) + q E_{\text{amb}}(r_c) \frac{r_c^2}{r} = \frac{1}{2} m u_n^2 + q E_{\text{amb}}(r_c) \frac{r_c^2}{r_s} $$

that is

$$ u_j(r) = u_n \sqrt{1 + \frac{2qE_{\text{amb}}(r_c)r_c}{mu_n^2} \left( \frac{r_c}{r_s} - \frac{r_c}{r} \right)} \qquad (C1) $$
$$ \frac{dy_2}{dx} = \frac{Cy_1}{x^2} - \frac{Dy_2}{x^2} \qquad (B5) $$

$$ \frac{dy_3}{dx} = \frac{Ey_1}{x^2} + \frac{Fy_2}{x^2} \qquad (B6) $$

which can be rewritten as

$$ \frac{dy_1 \exp(-B/x)}{dx} = \exp(-B/x), \qquad (B7) $$

$$ \frac{dy_2 \exp(-D/x)}{dx} = \frac{Cy_1}{x^2} \exp(-D/x), \qquad (B8) $$

$$ \frac{dy_3}{dx} = \frac{Ey_1}{x^2} + \frac{Fy_2}{x^2}. \qquad (B9) $$

By setting $Y_1 = y_1 \exp(-B/x)$ and $Y_2 = y_2 \exp(-D/x)$, and using $Ei(-x) = -E_1(x)$ and $E_2(x) = \exp(-x) - E_1(x)$, these can now be solved for $Y_1$, $Y_2$, and $Y_3$:

$$ Y_1 = xE_2(B/x) - E_2(B) \qquad (B10) $$

$$ Y_2 \exp(D/x) = - \frac{C}{B-D} [xE_2(B/x) - E_2(B)] \exp(B/x) + \frac{C}{B-D} [xE_2(D/x) - E_2(D)] \exp(D/x) \qquad (B11) $$

With photoabsorption (equations 21–23), one has to replace $B$ with $B + \tau$ and $D$ with $D + \tau$ in equations (B10) and (B11). Putting it all together, we obtain:

$$ n_i = \frac{\nu Q}{4\pi u_n^2 r_c} \frac{r_c^2}{r^2} \left[ \frac{r}{r_c} E_2 \left( \tau \frac{r_c}{r} \right) - E_2(\tau) \right] \qquad (B12) $$

$$ n_{H_2O^+} = \frac{\nu Q}{4\pi u_n^2 r_c} \frac{r_c^2}{r^2} \left[ \frac{r}{r_c} E_2 \left( (B + \tau) \frac{r_c}{r} \right) - E_2(B + \tau) \right] \exp \left( B \frac{r_c}{r} \right) \qquad (B13) $$

$$ n_{H_3O^+} = \frac{\nu Q}{4\pi u_n^2 r_c} \frac{r_c^2}{r^2} \frac{k_1 f_{H_2O}}{k_1 f_{H_2O} + (k_2 - k_3) f_{NH_3}} \left( \left[ \frac{r}{r_c} E_2 \left( (B + \tau) \frac{r_c}{r} \right) - E_2(B + \tau) \right] \exp \left( B \frac{r_c}{r} \right) - \left[ \frac{r}{r_c} E_2 \left( (D + \tau) \frac{r_c}{r} \right) - E_2(D + \tau) \right] \exp \left( D \frac{r_c}{r} \right) \right) \qquad (B14) $$

$$ n_{NH_4^+} = n_i - n_{H_2O^+} - n_{H_3O^+} \qquad (B15) $$

where the constants are given by

$$ B = \frac{(k_1 f_{H_2O} + k_2 f_{NH_3}) Q}{4\pi u_n^2 r_c} \qquad (B16) $$

$$ C = \frac{k_1 f_{H_2O} Q}{4\pi u_n^2 r_c} \qquad (B17) $$

$$ D = \frac{k_3 f_{NH_3} Q}{4\pi u_n^2 r_c} \qquad (B18) $$

$$ \tau = \frac{\sigma Q}{4\pi u_n r_c}. \qquad (B19) $$

## APPENDIX C: ANALYTICAL MODEL III

Assuming now that ions are still produced radially at $u_n$, the ions accelerate from $u_n$ at their birthplace $r_s$ to $u_j$ at the final location $r$. By the conservation of mechanical energy, $u_j$ is given by

$$ \frac{1}{2} m u_j^2(r) + q E_{\text{amb}}(r_c) \frac{r_c^2}{r} = \frac{1}{2} m u_n^2 + q E_{\text{amb}}(r_c) \frac{r_c^2}{r_s} $$

that is

$$ u_j(r) = u_n \sqrt{1 + \frac{2qE_{\text{amb}}(r_c)r_c}{mu_n^2} \left( \frac{r_c}{r_s} - \frac{r_c}{r} \right)} \qquad (C1) $$
The ion number density is then obtained from:

$$n_{\text{ion}}(r) = \frac{v_{\text{ioni}} Q}{4\pi u_n^2 r^2} \int_{r_c}^{r} \frac{1}{\sqrt{1 + \frac{2qE_{\text{amb}}(r_c)r_c}{mu_n^2} \left(\frac{r_c}{r_s} - \frac{r_c}{r}\right)}} dr_s \quad (C2)$$

The integral (C2) can be solved analytically. We define first two constants:

$$\xi = \frac{2qE_{\text{amb}}(r_c)r_c}{mu_n^2} \quad (C3)$$

$$C = \frac{\xi v_{\text{ioni}} Q r_c}{4\pi u_n^2 r^2} \quad (C4)$$

$\xi$ gives the global shape and the $r$ dependence whereas $C$ gives the amplitude. Two cases should be distinguished:

(i) $\xi \frac{r_c}{r} - 1 > 0$

$$A = \sqrt{\frac{2qE_{\text{amb}}(r_c)r_c r_c}{mu_n^2 r} - 1}$$

$$n_{\text{ion}}(r) = \frac{C}{2A^3} \left[ 2 \left( \arctan \left( \frac{\sqrt{\alpha - A^2}}{A} \right) - \arctan \left( \frac{1}{A} \right) \right) + \sin \left( 2 \arctan \left( \frac{\sqrt{\alpha - A^2}}{A} \right) \right) - \sin \left( 2 \arctan \left( \frac{1}{A} \right) \right) \right] \quad (C5)$$

(ii) $\xi \frac{r_c}{r} - 1 < 0$

$$B = \sqrt{1 - \frac{2qE_{\text{amb}}(r_c)r_c r_c}{mu_n^2 r}}$$

$$n_{\text{ion}}(r) = \frac{C}{2B^3} \left[ 2 \left( \arctan\left(\frac{B}{\sqrt{\xi + B^2}}\right) - \arctan\left(B\right) \right) + \sinh\left(2 \arctan\left(B\right)\right) - \sinh\left(2 \arctan\left(\frac{B}{\sqrt{\alpha + B^2}}\right)\right) \right] \quad (C6)$$

For comparison, other formula are available for a different ambipolar electric field profile. For $E_{\text{amb}}(r) = E_{\text{amb}}(r_c)(r_c/r)$, one can get Vigren et al. (2015):

$$n_{\text{ion}}(r) = \frac{v_{\text{ioni}} Q}{4\sqrt{\xi \pi u_n^2 r}} \left[ \erf\left(\sqrt{\frac{1}{\xi} + \log\left(\frac{r}{r_c}\right)}\right) - \erf\left(\sqrt{\frac{1}{\xi}}\right) \right] \exp\left(\frac{1}{\xi}\right) \quad (C7)$$

and for $E_{\text{amb}}(r) = E_{\text{amb}}(r_c)$, one can get

$$n_{\text{ion}}(r) = \frac{v_{\text{ioni}} Q r_c}{2\xi \pi u_n^2 r^2} \left( \sqrt{1 - \xi \left(1 - \frac{r}{r_c}\right)} - 1 \right) \quad (C8)$$

Other generic profiles can be derived with hypergeometric functions for different power laws $E_{\text{amb}} \propto r^{-\gamma}$, beyond the scope of this work. Once the density has been derived, it is straightforward to get the mean ion radial velocity. As the production term in the continuity equation is the same as for the Haser model (Section 2, equation 4) as well as the symmetry, the flux is equal between both models so that the mean velocity is

$$\langle v_{\text{ioni}} \rangle(r) = \frac{n_{\text{ion,H}}(r)}{n_{\text{ion}}(r)} u_n. \quad (C9)$$

This paper has been typeset from a TeX/LATEX file prepared by the author.
The ion number density is then obtained from:

$$n_{\text{ion}}(r) = \frac{v_{\text{ioni}} Q}{4\pi u_n^2 r^2} \int_{r_c}^{r} \frac{1}{\sqrt{1 + \frac{2qE_{\text{amb}}(r_c)r_c}{mu_n^2} \left(\frac{r_c}{r_s} - \frac{r_c}{r}\right)}} dr_i \quad (C2)$$

The integral (C2) can be solved analytically. We define first two constants:

$$\xi = \frac{2qE_{\text{amb}}(r_c)r_c}{mu_n^2} \quad (C3)$$

$$C = \frac{\xi v_{\text{ioni}} Q r_c}{4\pi u_n^2 r^2} \quad (C4)$$

$\xi$ gives the global shape and the $r$ dependence whereas $C$ gives the amplitude. Two cases should be distinguished:

(i) $\xi \frac{r_c}{r} - 1 > 0$

$$A = \sqrt{\frac{2qE_{\text{amb}}(r_c)r_c r_c}{mu_n^2 r} - 1}$$

$$n_{\text{ion}}(r) = \frac{C}{2A^3} \left[ 2 \left( \arctan \left( \frac{\sqrt{\alpha - A^2}}{A} \right) - \arctan \left( \frac{1}{A} \right) \right) + \sin \left( 2 \arctan \left( \frac{\sqrt{\alpha - A^2}}{A} \right) \right) - \sin \left( 2 \arctan \left( \frac{1}{A} \right) \right) \right] \quad (C5)$$

(ii) $\xi \frac{r_c}{r} - 1 < 0$

$$B = \sqrt{1 - \frac{2qE_{\text{amb}}(r_c)r_c r_c}{mu_n^2 r}}$$

$$n_{\text{ion}}(r) = \frac{C}{2B^3} \left[ 2 \left( \arctan h \left( \frac{B}{\sqrt{\xi + B^2}} \right) - \arctan h(B) \right) + \sinh(2 \arctan h(B)) - \sinh \left( 2 \arctan h \left( \frac{B}{\sqrt{\alpha + B^2}} \right) \right) \right] \quad (C6)$$

For comparison, other formula are available for a different ambipolar electric field profile. For $E_{\text{amb}}(r) = E_{\text{amb}}(r_c)(r_c/r)$, one can get Vigren et al. (2015):

$$n_{\text{ion}}(r) = \frac{v_{\text{ioni}} Q}{4\sqrt{\xi\pi u_n^2 r}} \left[ \erf\left(\sqrt{\frac{1}{\xi} + \log\left(\frac{r}{r_c}\right)}\right) - \erf\left(\sqrt{\frac{1}{\xi}}\right) \right] \exp\left(\frac{1}{\xi}\right) \quad (C7)$$

and for $E_{\text{amb}}(r) = E_{\text{amb}}(r_c)$, one can get

$$n_{\text{ion}}(r) = \frac{v_{\text{ioni}} Q r_c}{2\xi\pi u_n^2 r^2} \left( \sqrt{1 - \xi \left( 1 - \frac{r}{r_c} \right)} - 1 \right) \quad (C8)$$

Other generic profiles can be derived with hypergeometric functions for different power laws $E_{\text{amb}} \propto r^{-\gamma}$, beyond the scope of this work. Once the density has been derived, it is straightforward to get the mean ion radial velocity. As the production term in the continuity equation is the same as for the Haser model (Section 2, equation 4) as well as the symmetry, the flux is equal between both models so that the mean velocity is

$$\langle v_{\text{ioni}} \rangle (r) = \frac{n_{\text{ion,H}}(r)}{n_{\text{ion}}(r)} u_n. \quad (C9)$$

This paper has been typeset from a TeX/LATeX file prepared by the author.