# A collisional test-particle model of electrons at a comet

Peter Stephenson, M. Galand, Jan Deca, P. Henri, G. Carnielli

► To cite this version:

Peter Stephenson, M. Galand, Jan Deca, P. Henri, G. Carnielli. A collisional test-particle model of electrons at a comet. Monthly Notices of the Royal Astronomical Society, 2022, 511, pp.4090-4108. 10.1093/mnras/stac055. insu-03611698

HAL Id: insu-03611698

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

Submitted on 17 Mar 2022

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
<figure><img src="image_2.png" /></figure>

# A collisional test-particle model of electrons at a comet

Peter Stephenson<sup>1,1*</sup> M. Galand<sup>1</sup> J. Deca,<sup>2</sup> P. Henri<sup>3,4</sup> and G. Carnielli<sup>1</sup>

Accepted 2021 December 14. Received 2021 December 10; in original form 2021 July 19

## ABSTRACT

We have developed the first 3D collisional model of electrons at a comet, which we use to examine the impact of electron-neutral collisions in the weakly outgassing regime. The test-particle Monte Carlo model uses electric and magnetic fields from a fully kinetic Particle-in-Cell (PiC) model as an input. In our model, electrons originate from the solar wind or from ionization of the neutral coma, either by electron impact or absorption of an extreme ultraviolet photon. All relevant electron-neutral collision processes are included in the model including elastic scattering, excitation, and ionization. Trajectories of electrons are validated against analytically known drifts and the stochastic energy degradation used in the model is compared to the continuous slowing down approximation. Macroscopic properties of the solar wind and cometary electron populations, such as density and temperature, are validated with simple known cases and via comparison with the collisionless PiC model. We demonstrate that electrons are trapped close to the nucleus by the ambipolar electric field, causing an increase in the efficiency of electron-neutral collisions. Even at a low-outgassing rate ($Q = 10^{26} \text{ s}^{-1}$), electron-neutral collisions are shown to cause significant cooling in the coma. The model also provides a multistep numerical framework that is used to assess the influence of the electron-to-ion mass ratio, enabling access to electron dynamics with a physical electron mass.

**Key words:** comets; general – comets; individual: 67P/CG.

## 1 INTRODUCTION

A comet is formed by a solid nucleus, made-up of ices, dust, and rock, surrounded by an extensive envelope of neutral gas and plasma. As the comet approaches the Sun, the outgassing from the icy surface volatiles increases and the neutral coma becomes denser. The plasma within the coma originates either through transport of the solar wind (SW) or from ionization of the cometary molecules. Photoionization by solar extreme-ultraviolet (EUV) photons and electron-impact ionization are the two main sources of cometary electrons at a weakly outgassing comet (Galand et al. 2016; Heritier et al. 2018).

The electrons within the coma are not well described as a single population. They are typically divided into a hot, a warm, and a cold population. During the flyby of 21P/Giaccobini-Zinner by the International Cometary Explorer, distinct electron populations were observed by the electron spectrometer at 4, 10, and 80 eV (Zwickl et al. 1986). Thermal noise spectroscopy during the same flyby revealed that the electron temperature dropped to 1 eV near closest approach (7800 km; Meyer-Vernet et al. 1986) to the highly outgassing comet ($Q = 3 \times 10^{28} \text{ s}^{-1}$; Weaver et al. 1999).

More recently, the Rosetta mission to comet 67P/Churyumov Gerasimenko (hereafter 67P; Glassmeier et al. 2007) has probed the cometary plasma throughout a 2 yr escort phase. The Rosetta Plasma Consortium (RPC; Carr et al. 2007) was a package of five instruments to monitor the local plasma, three of which measured properties of

cometary (and SW) electrons. The Ion and Electron Sensor (Burch et al. 2007) measured the energy distribution of electrons at the spacecraft above energy of at least 4.32 eV (depending on the spacecraft potential). In the coma, this electron population was not thermalized (Clark et al. 2015), but was better described by two kappa distributions (Broiles et al. 2016). Myllys et al. (2019) looked at statistics of the kappa distributions, which fit reasonably well close to perihelion but are less applicable to the larger heliocentric distances of interest here.

The Langmuir Probe (LAP; Eriksson et al. 2007) and Mutual Impedance Probe (MIP; Trotignon et al. 2007) on Rosetta probed the warm population (~5–10 eV) throughout the mission. This was largely made-up of the newly born cometary electrons with a density observed between 10 and $1000 \text{ cm}^{-3}$ (Edberg et al. 2015; Hajra et al. 2020). When comet 67P was close to perihelion, the cometary electrons were mostly produced through photoionization but at larger heliocentric distances electron-impact ionization was at least as important, if not more (Heritier et al. 2017, 2018). There may also be a contribution from core SW electrons (Deca et al. 2017), but the extent of this is not well known. The hot electrons (>40 eV) were much more diffuse ($0.1–1 \text{ cm}^{-3}$) and have been attributed to SW electrons (Broiles et al. 2016; Myllys et al. 2019).

RPC/LAP and MIP observed a cold population (<0.1 eV) of electrons at comet 67P (e.g. Eriksson et al. 2017; Gilet et al. 2017), formed via cooling of the warm population. The cooling was typically caused by collisions with the neutral gas which was at ~0.01 eV (Gulkis et al. 2015; Mandt et al. 2016) and therefore occurred when the electrons were collisionally coupled to the neutral coma. This

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

<footnote n="2"><sup>2</sup>Laboratory for Atmospheric and Space Physics, University of Colorado, Boulder, Colorado, 80303, USA</footnote>

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

<footnote n="4"><sup>4</sup>LPC2E, CNRS, Orleans, 45100, France</footnote>

<footnote>*E-mail: pps18@ic.ac.uk</footnote>
was the case in highly active comets ($Q > 5 \times 10^{27} \text{ s}^{-1}$) such as 21P/G-Z, 1P, and 67P close to perihelion. However, the cold electron population was seen throughout the escort phase even at large heliocentric distances (>3.5 au) and at low-outgassing rates ($Q < 10^{26} \text{ s}^{-1}$) towards the end of mission (Gilet et al. 2020; Wattieaux et al. 2020).

The electron exobase has conventionally been taken to define the boundary below which electron-neutral collisions occur frequently (Mandt et al. 2016) and therefore where the electrons are coupled to the neutral gas. Assuming the electrons flow radially outward from the nucleus, the exobase is given by the distance at which the scale height of the coma is equal to the mean free path of electrons within the coma. Rosetta was rarely below the electron exobase (defined using electron collision cross-section and a fixed assumed temperature) and most observations of cold electrons occurred throughout the mission outside of this collisional region (Engelhardt et al. 2018; Gilet et al. 2020). For large parts of the mission, the predicted electron exobase was below the surface of the nucleus and no collisional region should have formed (Engelhardt et al. 2018; Gilet et al. 2020). However, even during these times, the cold population of electrons was persistent and observed to have high densities (up to 90 per cent of total density, Gilet et al. 2020; Wattieaux et al. 2020). As such, the assumption of radially outflowing electrons cannot explain the formation of the cold population.

An ambipolar field around the nucleus, set up by the cometary electron pressure gradient, may be a key mechanism in the formation of the cold electron population. Particle-in-Cell (PiC) simulations have shown that an ambipolar field arises around a comet (Deca et al. 2017, 2019; Divin et al. 2020), which can accelerate SW electrons towards the nucleus (Madanian et al. 2016; Galand et al. 2020). This also confines cometary electrons to a region close to the comet (Engelhardt et al. 2018; Sishtla et al. 2019; Divin et al. 2020), where the coma is densest. The trapping increases the column density traced by the electrons and therefore increases the chance that electrons undergo substantial cooling. However, it is not yet known how significant and efficient this process is in the cooling of electrons under low-outgassing conditions.

Gan & Cravens (1990) used a two-stream model to assess cooling of the suprathermal electron population, while applying a fluid approach to assess the temperature of the bulk of the cometary electron population at 1P/Halley. They solved the Boltzmann equation for the suprathermal electrons along a magnetic field line (1D in position space and 2D in velocity). The thermal electrons, which made up the bulk in the inner coma at Halley, are approximated as a fluid with a fixed density and are only treated through the energy equation. This is not applicable to the electron behaviour at a weakly outgassing comet as the bulk of the electrons is made up of several populations with distinct sources and temperatures. Instead, it is necessary to model all three populations of electrons (cold, warm, and hot) kinetically, to capture them in 3D position space and allow for the transition of particles between populations. Madanian et al. (2016) applied a kinetic approach, similar to Gan & Cravens (1990), to comet 67P to look at the formation of the suprathermal electron population, using an ambipolar field derived from a generalized Ohm's law. However, it is not computationally feasible to solve the Boltzmann equation in three spatial dimensions. Without a set of self-consistently generated 3D fields, it is not possible to assess the strength of the cooling because of the coupling between the plasma density and the fields. Self-consistent fields provide the most realistic calculation of the field structure and the depth of the potential well, which are key to the trapping process.

Deca et al. (2017) used a PiC model to treat cometary electrons and ions kinetically, which captures the formation and the dynamics of both the warm and hot electron populations. This also provides a set of electric and magnetic fields, which are self-consistently coupled to the plasma distribution and dynamics, as an output of the simulation. The PiC simulations show the formation of an ambipolar field (Deca et al. 2019; Divin et al. 2020) and subsequent electron trapping (Sishtla et al. 2019). However, they are collisionless simulations and cannot yet be used to model the energy degradation of warm electrons and the subsequent formation of a cold electron population.

In order to self-consistently model the formation and dynamics of all three populations, we require a kinetic approach with a 3D set of electric and magnetic fields. This must be combined with electron-neutral collision processes which can degrade warm electrons in energy, allowing the cold population to form. The electron cooling has been modelled along a magnetic field line by Gan & Cravens (1990) and the hot and warm electron populations have been self-consistently modelled in 3D by Deca et al. (2017). However, a single model has not yet incorporated both the collisional processes and complex electromagnetic environment required to sufficiently capture this cooling process.

Fully kinetic PiC simulations commonly use a reduced ion-electron mass ratio because of tight computational constraints (Bret & Dieckmann 2010). This speeds up numerical simulations while retaining the necessary separation of spatial and temporal scales between the electron and ion dynamics. However, the impact of using a reduced mass ratio (e.g. $m_e = m_p/100$) instead of the realistic mass ratio (with $m_e = 9.11 \times 10^{-31} \text{ kg}$) is a source of concern as it has not been fully investigated. Several studies have considered the effects of using a reduced mass ratio in 2D kinetic simulations for specific applications, such as magnetic reconnection (e.g. Jun & Quan-Ming 2007; Pritchett 2010; Divin et al. 2012), ion beam instabilities (Hong et al. 2012), or drift instabilities (Lavorenti et al. 2021). However, it is computationally unreasonable to perform significant 3D kinetic studies with a physical mass ratio so that this has not been investigated in 3D simulations, although it is key to understanding the limitations of PiC models.

Here, we present a new collisional test-particle model of a weakly outgassing comet. This is the first collisional 3D model of electrons at a comet. A similar Monte Carlo approach was used to model multiple species in Ganymede's ionosphere by Carnielli et al. (2019). The model uses the self-consistently calculated electric and magnetic fields obtained from a fully kinetic, collisionless PiC code (Deca et al. 2017), as an input. This provides realistic fields within the coma, which have previously been used to demonstrate electron trapping by the ambipolar field (Sishtla et al. 2019).

We only consider the weakly outgassing regime as this corresponds to the case of a weakly collisional coma, where collisions do not significantly modify the large-scale plasma distributions or dynamics. This is required to ensure that the fields, which are used as a stationary input to the test-particle model, would not be significantly modified when collisions are included in the simulations.

We combine the fields with electron-neutral collision cross-sections to model elastic and inelastic (ionization, excitation, and dissociation) collisions. This also incorporates electrons from multiple sources: photoionization of the coma, the SW, and secondary electrons from electron-impact ionization. The test-particle model can be used to quantitatively examine the formation of the cold electron population at a weakly outgassing comet.

In Section 2, we present the new test-particle model. The validation of the particle motion in terms of drifts and energy degradation by collisions are outlined in Appendices A and C. We also compare
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic representation of the processes included in the test-particle model. The purple boxes indicate inputs of the simulation. The secondary electron properties are an output of the previous generation, $G_n-1$, but are used as an input for the following generation, $G_n$ ($n > 1$). The blue line indicates the processes that occur within the model, whereas the boxes outside are inputs to, or outputs of, the model.</figcaption></figure>

moments of the electron distributions from our simulation to a PiC simulation driven by the same set of parameters in Section 2.4. Next, we present several early results: in Section 3.1, we look at the trajectory and energy variation of several electrons within the collisional coma; in Section 3.2, we demonstrate the impact of collisions on the wider coma, analysing the electron distribution and some of its moments. In Section 3.3, we compare two test-particle simulations using different electron masses to investigate the effect of changing the ion-electron mass ratio on the electron dynamics in a 3D kinetic model. Finally, in Section 4, we discuss the implications of our results on the current understanding of electron cooling and potential future applications of the model.

## 2 PRESENTATION AND VALIDATION OF THE TEST-PARTICLE MODEL

We have developed a collisional test-particle model of electrons in the cometary environment. The simulation domain extends 2200 km in each direction. This is the domain where the electromagnetic fields are available from the PiC simulations, which we use as an input. The x-axis points along the Sun-comet line towards the nucleus and is aligned with the bulk flow of the SW. At the large heliocentric distances we are considering in this study (~2.3 au), the Parker angle is approximately 90° so it is reasonable to assume that the SW flow and the magnetic field are perpendicular. Therefore, the y-axis points along the direction of the interplanetary magnetic field. The z-axis completes the set. The comet nucleus lies at the centre of the domain in the $\hat{y}$ and $\hat{z}$ directions but is offset from the centre along the x-axis. The origin is defined at the centre of the comet nucleus and the boundaries of the domain along $x$ are at -770 and 1430 km.

We present the mechanics of the model in three parts. In Section 2.1, we discuss the methodology by which the macroparticle

trajectories are calculated, including how electron-neutral collisions have been modelled and implemented. The processes involved in the calculation of the particle trajectories are summarized in Fig. 1. In Section 2.2, we describe how multiple particle trajectories are collated and converted into moments of the particle distributions. Finally, in Section 2.3, the initialization and boundary conditions of the model are presented. This covers the injection of particles into the simulation, as well as the choice of stationary coma and electromagnetic fields. The creation of a particle is dependent on what population it belongs to: photoelectrons (Section 2.3.4), SW electrons (Section 2.3.5), or secondary electrons produced from electron-impact ionization (Section 2.1.3).

Fig. 1 shows the key processes used in the test-particle simulation. The electron macroparticles are injected into the simulation one generation ($G_n$) at a time. That is to say all the photoelectrons (or SW electrons) are created and traced through the simulation first (top-left panel, Fig. 1). Following the completion of the first generation ($G_1$), any further electrons produced through electron-impact ionization are simulated in a second generation, $G_2$. As these secondary electrons can also ionize, further generations of electrons may follow until either no more electrons are produced through collisions or until a maximum number of generations $n_{Max} = 7$ is reached.

### 2.1 Single particle trajectories

#### 2.1.1 Description of the model

The life of a particle within a simulation can be split into three parts: its creation, propagation, and termination. The creation of a particle requires the definition of an initial position and velocity. The choice of position and velocity at its creation is dependent on the population from which the particle is sampled (left-hand panel, Fig. 1) and is discussed in Sections 2.3.4, 2.3.5, and 2.1.3.
The particle is then pushed through the simulation domain as it responds to the electric and magnetic fields at its position (top-middle panel, Fig. 1; Section 2.3.3). The particle may also undergo collisions with the neutral coma (Section 2.1.3), causing it to lose energy or be scattered (right-hand panel, Fig. 1). The macroparticle position and velocity are updated at each step using the Boris integrator (Boris 1970). This continues until the particle leaves the simulation domain at which point it is terminated. This is achieved by passing the outer domain boundaries or by entering the comet nucleus.

The calculation of single-particle trajectories is validated through comparison to standard particle drifts which are known analytically (see Appendix A).

## 2.1.2 Time-step selection

There are a number of conditions on the time-step in order to ensure that the particle trajectories are well resolved. We use a constant time-step to ensure the pusher behaves as expected. In addition, it has to fulfil the following conditions:

(i) $dt < dt_{\text{gyro}} = T_{\text{gyro}}/N_{\text{Cycle}}$. To resolve the gyromotion of the electrons, the time-step size must be small compared to the electron gyroperiod at the particle location, $T_{\text{gyro}} = \frac{2\pi m_e}{qB}$. We choose the parameter $N_{\text{Cycle}} = 20$, which has been shown to closely reproduce the gyromotion of charged particles (Carnielli 2019).

(ii) $dt < dt_{\text{EB}} = \frac{dx_{\text{EB}}}{2|v|}$. It is imperative that the particle does not move further than the resolution of the input electric and magnetic fields, $dx_{\text{EB}}$ (see Section 2.3.3), in a single step. This is necessary to prevent the particle from tunnelling through cells of the field.

(iii) $dt < dt_{\text{Grid}} = \frac{\min(dx_i, dy_j, dz_k)}{2|v|}$, where $dx_i, dy_j$, and $dz_k$ are the dimensions of the grid cell $(i, j$, and $k)$. This ensures that the particle does not jump past cells of the spatial grid which is used to calculate the moments of the particle distributions (see Section 2.3.1).

(iv) $dt < dt_{\text{Eng}} = \frac{d\epsilon_I}{2\epsilon E v}$. We require that the particles do not skip cells in the energy dimension due to acceleration by the electric field. The time-step should ensure that the work done by the fields does not exceed the width of the energy bin, $d\epsilon_I$.

(v) $dt < dt_{\text{prob}} = -\frac{\log(1-p_{\text{max}})}{n(x)\times\sigma_{\text{Tot}}\times|v|}$, where $\sigma_{\text{Tot}}$ is the total electron-neutral collision probability and $n(x)$ is the neutral density at the particle position. This ensures that the probability of an electron-neutral collision is less than $p_{\text{max}}$. We define the parameter $p_{\text{max}} = 0.01$ in our model which makes it very unlikely ($< 0.01$ per cent) that two collisions will occur in the same step. This requirement is discussed further in Section 2.1.3.

## 2.1.3 Electron-impact collisions

A number of electron-neutral collision processes are included in the model. Currently, the model of the coma includes only H$_2$O (see Section 2.3.2); so only collisions between electron and water are accounted for. The addition of other species, such as CO$_2$, would be a straightforward generalization of this model and might be considered in the future. The collision processes include elastic scattering and inelastic collisions (excitation, dissociation, and ionization). The inelastic collision cross-sections for electron impact on water are from Itikawa & Mason (2005).

Inelastic electron-neutral collisions are strongly forward peaked in the differential cross-sections, especially for larger electron energies. Therefore, we approximate all electron scattering into a single elastic process which is outlined in Appendix B. This incorporates the purely elastic scattering as well as scattering from rotational (vibrationally elastic) transitions.

For excitation and dissociation of a water molecule, the energy lost by the incident electron is equal to the threshold energy of the process. The threshold energies of the processes cover several orders of magnitude from 0.004 eV for the rotational transition to 10 eV for the electronic transition (Chutjian, Hall & Trajmar 1975). In the case of ionization, the remaining energy of the incident electron, once the threshold energy is subtracted, is split equally between the outgoing electrons (i.e. two electrons for single ionization, three for double ionization). In the future, it may be useful to sample the energy partition between the outgoing electrons from a distribution (Lummerzheim & Lilensten 1994).

In each step of the simulation, it is determined whether or not a single collision occurs. The probability of each collision process, $i$, occurring in a given time-step $\Delta t$ is

$$P_{\text{coll}, i}(\Delta t) = 1 - \exp[-n(\mathbf{x})\sigma_i v \Delta t]. \quad (1)$$

As the possibility of two collisions occurring within a single step is unlikely (see Section 2.1.2), the total collision probability is given by

$$P_{\text{coll}}(\Delta t) = \sum_{i} P_{i} \approx 1 - \exp[-n(\mathbf{x})\sigma_{\text{Tot}} v \Delta t]. \quad (2)$$

This relation leads to the expected energy loss and number of collisions, although the total probability that a collision occurs in a single step is slightly overestimated. A random number, sampled from a uniform distribution, is used to determine whether a collision occurs and, if so, which type of collision occurs for each step of the particle trajectory.

The energy degradation through electron-neutral collisions is validated by comparison to the continuous slowing down approximation in Appendix C.

### 2.1.3.1 Electron-impact ionization and generation of secondary electrons
When a collision results in electron-impact ionization, all of the outgoing electrons are traced through the simulation. Electron-impact ionization collisions are also forward peaked (Champion, Hanssen & Hervieux 2002) so the incident electron is not scattered, but continues as for any other inelastic collision. In addition, the secondary electrons are also created. The velocity vectors of the new secondary electrons are assigned a random direction, sampled from a uniform distribution. The energies of the new secondary electrons are derived from the energy of the incident electron (see Section 2.1.3). A more refined treatment of energy distribution for secondary electrons could be developed in the future. However, when testing extreme cases of the energy distribution (50 per cent or 90 per cent of the kinetic energy is retained by the incident electron) there was no major impact on the electron energies in the coma.

The positions, energies, and weights (see Sections 2.3.4 and 2.3.5) of the secondary electrons that are produced by the collisions of the previous generation of electrons are recorded. Following the completion of the generation, they are used to create new macroparticles. The velocity vectors of the new macroparticles are assigned a random direction, sampled from a uniform distribution. In order to improve the statistics of the secondary electrons, each ionization spawns many macroparticles ($N_{\text{ioni}} = 10$) so the weight of each secondary electron macroparticle is given by $W_{\text{ioni}}/N_{\text{ioni}}$, where $W_{\text{ioni}}$ is the weight of the ionizing electron. The aggregate gives the required total weight of electrons leaving the ionization site. We consider any secondary electrons produced through electron-neutral collisions up to a maximum of seven generations. After seven generations, the weights of the macroparticles are $10^6$ times smaller than the initial photoelectron or SW electron from which they are descended.
The choice of $N_{\text{ioni}} = 10$ also results in fewer secondary electron macroparticles being produced with each generation. Therefore, particles beyond the seventh generation contribute a negligible amount to the electron moments.

## 2.2 Moments of the electron distributions

The energy distribution function (EDF) of the electrons is calculated on a grid with three spatial dimensions (x, y, and z) and one energy dimension (ε, see Section 2.3.1). The EDF, $f(x_i, y_j, z_k, \epsilon_l)$ (cm⁻³ eV⁻¹) in a given cell (i, j, k, and l) is given by

$$f(x_i, y_j, z_k, \epsilon_l) = \sum_{p, \alpha} \frac{W_p \times dt_{\alpha}}{V_{i,j,k} \times d\epsilon_l}, \qquad (3)$$

where the summation is over all particles, p, that pass through the cell and all points in time, α, that the particle spends in the cell. $W_p$ is the weight of the macroparticle in electrons s⁻¹ (see Section 2.3). $V_{i,j,k}$ is the spatial volume of the cell and $d\epsilon_l$ is the size of the energy bin. The electron density is the integral of the EDF over energy:

$$n_e(x_i, y_j, z_k) = \sum_l f(x_i, y_j, z_k, \epsilon_l) d\epsilon_l. \qquad (4)$$

The components of the bulk electron velocity, $u_s$ for $s = x, y$, and $z$, are calculated only on the spatial grid and are given by

$$u_s(x_i, y_j, z_k) = \frac{1}{n_e(x_i, y_j, z_k) V_{i,j,k}} \sum_{p, \alpha} W_p \times dt_{\alpha} \times v_s. \qquad (5)$$

A similar expression is used to calculate the square velocity of the electrons, which can then be converted into an electron temperature, $T_e$, through equation (7):

$$\vec{v}_s^2(x_i, y_j, z_k) = \frac{1}{n_e(x_i, y_j, z_k) V_{i,j,k}} \sum_{p, \alpha} W_p \times dt_{\alpha} \times v_s^2 \qquad (6)$$

$$\frac{3k_B T_e}{2} = \frac{m_e}{2} \left( \sum_s \vec{v}_s^2(x_i, y_j, z_k) - |\mathbf{u}|^2 \right). \qquad (7)$$

To ensure the moment results are statistically significant, we require that at least 20 particles pass through a given cell to include it in our results.

## 2.3 Initialization of the model

### 2.3.1 Grid

The grid on which the moments are calculated has three Cartesian spatial dimensions centred on the comet nucleus. It is independent of the resolution of the electric and magnetic field. However, the resolution of the fields (see Section 2.3.3) does constrain the resolution of the moment grid, as structures on scales much smaller than the field resolution may not be physical. As such, the minimum spatial size of the grid should not be much smaller than the field resolution. In addition to the spatial grid, a grid in energy space is used to calculate the EDF.

We use two different methods to set the spatial grid size. In the first, we use the same constant grid spacing as the resolution of the electromagnetic fields. This is particularly useful for the direct comparison to the PiC model (see Section 2.4). The second method uses a variable grid spacing across the domain, with higher resolution close to the comet. This allows closer investigation of electron properties within the inner coma, where collisional processes are strongest. Near the comet nucleus, a constant grid spacing of 4 km is used for the 10 closest cells to the comet in each direction.

Beyond this, the length of the cells in each dimension increases exponentially with the distance from the nucleus along each axis. The cell dimensions are only a function of the corresponding component of the position, i.e. $dx_i = dx_i(x_i)$, where $dx_i$ is the side length of a cell with midpoint $x_i$. The same variation is applied in both the y and z directions for $dy_j = dy_j(y_j)$ and $dz_k = dz_k(z_k)$. The grid cells extend up to a maximum distance from the comet of 1000 km, beyond which the moments are not recorded. This allows for a higher spatial resolution close to the nucleus, where collisional processes are strongest. A similar distribution of cell sizes is used for the energy scale. At low energies, the cell size is constant at 0.5 eV for 10 cells before increasing exponentially, until a maximum energy of 300 eV is reached. This energy is significantly larger than the energy of any electrons seen in the model, so it is reasonable to neglect any contribution from higher energies.

### 2.3.2 Neutral coma

We use a simple model of the neutral coma. The modelled coma is spherically symmetric and follows a $1/r^2$ profile in the neutral density, with r the cometocentric distance. We consider only a pure H₂O coma, although the model could be expanded to include CO₂ in the future. In this paper, we set the outgassing rate as $Q = 10^{26}$ s⁻¹ and the outflow velocity of the neutral gas is $u_{\text{gas}} = 1$ km s⁻¹, which are the same parameters used in the PiC simulation (Deca et al. 2017).

### 2.3.3 Electromagnetic fields

The electromagnetic fields are taken from the time-averaged electric and magnetic fields self-consistently obtained from a fully kinetic, collisionless PiC simulation (Deca et al. 2017, 2019). These are a stationary input to the model, so there is no feedback from the simulated particles on the fields. Therefore, it is important that the bulk of the electron population is the same as that modelled in the PiC simulation. As such, we use the same neutral coma and upstream SW conditions as implemented in the PiC model. Other parameters, such as the photoionization rate, are also the same to ensure our approach is consistent.

The fields are evaluated at 288 points in each dimension over a 2200 x 2200 x 2200 km domain. The fields at the upstream SW boundary are $\mathbf{B} = 6$ nT $\hat{\mathbf{j}}$ and $\mathbf{E} = -2.4 \times 10^{-3}$ V m⁻¹ $\hat{\mathbf{z}}$, which result in an $\mathbf{E} \times \mathbf{B}$ drift velocity of $v_{\text{SW}} = 400$ km s⁻¹ downstream.

### 2.3.4 Photoelectrons

Photoelectrons are generated throughout the coma through the ionization of the neutral molecules by EUV photons. We use a constant photoionization rate throughout the coma. This is a reasonable approximation, as we are considering an outgassing rate sufficiently low to have an optically thin coma in the EUV. The photoionization frequency of $v_{\text{ioni}}^{\text{hv}} = 1.32 \times 10^{-7}$ s⁻¹, corresponds to a heliocentric distance of approximately 2.3 au (Huebner & Mukherjee 2015). At 2.3 au, comet 67P had a H₂O production rate of ~3 × 10²⁶ s⁻¹ (Läuter et al. 2018; Biver et al. 2019). This is a few times higher than what we assume, but is still in the low-outgassing regime and therefore applicable to understanding observations from the Rosetta mission.

Many photoelectron macroparticles are generated in each spatial grid cell within 300 km of the comet, but the weight of each macroparticle depends on the size of the grid cell and the neutral
density in the cell. The weight of a photoelectron macroparticle produced within a grid cell is given by

$$W_p = \frac{1}{N_p} \int_{Cell} v_{\text{ioni}}^{hv} n(\mathbf{x}) d^3 \mathbf{x}, \qquad (8)$$

where $N_p$ is the number of macroparticles produced within the cell. This is calculated using a Monte Carlo integration for each grid cell.

The initial position of the macroparticles is distributed uniformly throughout the grid cell. Although this does not reflect the exact distribution of the photoelectron production within one grid cell, the cells are sufficiently small that the $1/r^2$ dependence is reconstructed when looking across many cells. The radial dependence of the production is inherited from the neutral density, as the photoionization rate is constant throughout the coma when the coma is optically thin.

The energy distribution of the photoelectrons is approximated as a Maxwellian distribution with a temperature of 10 eV (Cravens 1987; Galand et al. 2016), which is consistent with the PiC simulation. The photoelectron macroparticles are generated with a random velocity from an isotropic Maxwellian distribution using the Box–Mueller algorithm (Box & Muller 1958).

In Appendix E1 the photoelectron production scheme is validated, using the case of radially outflowing electrons.

### 2.3.5 SW electrons

The second primary source of electrons is the SW. We assume that the SW is undisturbed at the boundaries of the domain, such that the SW density and the velocity distribution are uniform across the faces. We approximate the SW as a plasma of uniform density and a single isotropic temperature, $T_{SW}$. The velocity distribution function, $f_{SW}(\mathbf{x}, \mathbf{v})$, is shifted to the bulk velocity of the SW $u_{SW}$, and is given by

$$f_{\text{SW}}(\mathbf{x}, \mathbf{v}) = n_{\text{SW}}(\mathbf{x}) \left( \frac{m_e}{2\pi k_B T_{\text{SW}}} \right)^{\frac{3}{2}} \exp \left( -\frac{m_e |\mathbf{v} - \mathbf{u}_{\text{SW}}|^2}{2k_B T_{\text{SW}}} \right). \quad (9)$$

This is the distribution which is sampled when generating the initial velocities of the SW macroparticles in the same way as outlined for the photoelectrons (see Section 2.3.4). The total weight attributed to the SW macroparticles is given by the particle flux into the domain. The particle flux into the simulation domain is calculated for five faces of the domain. We do not consider the electrons entering at the downstream face of the domain as the SW is not pristine at this boundary so the properties of the electron population are unknown. Electrons at the downstream boundary are also unlikely to enter the inner coma, as they are dragged downstream by the $E \times B$ field. This has been verified through simulation of many SW particles produced at the downstream boundary (assuming pristine conditions), none of which passed within 300 km of the nucleus. Throughout this study, we have used $n_{SW} = 1 \text{ cm}^{-3}$, $v_{SW} = 400 \text{ km s}^{-1}$, and $T_{SW} = 10 \text{ eV}$ as the parameters of the SW, which is in agreement with the SW parameters used by Deca et al. (2017).

The weight of a SW macroparticle, $W$, is given by

$$W = \frac{A_{\text{Face}} \times n_{\text{SW,in}} \times v_{\perp}}{N_{\text{Face}}}, \qquad (10)$$

where $N_{\text{Face}}$ is the number of macroparticles created at each face of the domain and $A_{\text{Face}}$ is the area of the face. $v_{\perp}$ is the velocity of the macroparticle perpendicular to the face through which it enters the domain. $n_{\text{SW,in}}$ is the density of SW electrons at each boundary that will enter the domain, i.e. the zeroth-order moment of particles

travelling into the box. The SW weighting and velocity distributions are validated in quiet SW conditions in Appendix E2.

## 2.4 Comparison to collisionless PiC models

### 2.4.1 Photoelectrons

The test-particle model is further validated by comparison of photoelectron densities from collisionless test-particle simulations (Figs 2a, b, d, and e) and PiC simulations (Figs 2c and f). The PiC code is the one used for defining the electromagnetic field in the test-particle code (see Section 2.3.3). We first consider a test-particle simulation using the same non-physical electron mass ($m_e = m_p/100$; see Figs 2b and e) used by the PiC simulations. In addition to the adjustment of the electron mass, the production of photoelectrons is capped close to the comet ($r < 11$ km) in the PiC model to improve performance. To ensure consistency between approaches for this validation, we have also capped the production rate in the test-particle simulation.

When the same (artificial) electron mass is used, the test-particle and PiC simulations agree very closely in both shape and magnitude in the wider coma. A broad near-tail is seen extending 400 km downstream of the nucleus in both the test-particle and PiC simulations (see Figs 2e and f). The near-tail exhibits the same shape in all three cases and there is close agreement in density between approaches (<15 per cent deviation from the PiC density).

A narrow tailward offshoot is seen extending downstream from the inner coma at $z = 100$ km in the $xz$ plane (see Figs 2a–c). The near-tail in the test-particle simulation (Fig. 2b) has the same width as in the PiC simulation (see Fig. 2c) but the density is 25 per cent larger in the test-particle simulation up to 100 km downstream of the nucleus. Beyond 100 km downstream, the electron density is 80 per cent larger in the test-particle simulation. The enhanced density seen in Fig. 2b is a result of the tailward branch being more tightly confined along $y$ in the test-particle simulation than in the PiC simulation. When the width of the branch (up to ±200 km along $y$) is accounted for, the total electron density in the region is similar in both the test-particle and PiC simulations. The difference in the spatial extent along $y$ may be caused by the fact that stationary fields are used in the test-particle model (by construction) instead of dynamical fields in the PiC simulation. In the inner coma ($r < 50$ km), the photoelectron density is up to 3.5 times larger in the test-particle simulation when compared to the PiC simulation. This is caused by differences in the photoelectron production close to the nucleus between the two models. The PiC simulation uses a larger cometoradius than the test-particle model ($r = 1.7$ km) which reduces the number of electrons produced within the coma. When we increase the radius of the nucleus, the electron densities in the inner coma decrease to within a factor of 2 of the PiC densities. The electron density in the inner coma is particularly sensitive to the production close to the nucleus. This is where the photoelectron production peaks and is also where the potential well, formed by the ambipolar field, is at its deepest. Consequently, electrons produced very close to the nucleus spend more time trapped in the region, which increases the electron density. This enables us to point out the importance of plasma production in the near-nucleus region, which should be taken into account in future kinetic (full kinetic or hybrid) cometary plasma simulations.

The test-particle model is significantly less computationally intensive than the PiC model, which allows higher ion-electron mass ratios to be explored. Actually, large-scale, fully kinetic, 3D PiC simulations using a realistic ion-electron mass ratio are still out of
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Photoelectron density in the xz (a–c, at y = 0 km) and xy (d–f, at z = 0 km) planes. (a) and (d) show the electron density from a test-particle simulation with electron mass $m_e = 9.11 \times 10^{-31}$ kg. (b) and (e) show the densities from a test-particle simulation with an increased electron mass $m_e = m_p/100$. (c) and (f) show the photoelectron density from a PiC simulation which uses the increased electron mass.</figcaption></figure>

reach of state-of-the-art computational facilities. Therefore, reduced ion-to-electron mass ratios are used in such models, even if the impact of using such reduced mass ratios is not assessed. To do so, we hereafter perform a test-particle simulation of photoelectrons using the realistic electron mass (see Fig. 2a and d). The electron density is largely unchanged throughout the coma from the case using an increased electron mass. There is a notable deviation in the narrow tailward region, seen extending from the coma at z = 100 km in the xz plane. When using a smaller electron mass, the near-tail is lower in density (by 30 per cent with respect to the PiC simulation) and has a thinner profile in the xz plane. The lower density and narrower projection in the mass case is compensated by a wider extent along y (not shown). Using the realistic electron mass, the width of the near-tail along y is the same as the width in the PiC simulation. A low-density region is seen below the tailward branch, due the thinner projection in the xz plane (see Fig. 2a). This is a consequence of the disparity in velocity scales between the high- and low-mass simulations, which drives differences in the electron transport.

The test-particle and PiC simulations are largely in good agreement when the same electron mass ($m_e = m_p/100$) is used. The discrepancy in the magnitude in the inner coma is driven by differences in the production scheme close to the nucleus between the two simulations. The difference in the near-tail is caused by a difference in the transport from the inner coma between the simulations, which may be caused by the stationary nature of the fields used in the test-particle model. When the realistic electron mass is used, the wider coma is largely unchanged from the case with the increased mass. Differences in the shape of the tailward region are attributed to the change in transport from the inner coma caused by the lower electron mass.

### 2.4.2 SW electrons

The final validation of the SW implementation is via comparison with the electron densities from a PiC simulation (see Figs 3c and

f). The PiC simulation is the same used to calculate the electric and magnetic fields (see Section 2.3.3). We use the same electron masses considered in the photoelectron validation described in Section 2.4.1, that is $9.11 \times 10^{-31}$ and $1.67 \times 10^{-29}$ kg. When using the enhanced electron mass, which is also used in the PiC simulation, the test-particle simulations (see Figs 3b and e) and PiC simulations (see Figs 3c and f) agree closely. The shape of the electron cloud is very similar in both the xy and xz planes. Both show a region of high density around the nucleus, with a roughly spherical shape, which leads to a narrow, near-tail structure directly anti-sunward in y and moving towards -z for increasing x. A high-density region of the near-tail is also observed around (300,0, -200) km in both the test-particle and PiC simulations (see Figs 3b and c).

However, the magnitude of the electron density is approximately 40 per cent larger close to the nucleus in the test-particle simulation. The SW electron density peaks at $7 \text{ cm}^{-3}$ in the test-particle model, as opposed to $5 \text{ cm}^{-3}$ in the PiC simulation. This may be caused by the stationary nature of the fields in the test-particle model, whereas they are time dependent in the PiC simulation.

When the physical electron mass is used in the test-particle simulation (see Figs 3a and d) there are some noticeable differences from the higher mass case (see Figs 3b and e). The shape of the electron cloud is less spherical and peaks at lower densities ($4.5 \text{ cm}^{-3}$ when using $m_e = 9.11 \times 10^{-31}$ kg). The narrow tailward branch seen in the higher mass simulations is less well defined, although in both cases a region depleted of electrons is seen at z = 50 km from the comet position and further downstream.

The major driver of the differences in the high- and low-electron mass cases is the set of conditions at the upstream boundary. The SW electrons entering the domain from upstream are sampled from a 10 eV Maxwellian with a bulk velocity of 400 km s$^{-1}$ (see Section 2.3.5). Increasing the mass of the electron causes a reduction in the thermal velocity. As the thermal velocity is diminished with respect to the bulk velocity, a larger proportion of the SW electrons have a positive
<figure><img src="image_5.png" /><figcaption>**Figure 3.** SW electron density in the xz (a–c, at y = 0 km) and xy (d–f, at z = 0 km) planes. (a) and (d) show the electron density from a test-particle simulation with electron mass $m_e = 9.11 \times 10^{-31}$ kg. (b) and (e) show the results from a test-particle simulation with an increased electron mass $m_e = m_p/100$. (c) and (f) show the SW electron density from a PiC simulation which uses the increased electron mass. The electron densities have been smoothed with a 3x3x3 sliding window.</figcaption></figure>

$x$ velocity and therefore more electrons enter from upstream. When using the lower, physical mass, only 62 per cent of SW electrons at the upstream boundary enter the domain, whereas this reaches 90 per cent in the high-mass case. There is also a noticeable difference in the energy distributions of electrons entering at the domain boundaries, due to the kinetic energy associated with the bulk SW flow (0.46 and 8.4 eV in the low- and high-mass cases).

Although there are some differences between the test-particle and PiC simulations, the electron behaviour is in agreement on large scales. Therefore, the electric and magnetic fields from the PiC simulation are a good set of fields which capture all the relevant plasma processes in three dimensions while being largely consistent with the densities calculated in the test-particle simulations.

# 3 RESULTS

## 3.1 Single particle trajectories

We present two typical examples of electron trajectories which serve to demonstrate many aspects of the model: the production of SW and photoelectrons; elastic and inelastic collisions; and electron-impact ionization. The first trajectory is a photoelectron produced within the coma (see Fig. 4a). The electron is produced at (−49.5, −2.7, and −72.6) km with an energy of 20.56 eV, with 12.3 eV of this perpendicular to the local magnetic field (see Fig. 4b). As the initial velocity along the magnetic field line (approximately aligned with the y-axis) is towards the nucleus, the electron accelerates as it falls deeper into the potential well created by the ambipolar field (Deca et al. 2017). This initial acceleration increases its parallel energy (yellow, Fig. 4b) to 12 eV. The parallel energy then decreases to 0 eV as the electron climbs the potential barrier along the magnetic field

line. The electron continues to be trapped within the potential well and the parallel energy oscillates between 0 eV and a maximum up to 24 eV. The perpendicular energy (orange, Fig. 4b) is seen as a band due to the fast gyromotion of the electron in the presence of a perpendicular electric field.

Theoretical drift velocities have been calculated at each point along the trajectory using the energy of the particle and the fields at its location (see Appendix A). The gyroaveraged velocity perpendicular to the magnetic field follows closely the expected behaviour from particle drifts. The gyroaveraged motion of the electron is largely attributed to an $E \times B$ drift with a magnitude up to $6 \times 10^5$ m s$^{-1}$. Within the inner coma, the particle also undergoes curvature drift of up to $6 \times 10^4$ m s$^{-1}$ due to the draping of the magnetic field lines around the cometary ionosphere. Once the particle leaves the inner coma at $t = 1.2$ s, the velocity closely follows the theoretical $E \times B$ drift until the particle leaves the domain at the downstream boundary.

The photoelectron undergoes two collisions with the neutral coma which can be identified by the jumps in the perpendicular energy (orange, Fig. 4b) at 0.54 and 0.67 s. The first collision to occur is elastic so the total energy of the particle is unchanged but the velocity is rotated by a scattering angle of 122°. The parallel energy decreases from 15.04 to 7.42 eV, while the perpendicular energy has a corresponding increase from 12.6 to 20.22 eV. The second collision to occur is an electronic transition ($E_{Th} = 10$ eV), which degrades the electron in energy from 35.9 to 25.9 eV. The velocity vector is not scattered in this case, so the changes in the parallel and perpendicular energies are caused only by the energy loss.

The second example trajectory, as shown in Fig. 5, is a SW electron that passes through the inner coma and causes ionization of a water molecule through electron impact. The trajectory of the secondary electron is shifted by 200 km along y for visibility. The electron enters
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) Trajectory of a photoelectron through the inner coma. The comet nucleus (black) has been scaled up for visibility. (b) Variation of the particle energy as it travels through the cometary environment. The energies parallel ($\epsilon_\parallel$) and perpendicular ($\epsilon_\perp$) to the local magnetic field are shown in yellow and red, respectively. The colour bar at the top maps the line colour in panel a to the total energy ($\epsilon_{tot}$; variation in time shown in b).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Trajectory of a SW electron (bottom panel) that causes ionization within the coma. The secondary electron trajectory (top panel) is plotted with a shift of 200 km in the y direction.</figcaption></figure>

the simulation at the upstream x boundary with energy of 27.6 eV. The particle drifts downstream towards the cometary nucleus and is accelerated from 30 to 110 eV as it falls into the potential well around the nucleus. The acceleration of the electron is primarily seen as an increase in the perpendicular energy (from 25 to 100 eV).

The perpendicular energy gained by the SW electron results from a combination of kinetic processes. First, the electron undergoes curvature drift as it nears the nucleus, which has a large component anti-parallel to the electric field ($v_{curv} \cdot \hat{E} \sim -2 \times 10^4$ m s$^{-1}$), which accelerates the electron. Alongside this acceleration, the strength of the magnetic field increases from 6 to 13.6 nT as the SW electron nears the nucleus. Under adiabatic conditions, the perpendicular energy scales linearly with the magnetic field strength. Although the electron motion is not adiabatic in the presence of a potential well, it is still useful to consider the impact of this relationship. The combination of the increase in magnetic field and the curvature drift account for an increase of 35 eV in the perpendicular energy. During the acceleration, the electron gyroradius is between 2 and

3 km, which is several times smaller than the resolution of the fields (7.67 km). The remaining perpendicular acceleration (~40 eV) may result from large gyroradius effects over ~400 gyrations.

This contrasts with the photoelectron in Fig. 4, which undergoes very little acceleration perpendicular to the field. The magnetic field strength at the photoelectron does not vary substantially in the inner coma and the magnetic moment is approximately constant, except when collisions occur. The electron drifts ($E \times B$ and curvature) are closely aligned with the $E \times B$ direction, resulting in movement along an equipotential. The energy variation is dominated by the motion parallel to the magnetic field, which leads to the oscillation in the potential well.

The SW electron passes within 3 km of the nucleus as it oscillates within the ambipolar electric potential well. The particle is confined as the parallel energy (<25 eV) is smaller than the potential barrier along the magnetic field line. In the densest region of the coma, the SW electron collides several times, one of which generates an $OH^+$ and a secondary electron through ionization of a water
<figure><img src="image_8.png" /><figcaption>**Figure 6.** EDFs of cometary photoelectrons and their secondaries along the y-axis going through the nucleus (0,0,0) (a) the non-collisional case and (b) the collisional simulation. Cells passed through by fewer than 30 macroparticles are shown in white.</figcaption></figure>

molecule. The incident electron energy drops from 104.4 to 43.14 eV and a secondary electron is produced with an initial energy of 43.14 eV.

The secondary electron is also confined to the inner coma by the ambipolar electric potential well, where it also undergoes several collisions. The secondary electron goes on to ionize a water molecule producing a tertiary electron (not shown) and is degraded in energy from 43.9 to 12.9 eV. The first secondary electron continues to be trapped in the well after the ionization as it drifts in the +z direction until it flows tailward. The behaviour and dynamics of the secondary electrons are similar to those exhibited by photoelectrons (e.g. Fig. 4).

## 3.2 Impact of collisions on electron populations in the coma

To examine the impact of electron-neutral collisions on the electrons in the coma, we compare two test-particle simulations: one without electron-neutral collisions and another where they are included.

The electron-distribution function at the comet position is shown for a slice along y in Fig. 6 for both the non-collisional and collisional simulations. Away from the inner coma, the shape and magnitude are almost identical in both simulations. However, there is clearly a significant difference close to the comet nucleus. In the non-collisional simulation, a population of high-energy electrons (up to 160 eV) is seen within 25 km of the nucleus. These are photoelectrons generated further away from the nucleus and then accelerated by the ambipolar field. As they are produced away from the nucleus, they behave quite similarly to SW electrons, as shown in Fig. 5. In the collisional case, the same acceleration processes are present, but the high-energy population is efficiently quenched by electron-neutral collisions. Electron-impact ionization is very efficient at cooling the highest energy electrons, due to the large associated energy loss and cross-sections peaking around 100 eV. The high-energy tail (>100 eV) is three orders of magnitude weaker when collisions are included but the density of lowest energy electrons (<8 eV) is enhanced by a factor of 3.

Between 20 and 40 km from the nucleus, there is still a reduction in the higher energy electrons from 20 to 80 eV by up to 30 per cent and a slight enhancement in electrons below 10 eV when collisions are considered. Beyond 40 km from the nucleus, there is little impact of collisions on the electron distribution function along the y direction. The impact of electron-neutral collisions that occur in the inner part of the coma, can also be observed in more distant regions of the coma, such as in the tailward direction. Indeed, in some regions, electrons that have suffered collisions in the inner coma are transported from the inner coma to the outer coma, where no significant local collisionality is expected anymore. This emphasizes how local cooling, associated with transport processes, can impact the coma at large scales.

In the collisional simulation, the electron temperature (see Fig. 7a) drops as low as 2 eV within 10 km of the nucleus where collisional processes are their strongest. The near-tail region extending from the inner coma also has electron temperatures of ~4 eV as the low-energy electrons are transported away from the coma. This low-temperature region is also seen in the non-collisional case, but the electron density is much higher when electron-neutral collisions are included. When collisions are added, a significant decrease in the electron temperature is observed within 60 km of the nucleus in +z and then stretching towards the tail (see Fig. 7b). When collisions are neglected, the highest energy electrons seen near the nucleus (see Fig. 6a) are transported into this region and then tailward. As such, the temperature in this region peaks at 17 eV. However, as the high-energy electrons are efficiently cooled by collisions with water in the coma, the temperature is much lower here in the collisional case at 6 eV.

The temperature is largely unchanged by collisions for $z < -20$ km (see Fig. 7) as the coma is sparse and there is little transport from the inner coma. The high-electron temperature (~19 eV) seen below the coolest region of the coma is driven by the potential well. The increase in temperature is primarily along the y-axis and therefore can be seen as an increase in the parallel energy. The high temperature
<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a) Photoelectron temperature in the xz plane from a collisional simulation. (b) Ratio of electron temperatures from a collisional and non-collisional simulation. Cells passed through by fewer than 200 macroparticles are shown in white.</figcaption></figure>

originates from photoelectrons, produced far from the comet, falling into the potential well and being accelerated along the magnetic field lines.

Despite these large differences in temperature, the large-scale behaviour of the cometary electrons is very similar in both the collisional and non-collisional simulations. The electron density exhibits the same shape in both simulations and is of the same order of magnitude throughout the coma (see Figs 2a and d for collisionless case). Away from the inner coma ($r > 40$ km), the plasma density is within 5 per cent of the non-collisional case. Close to the nucleus, the density does deviate by up to 40 per cent. This is reflective of a slight shift in the trajectories preferred by electrons from the innermost coma, arising from the disparity in the typical electron energy in different regions of the coma. The electron density is enhanced in the collisional case, compared with the collisionless simulation, in a region where the temperature is low. This region extends upstream from the inner coma as seen in Fig. 7a. However, there is a decrease in density of similar size just downstream of the low-temperature region, where higher energy electrons make up the population. These two paths converge between 60 and 80 km from the comet in the +z direction. The region of the coma magnetically connected to the nucleus also sees a 40 per cent drop in electron density, in the collisional case, between ±10 and ±60 km along y (not shown). The electrons close to the nucleus are likely to undergo collisions and therefore become degraded in parallel energy. As a result, the electrons are more tightly confined by the ambipolar field and cannot travel as far along the magnetic field lines. The validity of the electric and magnetic fields in a collisional simulation is discussed in Section 4.

### 3.3 Comparison of physical and reduced ion-electron mass ratio

A reduced ion-electron mass ratio ($m_p/m_e = 100$) is often used in PiC simulations due to computational constraints. There have been several studies investigating the impact of the ion-electron mass ratio (e.g. Jun & Quan-Ming 2007; Hong et al. 2012) but not in a 3D simulation. Here we consider two 3D collisionless test-particle simulations with a realistic ($m_e = 9.11 \times 10^{-31}$ kg; see Fig. 8a) and an increased ($m_e = 1.67 \times 10^{-29}$ kg; see Fig. 8b) electron mass. The simulations are applied to SW electrons evolving in electromagnetic fields, as described in Section 2.3.3.

The energy distributions in Fig. 8 are similar between simulations, using the two electron masses, near the upstream boundary (at $x < -300$ km) as there has been little perturbation from the 10 eV Maxwellian initiated at the domain boundaries (see Section 2.3.5). Close to the nucleus, the maximum electron energy increases to ~250 eV which is driven by the potential well which exceeds 200 eV in depth in this region. Notably, there is a significant difference at lower energies (<70 eV) between simulations with the two electron masses. When using the realistic electron mass, few SW electrons occupy energies below 70 eV around $x = 0$ km. However, SW electrons are abundant in this energy range when a larger electron mass is used. The electrons move through the same potential well in both cases which is reflected by the upper energy bound, so the disparity at low energy cannot be attributed to this factor. The distinction is due to the difference in the gyromotion when the particle mass is changed. During a single gyration, the energy of the electron varies due to the presence of an electric field, which leads to phenomena such as the $E \times B$ drift. The size of the gyroradius drives the scale of the energy variation given a set of fields and increases with the electron mass (by a factor of $\sqrt{m_{\text{heavy}}/m_{\text{real}}} = 4.3$ for the same perpendicular energy).

We consider uniform magnetic and electric fields of $B = 6$ nT $\hat{y}$ and $E = -2.4 \times 10^{-3}$ Vm$^{-1}\hat{z}$, which are typical values upstream in the SW. The electrons are pushed for several gyroperiods (see Section 2.1.1) for a range of initial perpendicular energies, where $v_{\text{init}} = v_{\text{init},\hat{z}}$. The maximum energy variation is shown in Fig. 9 for both the heavy (orange) and realistic (blue) electron masses. An electron with an initial perpendicular energy of 100 eV loses only 25.1 eV for the realistic electron mass, whereas the heavier electron decreases in energy by 82.1 eV. From 50 to 100 eV initially, the heavy electrons lose between 72 per cent and 97 per cent of their perpendicular energy in a given gyration under SW conditions.

As such, a heavy electron may contribute to the distribution function at 100 eV during one part of its gyration and then at 30 eV half a gyration later. Through this mechanism the electrons, which have been accelerated in the potential well, can still appear at low energies in the inner coma. In the realistic mass case, the smaller gyroradius gives a smaller range of energies that can be reached by a single electron. Consequently, the accelerated electrons, with a realistic mass, cannot reach low energies (<70 eV) close to the nucleus, which explains the absence of low-energy electrons near $x = 0$ km in Fig. 8a.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** EDF of SW electrons using (a) $m_e = 9.11 \times 10^{-31}$ kg (b) $m_e = 1.67 \times 10^{-29}$ kg. Both test-particle simulations are non-collisional and shown for a slice along the x-axis through the comet nucleus.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Energy variation of an electron during a gyration given $B = 6$ nT $\hat{y}$ and $E = -2.4 \times 10^{-3}$ Vm$^{-1}\hat{z}$. Above 8.35 eV, the variation is a reduction in energy for the heavy electron mass, (i.e. falls from 100 to 17.9 eV), but below 8.35 eV the variation sees the electron increase in energy (i.e. 1 eV rises to 22.8 eV). For the physical electron mass, the boundary is at 0.46 eV.</figcaption></figure>

## 4 DISCUSSION AND CONCLUSION

We have presented the first 3D collisional model of a cometary ionosphere. The test-particle model incorporates electron-neutral collision processes into the cometary environment, which are seen to have a significant effect even at a very low-outgassing rate ($Q = 10^{26}$ s$^{-1}$). We validated the test-particle model with comparison to simple analytical models. Theoretical expressions for particle drifts and the energy loss from electron-neutral collisions are replicated well by the model (see Appendices A and C). The calculation of moments has been corroborated for both the SW and photoelectrons in simple cases (Sections 2.3.4 and 2.3.5).

In addition, we have compared electron moments from a collisionless simulation of our test-particle model and from the PiC model (see Section 2.4) which is used to generate the input electric and magnetic fields. When using the same (reduced) ion-electron mass ratio ($m_e = m_p/100$), the properties of the electron populations in the

PiC and collisionless test-particle simulations are similar throughout the coma for both the SW and photoelectron populations as expected. Differences in the photoelectron densities close to the nucleus are caused by differences in the production close to the nucleus. There are some smaller differences between the models, such as the high-energy tail of electrons observed close to the nucleus. This may result from the different model assumptions and resolutions as PiC simulations are less able to capture the extended tails of a population. However, the models generate the similar bulk electron properties throughout the coma, which further validates the new test-particle model. The large-scale structures in the coma are also consistent between the collisional and non-collisional test-particle simulations.

When the realistic electron mass is used, the difference in scale of the velocity causes some changes in the electron density. For the photoelectron population, this is seen as a shift in the narrow near-tail structures (see Figs 2a and b) due to the difference in transport from the inner coma. The change in shape and magnitude of the SW electrons is caused by the boundary conditions upstream in the SW.

Despite some differences in the electron density between the test-particle and PiC simulations, the self-consistently calculated electric and magnetic fields from the PiC simulation are the most realistic set of fields available to examine the effectiveness of the trapping process. While the electromagnetic fields generated by the PiC and used as input of the test-particle model are not consistent with the electron densities calculated with the test-particle model, the PiC and test-particle electron densities have similar large-scale structures and are of similar order. This agreement is close enough to provide full confidence in the outcomes of this study, that is, close enough to demonstrate the efficiency of the trapping in the potential and the impact it has on collisional processes.

That said, it is worth looking at the potential feedback of the test-particle outcome on the electric and magnetic fields. When collisions are included, the electron density increases and temperature decreases close to the nucleus. These changes would set up a feedback on the fields in the inner coma in a self-consistent simulation. We cannot directly compute the total electric and magnetic fields in each test particle simulation as we only model the electron population. We can only make an assessment of the impact collisional processes have on the electric and magnetic fields, and to propose how these would feed back in a self-consistent model. The magnetic field is
primarily driven by the SW, which is compressed as it interacts with the cometary ionosphere. The higher electron density in the inner coma would lead to increased pile-up of the SW and, consequently, some increase in the magnetic field strength. However, the magnetic field quickly saturates during pile-up (Goetz et al. 2017) so the resultant change in field strength is not expected to be significant.

The electric field within the coma is impacted by collisions in several ways. First, the increase in electron density without change in the ion population violates quasi-neutrality. However, with the introduction of collisions there would also be some enhancement of the ion density in the inner coma through electron-impact ionization, which would mitigate the charge imbalance. Beyond this, it is likely that quasi-neutrality would be restored by transporting any excess of cold electrons away from the inner coma. Assuming quasi-neutrality is fulfilled, the electric field is assessed using the generalized Ohm's law for a collisionless and a collisional test-particle simulation. The ambipolar field is derived from simulated densities and temperatures, such as those shown in Figs 2a and 7a. The electron contribution to the Hall field is computed from similar measurements of the electron density and bulk velocity. Typically, ion speeds are much slower than electron speeds, but the bulk velocities of both species are primarily driven by the field conditions so are likely of similar magnitude. The ion properties are not calculated in the test-particle model, we can only compare the electron components of the Hall field in collisional and collisionless simulations. This is still indicative of regions where the bulk velocity is substantially impacted by collisional processes. As we do not expect significant differences in the magnetic field, we have assumed it is unchanged from the collisionless PiC simulation for calculation of the Hall field.

Away from the inner coma ($r > 30$ km), we found that the ambipolar and Hall fields show no significant difference between the collisional and collisionless cases. The ratios in field strengths through the wider coma ($30 < r < 300$ km) is close to 1 for both the ambipolar field ($0.9 \pm 0.26$) and the electron component of the Hall field ($1.00 \pm 0.07$). However, close to the nucleus ($r < 30$ km) both the ambipolar and electron Hall fields become weaker with the introduction of electron-neutral collisions. Within 30 km of the nucleus, the ambipolar (electron Hall) field calculated in the collisional simulation is reduced by a factor $0.28 \pm 0.33$ ($0.41 \pm 0.38$) with respect to the collisionless simulation. This is driven by the reduction in the electron temperature (see Fig. 7) and bulk velocity (not shown) in the collisional region of the coma.

The ambipolar field sets up the potential well around the inner coma and is the cause of electron trapping. As a result of the weakened ambipolar field, the depth of the parallel potential well also decreases when collisions are added but only by a factor of 0.85 near the nucleus. The potential barrier is also built up over a much larger distance ($y \approx \pm 150$ km) than spanned by the collisional region. Consequently, the electron trapping is not greatly diminished despite the substantial weakening of the fields close to the nucleus. As the trapping weakens with collisions, a negative feedback is set up, suggesting a self-consistent model may lie between the collisionless case and the collisional cases presented in this study.

Single particle trajectories demonstrate several key processes within the collisional coma. First, the confinement of both cometary and acceleration of SW electrons by an ambipolar electric field (Deca et al. 2017, 2019; Divin et al. 2020) is observed (see Figs 4a and 5) in line with the proposal by Madanian et al. (2016). Trapping of electrons in the fields generated by Deca et al. (2017) has also been demonstrated by Sishtla et al. (2019) but only in the collisionless case. The parallel energy of the electron oscillates between a maximum and 0 eV as the electron moves between the centre of the well and

an extreme along the magnetic field lines. This electron trapping process greatly increases the time spent by the electrons within the inner coma, and consequently greatly increases the electron depth built up along the trajectory. This greatly increases likelihood of an electron undergoing a collision even when the traditionally defined electron exobase (Lemaire & Scherer 1974; Mandt et al. 2016) is well below the nucleus surface. The electrons within the coma are not trapped perpetually within the inner coma as the drifts associated with the local fields drive the electrons towards the comet tail.

As well as electron trapping, we have demonstrated electron-impact processes (excitation, ionization) within the coma, which have not been incorporated into other 3D models of the cometary environment. This contributes to the coma both as a source of electrons and as a method of energy loss. The energy loss is particularly key for the highest energy electrons within the coma (see Section 3.2). As the electron energy (less $E_{th}$) is distributed between the multiple electrons leaving an ionization site, the high-energy electron tail is efficiently quenched. Currently, we use an extreme case of ionization with the energy equally distributed between the incident and secondary electrons. To test the other extreme, we ran a simulation with 90 per cent of the energy retained by the incident electron and the high-energy tail of the electrons was still efficiently quenched.

The impact of collisions on the wider electron population is also addressed. Even at the low-outgassing rate of $Q = 10^{26}$ s$^{-1}$, electron-neutral collisions have a significant impact on the cometary electron population within the coma. The highest energy electrons are efficiently quenched in the collisional case and the density of low-energy ($<5$ eV) electrons increases close to the nucleus. A noticeable reduction in the electron temperature is also observed close to the nucleus and propagates into the near-tail region.

The low-outgassing rates considered in this paper are a factor of 5 smaller than the lowest outgassing case considered by Engelhardt et al. (2018) for radial outflow. Negligible cooling was expected by Engelhardt et al. (2018), given the radial trajectories and the exclusion of electron excitations and electron-impact ionization. The trapping of the cometary electrons by the ambipolar field boosts the efficiency of the collisional processes and causes a noticeable difference in the cometary electron behaviour even in a very sparse coma.

This combination of trapping and electron-neutral collisions is a good candidate to explain the cold electron observations throughout the Rosetta mission when the electron exobase was below the surface of the nucleus (Gilet et al. 2020). While the test-particle model outputs are not directly comparable to observations from the Rosetta mission, it does give a useful insight into the collisionality of a weakly outgassing coma and a quantitative assessment of the impact of electron trapping. The acceleration of a SW electron and its collisionality as it enters the coma also lend support to the categorization of far-ultraviolet emissions from the coma of comet 67P as aurora (Galand et al. 2020; Stephenson et al. 2021).

At lower outgassing rates, the simulation would tend towards the collisionless regime. No significant electron cooling would be observed, although electron-impact ionization would remain a significant source of electrons. In extensions to higher outgassing rates, we would expect cometary electrons to exhibit similar behaviour. Significant cooling of the electrons should persist and become more efficient, as well as extending further from the nucleus. The importance of electron-impact ionization may decrease as photoelectrons become more dominant in the coma, especially when adjusted to a lower heliocentric distance for consistency. At high-outgassing rates, the model could be pushed to the limit of its domain of validity (not done here), where the cold electrons form the bulk
of the electron population. The large cold population would greatly affect the structure of the potential well around the nucleus.

A more detailed look at the cooling of electrons is beyond the scope of this paper but will be addressed in future studies. We also hope to apply this model to other electron populations within the coma going forward. In addition to these studies, this model has potential to feed back into the PiC models with regard to the incorporation of collisional processes, where significant simplifications of the collision mechanics would be required.

We have also provided a framework to investigate the impact of using a reduced mass ratio in full kinetic numerical simulation, by using a test-particle model. We have applied this method to cometary plasmas and performed a comparison of two test-particle simulations using a reduced ($m_e = m_p/100$) and a realistic ion-electron mass ratio. Electrons with a realistic mass are more restricted in energy space, most probably due to the smaller gyroradius, which affects the temperature of the electron population. The disparity in the energy variation may also have implications on the transport and density in the coma, which would feed back into the fields.

## ACKNOWLEDGEMENTS

Work at Imperial College London was supported by the Science and Technology Facilities Council (STFC) of the UK under grant ST/N000692/1 and ST/S505432/1. JD gratefully acknowledges support from NASA's Rosetta Data Analysis Program, grant no. 80NSSC19K1305. This work was supported in part by NASA's Solar System Exploration Research Virtual Institute (SSERVI): Institute for Modeling Plasmas, Atmosphere, and Cosmic Dust (IMPACT), and the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. We acknowledge Centre national d'études spatiales (CNES) fundings for the work performed at the Centre national de la recherche scientifique (CNRS). We would also like to acknowledge Arnaud Beth for his useful discussions on kinetic modelling.

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

## REFERENCES

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

* Boris J. P., Shanny R. A., 1970, Proceedings of the Fourth Conference on Numerical Simulation of Plasmas. Naval Research Laboratory, Washington, p. 3

* Box G. E. P., Muller M. E., 1958, Ann. Math. Stat., 29, 610

* Bret A., Dieckmann M. E., 2010, Phys. Plasmas, 17, 032109

* Broiles T. W. et al., 2016, J. Geophys. Res. Space Phys., 121, 7407

* Burch I. L., Goldstein R., Cravens T. E., Gibson W. C., Lundin R. N., Pollock C. J., Winningham J. D., Young D. T., 2007, Space Sci. Rev., 128, 697

* Carnielli G, Galand M., Leblanc F., Leclercq L., Modolo R., Beth A., Huybrighs H.L.F., Jia X., 2019, Icarus, 330, 42

* Carnielli G., 2019, PhD Thesis, Imperial College London

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

* Champion C., Hanssen J., Hervieux P. A., 2002, Phys. Rev. A, 65, 022710

* Cho H., Park Y., Tanaka H., Buckman S., 2004, J. Phys. B: At. Mol. Opt. Phys., 37, 625

* Chutjian A., Hall R. I., Trajmar S., 1975, J. Chem. Phys., 63, 892

* Clark G. et al., 2015, A&A, 583, A24

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

* 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

* Divin A., Deca J., Eriksson A., Henri P., Lapenta G., Olshevsky V., Markidis S., 2020, ApJ, 889, L33

* Divin A., Lapenta G., Markidis S., Semenov V. S., Erkaev N. V., Korovin- skiy D. B., Biernat H. K., 2012, J. Geophys. Res. Space Phys., 117, A06217

* Edberg N. J. T. et al., 2015, Geophys. Res. Lett., 42, 4263

* Edgar B. C., Miles W. T., Green A. E. S., 1973, J. Geophys. Res., 78, 6595

* Edgar B., Porter H., Green A., 1975, Planet. Space Sci., 23, 787

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

* Eriksson A. I. et al., 2007, Space Sci. Rev., 128, 729

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

* Faure A., Gorfinkiel J., Tennyson J., 2004, J. Phys. B: At. Mol. Opt. Phys., 37, 801

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

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

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

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

* Gilet N., Henri P., Wattieaux G., Cilibrasi M., Béghin C., 2017, Radio Sci., 52, 1432

* Glassmeier K.-H., Boehnhardt H., Koschny D., Kührt E., Richter I., 2007, Space Sci. Rev., 128, 1

* Goetz C., Volwerk M., Richter I., Glassmeier K.-H., 2017, MNRAS, 469, S268

* Gulkis S. et al., 2015, Science, 347, aaa0709

* Hajra R. et al., 2020, A&A, 635, A51

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

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

* Hong J., Lee E., Min K., Parks G. K., 2012, Phys. Plasmas, 19, 092111

* Huebner W., Mukherjee J., 2015, Planet. Space Sci., 106, 11

* Itikawa Y., Mason N., 2005, J. Phys. Chem. Ref. Data, 34, 1

* Jun G., Quan-Ming L., 2007, Chin. Phys. Lett., 24, 3199

* Läuter M., Kramer T., Rubin M., Altwegg K., 2018, MNRAS, 483, 852

* Lavorenti F., Henri P., Califano F., Aizawa S., André N., 2021, A&A, 652, A20

* Lemaire J., Scherer M., 1974, Space Sci. Rev., 15, 591

* Lummerzheim D., Lilensten J., 1994, Ann. Geophys., 12, 1039

* Madanian H. et al., 2016, J. Geophys. Res. Space Phys., 121, 5815

* Mandt K. E. et al., 2016, MNRAS, 462, S9

* Meyer-Vernet N., Couturier P., Hoang S., Perche C., Steinberg J., Fainberg J., Meetre C., 1986, Science, 232, 370

* Myllys M. et al., 2019, A&A, 630, A42

* Pritchett P. L., 2010, J. Geophys. Res. Space Phys., 115, A10208

* Rego D., Prangé R., Gérard J.-C., 1994, J. Geophys. Res. Planets, 99, 17075

* Singhal R. P., Chakravarty S. C., Bhardwaj A., Prasad B., 1992, J. Geophys. Res. Planets, 97, 18245

* Sishtla C. P., Divin A., Deca J., Olshevsky V., Markidis S., 2019, Phys. Plasmas, 26, 102904

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

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

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

* Weaver H. et al., 1999, Icarus, 142, 482

* Zwickl R. D., Baker D. N., Bame S. J., Feldman W. C., Fuselier S. A., Huebner W. F., McComas D. J., Young D. T., 1986, Geophys. Res. Lett., 13, 401

## APPENDIX A: VALIDATION WITH PARTICLE DRIFTS

We validate the calculation of the particle trajectories, in the collisionless case, with a selection of standard particle drifts which are known analytically. These are the $E \times B$, curvature and grad-B drifts.
<figure><img src="image_12.png" /><figcaption>**Figure A1.** (a) Trajectory of an electron undergoing $E \times B$ drift in the xz-plane. (b) Percentage difference of the particle position in the $E \times B$ direction (x here) from the expected value.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure A2.** (a) Trajectory of an electron undergoing curvature drift in a cylindrical field. The colour bar corresponds to the time from the creation of the test particle. (b) Percentage difference of the particle position in the direction of the curvature drift (z here) from the theoretical value.</figcaption></figure>

### A1 E × B drift

The $E \times B$ drift test case is set up with electric and magnetic fields that are uniform and orthogonal. We have used fields of $B = 6$ nT $\hat{y}$ and $E = -2.4 \times 10^{-3}$ V m$^{-1}$ $\hat{z}$, which are typical of the quiet SW. These conditions result in a theoretical uniform drift velocity of

$$v_{E \times B} = \frac{E \times B}{B^2} = 400 \text{ km s}^{-1} \hat{x}. \qquad (A1)$$

The test particles were created with zero velocity in the y direction. This was required to confine the particles within the simulation domain for enough time to observe a significant particle drift. In Fig. A1, we present a particle with initial energy of 25 eV, but a range of particle energies corresponding to the energies of interest in this study have also been validated with a similar outcome. Fig. A1a shows the trajectory of the electron as it drifts along $+x$, with its motion confined to the xz-plane. The gyroaveraged velocity of the electron (see Fig. A1c) closely agrees with the expected drift velocity with variation of the order of 2 m s$^{-1}$ (1 in $2 \times 10^6$). We determine the expected position of the particle along the x-axis using $x_{\text{exp}}(t) = x_0 + v_{E \times B}t$. The deviation of the particle position from this

expected value is less than $4 \times 10^{-6}$ per cent across the 2000 km domain. The $E \times B$ drift is well replicated in the test-particle model.

### A2 Curvature drift

The curvature drift is examined by setting the electric field to zero and imposing a cylindrical magnetic field of $B = 6$ nT $\hat{\phi}$. We do not include a gradient in the magnetic field to validate the numerical treatment of the curvature drift alone, but as a result the field is not physical. The curvature drift velocity is given by

$$v_{\text{curv}} = \frac{v_{\parallel}^2}{\Omega R_{\text{curv}}^2} (R_{\text{curv}} \times \hat{b}), \qquad (A2)$$

where $\Omega$ is the electron gyrofrequency and $R_{\text{curv}}$ is the radius of curvature of the magnetic field line. In this cylindrical magnetic field, the radius of curvature is simply $\rho$, the polar radial vector.

Fig. A2a shows the trajectory of a test particle initiated with energy of 4.5 eV of which 4.02 eV is parallel to the magnetic field. The particle gyrofrequency is 1054 rad s$^{-1}$ and the average radius of curvature is 201.6 km. These result in an expected curvature drift velocity of $-6676$ m s$^{-1}$ $\hat{z}$. The gyroaveraged test-particle velocity
<figure><img src="image_14.png" /><figcaption>Figure A3. (a) Trajectory of an electron undergoing grad-B drift. (b) Percentage difference of particle position from the expected value in the direction of the gradient drift (x here).</figcaption></figure>

(blue, Fig. A2c) is centred on the expected drift velocity and varies by up to 25 m s⁻¹ from the average. This variation is small compared to the particle velocity which is of the order of 10⁶ m s⁻¹. The expected particle position along the z-axis is given by $z_{exp}(t) = z_0 + v_{curv}Bt$. The difference between the gyroaveraged particle position and the expected z position does not exceed 0.02 per cent over 20 s of particle drift, which is similar to the lifetime of the electrons in the cometary environment. Therefore, the test-particle model has been shown to accurately reproduce the behaviour expected in a curved magnetic field.

## A3 Grad B drift

The drift resulting from a gradient in the magnetic field is given by

$$v_{\nabla B} = \frac{mv_{\perp}^2}{2qB} \frac{B \times \nabla B}{B^2}. \qquad (A3)$$

For validation of the drift, the magnetic field is along the z direction with the field gradient in the y direction: $B = B(y)\hat{z}$. The background field strength is $B_0 = 4.85$ nT and the field gradient is $5 \times 10^{-3}$ nT km⁻¹.

Fig. A3 shows the behaviour of a test particle undergoing grad-B drift. The particle has an initial energy of 25 eV all of which is perpendicular to the magnetic field so the particle motion is confined to the xy plane (see Fig. A3a). This was required to ensure the test particle remained inside the simulation domain for long enough to observe a significant drift. The expected grad-B drift velocity of $5276$ m s⁻¹ is accurately reproduced in the model. The gyroaveraged particle velocity varies by up to $20$ m s⁻¹ from the theoretical drift but this has little impact on the position of the particle as the error in the drift direction tends to $5 \times 10^{-4}$ per cent (see Fig. A3b).

## APPENDIX B: ELASTIC SCATTERING

Electrons may undergo substantial scattering during elastic and weakly inelastic collisions. Although we separate collisional scattering of electrons from the energy-loss process (see Section 2.1.3), we account for the scattering from both rotationally elastic ($E_{Th} = 0$ eV) and vibrationally elastic scattering ($E_{Th} = 0.004$ eV).

Differential scattering cross-sections are given from 0.25 to 50 eV for scattering angles between 10 and 180°. Any scattering of less than 10° is neglected as this is approximately equivalent to no collision occurring. The differential cross-sections for electron impact on

water are from Cho et al. (2004) and Faure, Gorfinkiel & Tennyson (2004). Differential cross-sections are only given down to 20° in Cho et al. (2004) at 4, 6, and 8 eV, so we have extrapolated to 10° by multiplying the cross-section at 20° by 2.2. This is consistent with the ratio seen up to 30 eV.

Elastic collisions are treated in two steps. First, the total elastic collision cross-section, $\sigma^{el}(E)$, is used to determine whether a collision occurs (see equation 1):

$$\sigma^{el}(E) = 2\pi \sum_{X=vib,rot} \int_{10^{\circ}}^{180^{\circ}} \frac{d\sigma^{el,X}}{d\Omega}(E, \theta) \cdot \sin\theta d\theta. \qquad (B1)$$

This approach to compute $\sigma^{el}(E)$ reflects the probability of undergoing significant scattering, whereas the total vibrational cross-section in Itikawa & Mason (2005) includes collisions with little scattering (<10°). Below 0.03 eV and above 100 eV, the scattering cross-section is set equal to the elastic scattering ($E_{Th} = 0$) from Itikawa & Mason (2005), as the rotational transition is weak beyond these energies. Between these bounds and the energy ranges covered by Faure et al. (2004) and Cho et al. (2004), the cross-section is interpolated using a power law. After the cross-section is calculated, the scattering angle is determined from the shape of the differential cross-section using a random number. The azimuthal angle is randomly sampled from a uniform distribution. Between 0.25 and 50 eV, the shape of the differential cross-section is interpolated to the energy of the particle. For energies below 0.25 eV, we assume the scattering is isotropic; and above 50 eV, we assume the shape of the differential cross-section does not change.

## APPENDIX C: VALIDATION OF THE ENERGY DEGRADATION

We validate collisional algorithm used in the test-particle model by comparison to the energy loss calculated using the continuous slowing down approximation (CSDA). The theory of the CSDA and the derivation of the expected energy are outlined in Appendix D.

In this test case, we simplify the particle dynamics by setting the magnetic and electric fields to zero throughout the simulation domain. The particles are generated with a velocity purely in the x-direction and the neutral density is taken as uniform with $n_{H_2O} = 10^{10}$ cm⁻³.
<figure><img src="image_15.png" /><figcaption>**Figure C1.** (a) Comparison of electron energy variation in the test-particle model (Test PI X) to the CSDA (Theory X) for three test cases. (b) The energy difference between the test-particle model and CSDA for each of the test cases (see text for the conditions for each test case).</figcaption></figure>

We consider three test cases of increasing complexity to verify the energy-degradation implementation. The first includes only one energy-loss process with a constant cross-section of $\sigma = 5 \times 10^{-16} \text{ cm}^2$ and a threshold energy of $E_{\text{Th}} = 0.004 \text{ eV}$. In the second test case, there is still only one process but the cross-section varies with energy. We use the rotational transition $J = 0 \rightarrow 1$ with a threshold energy $E_{\text{Th}} = 0.004 \text{ eV}$. The final test case includes three collision processes with energy-dependent cross-sections. As well as the rotational transition, two vibrational excitations are included with threshold energies of $E_{\text{Th}} = 0.298$ and $0.458 \text{ eV}$. We exclude collision processes with higher threshold energies as the assumptions of the CSDA ($\Delta E \ll E$) would not be satisfied.

As the collisions are implemented as a stochastic process in the test-particle model, we have taken the average energy of 5000 particles at set distances over a 2000 km range. Each of the 5000 particles is created at the same energy and with the velocity along $+x$. These collisional tests have been carried out and validated at different initial particle energies; however, here we only show the test results for 25 eV electrons. The energy variation for each test case is plotted in Fig. C1a and the difference in energy between the test-particle model and CSDA is shown in Fig. C1b.

The first and second test cases show very close agreement between the CSDA theory and stochastic implementation over the whole energy range considered. The very low-threshold energy of the collision process (0.004 eV) in both of these tests means that the assumptions of the CSDA are well satisfied. The energy difference between the test-particle and CSDA models in both cases never exceeds 0.01 eV (see Fig. C1b) over a total energy loss of up to 8 eV.

In the third test case, the collision cross-section for the higher energy processes is small above 20 eV, so initially the energy loss is purely a result of the rotational collisions. When the vibrational processes become significant, the average macroparticle energy

deviates by up to 0.06 eV from the expected energy (light blue, Fig. C1b). The deviation results from the fact that the energy loss $\Delta E$ for the vibrational processes become more comparable to the energy $E$ of the particle (factor as low as 25), while the CSDA relies on $\Delta E \ll E$. This is still small compared both to the threshold energy of the vibrational collisions and the total energy loss of the particle ($\sim 10 \text{ eV}$). The test-particle model accurately captures the energy loss expected from the CSDA approach, while the CSDA is valid ($\Delta E \ll E$).

Generally, the CSDA would not be appropriate in the model as the energy loss in a single collision may be comparable to the electron energy (e.g. a 10 eV electron losing 7 eV by causing dissociation of a water molecule) and so a stochastic approach is preferred.

## APPENDIX D: CONTINUOUS SLOWING DOWN APPROXIMATION

The CSDA is used to model the energy degradation of a particle when the energy loss of a process, $E_{\text{loss}}$, is much smaller than the energy of the particle, $E$. In this regime, the loss of energy can be viewed as a continuous process rather than stochastic in nature. Such an approximation has been applied to proton aurora at Earth (Edgar, Miles & Green 1973; Edgar, Porter & Green 1975) as well as the deposition of energy by electrons and protons in Jupiter's atmosphere (Singhal et al. 1992; Rego, Prangé & Gérard 1994), considering energies in excess of 1 keV.

The rate of change of energy with distance is given by

$$ \frac{dE}{dx} = -n(x)L(E), \qquad (\text{D1}) $$

where $n(x)$ is the number density of the neutral gas the particle passes through along $x$ and $L(E)$ is the loss function. The loss function is
defined as

$$L(E) = \sum_{j} \sigma_{j}(E) E_{\text{loss},j}, \qquad (\text{D2})$$

where $E_{\text{loss},j}$ and $\sigma_j(E)$ are the energy loss per collision and the collision cross-section of the jth process, respectively. The particle energy at a position x is calculated by integrating the energy loss along the particle trajectory, evaluating the cross-sections at each point.

In this study, we apply the CSDA when the threshold energies of the processes are smaller than, but comparable to, the energy of the particle. This is at the limit of when the approximation is valid, but this is used only to validate the stochastic approach and so is acceptable for our purposes.

where the optical depth of the coma is small (Galand et al. 2016). The electron density is

$$n_e(r) = \frac{Q v_{\text{ioni}}^{\text{hv}} (r - r_{\text{com}})}{4\pi u_{\text{gas}}^2 r^2}. \qquad (\text{E1})$$

This has been replicated in the test-particle model by initializing all the photoelectrons with a uniform outflow velocity and setting the electromagnetic fields to zero. Fig. E1 shows the analytical and test-particle electron densities for a slice in the xy-plane containing the cometary nucleus ($z = 0$). The photoelectron density is within 15 per cent of the theoretical value throughout the coma. This indicates that the weighting and spatial distribution of the photoelectrons is implemented correctly.

# APPENDIX E: VALIDATION OF ELECTRON PRODUCTION SCHEMES

## E1 Validation of photoelectrons with radial outflow

We have verified the photoelectron weighting and production by considering the case of uniform radial outflow of electrons at velocity $u_{\text{gas}}$ from the nucleus. This can be solved analytically in the case

## E2 Validation of SW weighting

The weighting and initialization of the SW electrons is validated through comparison to pristine SW. We take parameters of the SW as $n_{\text{SW}} = 1 \text{ cm}^{-3}$, $V_{\text{SW}} = 400 \text{ km s}^{-1}$, and $T_{e,\text{SW}} = 10 \text{ eV}$. The electric and magnetic fields are uniform throughout the domain ($B = 6 \text{ nT}$, $\hat{y}$, $E = -2.4 \times 10^{-3} \text{ V m}^{-1}$, $\hat{z}$) and cause an $E \times B$ drift equal to the

<figure><img src="image_16.png" /><figcaption>Figure E1. Plasma density profiles for radial outflow of electrons from (a) the analytical solution and (b) the test-particle model. The lines of higher density in the test-particle simulation are generated by heavily weighted particles that originate in the inner coma.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure E2. Moments of the SW electrons for in pristine SW conditions. (a) Electron density in the $xz$ plane. (b) SW electron temperature calculated from equation (7) in the $xz$ plane.</figcaption></figure>
bulk SW velocity. For this test case, no neutral coma is considered
and therefore there are no electron-neutral collisions.

The electron density (see Fig. E2a) is close to 1 cm⁻³ throughout the domain, although it is enhanced by up to 3 per cent near the upstream boundary. This is a small error compared to the variation between cycles in the PiC simulations and will have little impact on any conclusions drawn from the model. The SW bulk velocity (not shown) and temperature (see Fig. E2) do not vary spatially and are

centred on the set parameters, 10 eV and 400 km s⁻¹, as required.
Therefore, the SW electrons are behaving as expected indicating that
the weighting and distribution of macroparticles across the domain
boundaries have been done correctly.

This paper has been typeset from a *TEX/LATEX* file prepared by the author.