© ESO, 2017

1. Introduction

From August 2014 to September 2016, the European Space Agency spacecraft Rosetta orbited comet 67P/Churyumov-Gerasimenko (67P, see Glassmeier et al. 2007), embarking a suite of instruments dedicated to the study of the comet’s plasma environment, the Rosetta plasma consortium (RPC). The Rosetta plasma consortium includes ion and electron spectrometers (RPC-ICA and RPC-IES, Ion Composition Analyser and Ion Electron Spectrometer), a magnetometer (RPC-MAG), and two instruments dedicated to the electron component of the plasma (RPC-LAP and RPC-MIP, LAngmuir Probes and Mutual Impedance Probe; Carr et al. 2007). Comet 67P has a complex plasma environment for which observations are a constant challenge. Most of the plasma boundaries and interaction regions seen with Rosetta result from an interplay between the neutral atmosphere, variations in the solar wind upstream flow, and particle ionisation processes, such as photoionisation, charge exchange, and electron ionisation, which lead to the creation of cometary ions. This study proposes to investigate this interplay in detail, from the point of view of a physical simulation, to aid in the interpretation of Rosetta’s datasets. The observational work performed with Rosetta on comet 67P’s neutral and plasma environment is first presented with a discussion of the characterisation of macroscopic boundaries, such as the bow shock and the cometopause, before the current modelling work is reviewed.

The neutral atmosphere was monitored in situ, for example, with the Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA, see Balsiger et al. 2007), which is composed of the comet pressure sensor (COPS) and of two mass spectrometers. The total neutral outgassing rate Q₀ from the nucleus was measured throughout the mission, showing asymmetries with respect to equinoxes and perihelion location; a few weeks after perihelion in mid-August 2015 (~ 1.3 astronomical units, or AU) a maximum of neutral outgassing activity was indeed observed (Hansen et al. 2016). The outgassing rates measured by ROSINA/COPS range from 10²⁵ to 10²⁸ molecules per second with a dependence for the inbound (outbound) legs of the comet’s orbit, where R_h is the heliocentric distance in AU. Estimates of the neutral outgassing rate, using RPC-ICA as a remote-sensing ion instrument of the charge exchange between solar wind ions and neutrals, are also in good agreement with ROSINA’s measurements (Simon Wedlund et al. 2016). Using neutral densities measured by ROSINA/COPS and a chemistry model, Vigren et al. (2016) compute electron densities in good agreement with in-situ RPC-LAP measurements, showing that ions and neutrals are strongly coupled.

Regarding the plasma environment during the early mission, Rosetta witnessed a transition from a solar wind-dominated region to a region dominated by the neutral coma coupled to ions and electrons (Yang et al. 2016). Using RPC-ICA, Nilsson et al. (2015a) and Nilsson et al. (2015b) monitored the birth and evolution of the comet’s induced magnetosphere. For the most part of the mission, RPC-ICA observe both solar wind ions at high energy, composed of H⁺, He²⁺, and often He⁺, and heavy cometary ions at low energy, themselves often composed of two populations, one cold below ~ 50 eV, one hot above ~ 100 eV (Nilsson et al. 2015a,b; Behar et al. 2016b; Simon Wedlund et al. 2016). Evidence of early ion pickup processes during the final approach towards the nucleus is also reported (Nilsson et al. 2015a; Goldstein et al. 2015; Coates et al. 2015). Using RPC-IES, Clark et al. (2015) show that the interaction of the comet with the solar wind is much more turbulent than anticipated, emphasising the role of mass loading in the formation of boundaries and distinct energy populations. Electron ionisation in the early coma is discussed by Galand et al. (2016) to explain the variations in RPC-LAP and RPC-MIP measurements over the winter hemisphere, which photoionisation alone can not account for. Broiles et al. (2016), modelling RPC-IES spectra with kappa distributions, find two cometary electron populations, one dense and warm, the other rarefied and hot, and both with temperatures of the order of 10⁵ K. They subsequently compare the evolution of each population at two different points in the mission lifetime, between the nearly inactive comet far from the Sun and the active one near perihelion. The authors conclude that these electron populations are most likely the rarefied hot solar wind halo electrons and some dense warm electrons of cometary origin.

In their multi-instrumental study, Edberg et al. (2015) find a highly structured low-energy plasma rapidly changing on scales of a few tens of kilometres during the orbit of Rosetta. Characterised by large density and electron temperature variations, this low-energy plasma was seen to evolve with respect to both the comet’s neutral atmosphere and the solar wind variations. These solar-wind variations may sometimes be driven by solar transients such as coronal mass ejections or co-rotating interacting regions (McKenna-Lawlor et al. 2016; Edberg et al. 2016a,b), the net effects of which are mainly to compress the comet’s plasma environment and increase the ionisation due to particle impact. Solar wind variations may also include those of the solar wind dynamic pressure. Using RPC instruments, Mandt et al. (2016) report the presence of plasma boundaries between April 2015 and February 2016. The authors in turn link these observations to the predominance of ion-neutral collisions and the so-called “collisionopause” deep in the comet’s atmosphere, modulated by the neutrals and the solar wind dynamic pressure.

Because the comet’s atmosphere typically extends in space on scales of hundreds of thousands of kilometres, mass loading of the solar wind is expected to play a major role in the formation of boundaries and plasma regions (Biermann et al. 1967; Szegö et al. 2000). Its effects, such as deflection and deceleration of the solar wind flow, are investigated with RPC-ICA by Behar et al. (2016b) and later by Behar et al. (2016a): they find that the observed solar wind deflection was the clearest indication of ongoing mass loading, with only little slowing down of the flow. Deflection of the solar wind from the Sun-comet line occurs because solar wind ions undergo a Lorentz force in an opposite direction to pickup ions. The pickup ions are accelerated mostly in the direction of the solar wind convection electric field E^sw, which is necessary to conserve energy and momentum between fast-moving solar wind ions and slow-moving cometary ions. As a result, a weak pile-up of the plasma may occur, which can be seen in the magnetic field data. This effect is reported at 67P by Koenders et al. (2016a), with strong signatures observed in RPC-MAG data at 2.0 AU that arise from the B-field piling up and draping. Using a hybrid model, these authors argue for a strongly deflected plasma flow in the first few tens of kilometres from the nucleus, in agreement with the particle instruments onboard Rosetta.

Ions that are newly picked up by the solar wind, originally at rest, will accelerate along E^SW to a maximum speed of 2 × U^SW in the situation where the Interplanetary Magnetic Field (IMF) is perpendicular to the flow of the solar wind, where U^SW is the solar wind speed; they will also gyrate about the magnetic field B^SW with a typical radius r_L = m_iU^SW/e | B^SW | , where m_i is the mass of the ion i and e is the elementary charge. A typical value of r_L at 1.3 AU is about 1.5 − 2 × 10⁴ km, which is of the order of magnitude of the macroscopic plasma boundaries found at weakly outgassing comets such as the bow shock and cometopause. Pickup ions adopt ring or ring-beam velocity distribution functions (Coates & Jones 2009), giving rise to instabilities which are in turn dissipated by the generation of ion cyclotron or mirror-mode waves (Volwerk et al. 2016). Though no ring-like distributions have so far been detected at 67P (Goldstein et al. 2015; Koenders et al. 2016b), Volwerk et al. (2016) report mirror-mode waves in the RPC-MAG data. Recently, ion acoustic waves were reported by Gunell et al. (2017) using RPC-LAP, RPC-ICA and RPC-MIP, calculating wave dispersion relations based on the ion distribution functions. Their excitation mechanism is still under debate, but such wave studies may contribute to our understanding of how phenomena on different scales, from the ion acoustic wavelength of ~ 10 m, via mean free paths for ion-neutral collisions of 1 − 10 km, to characteristic length scales for electric fields of a few hundred kilometres, interact to shape the plasma environment of a weakly active comet.

Due to the mass loading upstream of the nucleus, the super-magnetosonic solar wind may become sub-magnetosonic close to perihelion when the neutral outgassing is larger than about 10²⁷ s^-1 and a weak subcritical shock at low Mach numbers may form as a result, situated at a few thousand kilometres from the nucleus at the subsolar point (Coates & Jones 2009). This shock is distinct from other planetary bow shocks while most resembling that of Venus (Wallis 1973; Slavin et al. 1984; Russell 1985; Mendis et al. 1986; Coates & Jones 2009; Balogh & Treumann 2013) because the flow deceleration and energy dissipation may take place over an extended region of space. Since the first possible identification of a cometary bow shock at a comet (Russell et al. 1984), many observations of cometary bow shocks have been performed: for pre-Rosetta missions, the reader is referred to Coates (1995) and Coates et al. (1997). Cometary bow shock boundaries have also been modelled extensively (Biermann et al. 1967; Galeev & Khabibrakhmanov 1990) and more recently by Koenders et al. (2013, 2015). Despite Rosetta’s excursion on the dayside to nearly 1500 km, attempts to detect a bow shock forming upstream of the nucleus or indications of a change in cometary ion density were unsuccessful; this was in contrast to hybrid model expectations predicting such a boundary at about that cometocentric distance (Koenders et al. 2013). However, the absence of detection does not necessarily imply the absence of any bow shock, in particular as Rosetta appeared to be still inside a cometopause-like boundary (Edberg et al. 2016a). For comparison, Neubauer et al. (1993) report a standoff bow shock distance of ~ 1.7 × 10⁴ km at 1.0 AU at comet 26P/Grigg-Skjellerup, which has a similar outgassing rate at perihelion as 67P.

It is therefore more likely that Rosetta crossed another plasma boundary, such as the so-called cometopause, due to the spacecraft-nucleus distance of a few 100 km adopted throughout the mission. One definition of the cometopause is (i) the region where water-group cometary ions start to dominate the plasma over solar wind ions (Galeev et al. 1988). Gombosi (1987) and Cravens (1989) also defined the cometopause as (ii) the solar wind charge-exchange collisionopause. The nature of the cometopause has been vigorously analysed for comet 1P/Halley (see Gringauz et al. 1986; Gombosi 1987; Ip 1989; Reme et al. 1994; Tátrallyay et al. 1995; Coates & Jones 2009); at comet 26P, the inner boundary of the cometopause was found at about 1.8 × 10³ km (Johnstone et al. 1993). At 67P, Mandt et al. (2016) remarked that all of the Rosetta observations of collision-related boundaries occurred within the cometary ion-dominated cometopause shell, defined as in (ii), that is the collisionopause for charge exchange of solar wind protons. If defined as in (i), one may argue that the disappearance of the solar wind signature from the field of view of the RPC ion spectrometers in April 2015 (Nilsson et al. 2015b; Simon Wedlund et al. 2016) and its later reappearance in December 2015 (Nilsson et al., in prep.) may be the signature of the spacecraft gradually transiting in and out of a cometary ion-dominated region, and hence the cometopause. Although predicted by numerical simulations, a more refined analysis of the archived RPC data will likely decide this problem.

For a discussion of the plasma boundaries found so far at comet 67P by Rosetta and how they relate to earlier missions’ findings, the reader is referred to Mandt et al. (2016). For a discussion of more transient plasma boundaries as seen by the Rosetta mission, such as the diamagnetic cavity, the reader is referred to Goetz et al. (2016b), Goetz et al. (2016a) and Nemeth et al. (2016).

One key aspect involved in the formation of these transitions and boundaries is cometary ion production and how it relates to the solar wind plasma and the neutral atmosphere. Three main physico-chemistry processes are responsible for the creation of heavy cometary ions: photoionisation, charge exchange and electron ionisation, sometimes referred to as “electron impact ionisation” in older literature. They take the following form: Ambient neutral species X, mostly H₂O molecules at small velocities, become ionised in processes 1-3 (X → X⁺), after photons hν, solar wind ions Yⁿ⁺ of charge n +, or suprathermal electrons e⁻ collide with them. Cometary ions X⁺ are produced.

Loss processes of these cometary ions need also to be taken into account, namely, electron recombination with thermalised electrons: (4)This process results in a range of neutral atoms and molecules, through dissociative recombination.

Processes 1−3and their known impact on the cometary environment are described in more detail below. Because they change the overall ion composition and are efficient at different cometocentric distances, they constitute the backbone of any attempt to model the formation of the macroscopic plasma boundaries of a comet.

Photoionisation due to solar extreme ultraviolet (EUV) radiation, process 1, is one of the main contributors of water-group ions in the close vicinity of the nucleus because the atmosphere is nearly always optically thin for medium activity comets (Huebner 1985; Cravens 1987). In the process, the produced photoelectrons may be suprathermal, in turn ionising the local neutral species (Cravens 1991). A detailed 1D model of these aspects can be found in Vigren & Galand (2013) and Galand et al. (2016). Collision cooling may also occur, such that a population of cold electrons may arise from higher-energy photoelectrons.

Whereas photoionisation is efficient in the inner coma, the importance of charge exchange, process 2, in the upstream solar wind and near the cometopause was previously pointed out by Gringauz et al. (1986) and Gombosi (1987), motivated by Giotto and VEGA’s measurements at comet 1P/Halley. For instance, Gombosi (1987) used a two-fluid model including photoionisation and charge exchange to successfully simulate the observations made in and around the so-called cometopause boundary at 1P/Halley. At comet 67P, the continuous presence of high-energy He⁺ ions in the RPC-ICA spectra (Nilsson et al. 2015a; Simon Wedlund et al. 2016) is attributed to single charge transfer. Double charge transfer may also occur, as Burch et al. (2015) have shown early in the mission with the discovery of H⁻ negative ions in the RPC-IES electron spectra. Because charge-exchange reactions depend on the neutral atmosphere, it was shown by Simon Wedlund et al. (2016) that a remote-sensing estimate of the neutral outgassing rate could be retrieved from the in-situ ion measurements of He⁺/He²⁺ flux ratio. It is interesting to note that charge exchange of highly ionised solar wind ions (O⁶⁺, etc.) can lead to a bright soft X-ray and far-ultraviolet (FUV) emission in comets (Cravens 1997; Bodewits et al. 2004, 2012). Rosetta did not possess any X-ray measurement capability but embarked the ALICE UV spectrometer (Feldman et al. 2015). To the authors’ knowledge, no detection of charge-exchange induced FUV emissions have been detected with Rosetta so far.

Electron ionisation, process 3, occurs for electron energies larger than the ionisation potential at 12.65 eV in the case of a H₂O gas. It is one of the most important mechanisms in the inner coma of comets (Cravens et al. 1987; Gan & Cravens 1990; Galand et al. 2016). The origin of these suprathermal electrons is still debated: in the case of 67P, Clark et al. (2015) suggested that an admixture of photoelectrons and electrons that have been heated by pickup ion instability waves may best explain the RPC-IES observations of high-energy electrons in the 10 − 300 eV energy range.

As emphasised above, cometary plasma observations and their interpretations have been intimately connected to numerical or analytical models. Since the initial historical attempts (Alfvén 1957; Biermann et al. 1967), many modelling approaches have been conducted. Recently, these include purely hydrodynamic (Reyes-Ruiz et al. 2010), magnetohydrodynamic (MHD; Hansen et al. 2007; Rubin et al. 2014, 2015; Shou et al. 2015) or kinetic hybrid (Gortsas et al. 2010; Koenders et al. 2013, 2015, 2016a) approaches. Overviews of this vast subject, including observations, can be found for example in Szegö et al. (2000) and, more recently, in Gombosi (2015). The strength of a hybrid approach, with ions considered as kinetic particles and electrons behaving as a thermalised fluid, is to resolve the ion scales while keeping a reasonable computation time. It naturally includes kinetic effects such as thermal spreading, gyromotion effects, and asymmetries due to pick-up and J × B forces.

We present a new 3D cometary hybrid plasma model, capable of describing the stationary and dynamic environment of a comet. This model is firstly dedicated to the interpretation of ion measurements onboard Rosetta, such as those of RPC-ICA. It self-consistently takes into account photoionisation and charge-exchange processes, but also includes electron ionisation. We use this global model in stationary conditions to investigate quantitatively the comparative role of these ionisation processes in the formation of the dayside macroscopic plasma regions of comet 67P, such as the bow shock and the cometopause regions, at a periheliocentric distance of 1.3 AU. We first describe the model in detail, explain our methodology to interpret the simulations, and discuss the formation of the cometary ion environment, the bow shock and the cometopause with respect to each process separately and incrementally. Finally, we provide a discussion of our results in the context of Rosetta and the detection, or lack thereof, of such large plasma boundaries.

2. A new global 3D hybrid model

A 3D global hybrid model of solar wind-comet interactions is presented in this section. It is based on a hybrid-kinetic platform developed over the last fifteen years by our team to study the solar wind action on the plasma environment of various solar system bodies (Kallio & Janhunen 2002, 2003; Kallio et al. 2006a; Kallio & Jarvinen 2012; Jarvinen & Kallio 2014; Dyadechkin et al. 2015, and references therein). Models at Mars and Venus (Kallio et al. 2006b, 2012; Jarvinen et al. 2009) or Titan (Sillanpää et al. 2007) have given new insights regarding in-situ observations.

The present model is the adaptation of this 3D hybrid platform to a cometary environment. Recently, it was used to evaluate the mass-loading effects for the passage of comet C/2013 A1 “Siding Spring” near Mars (Gronoff et al. 2014) and to derive the large-scale electric and magnetic fields for the calculation of trajectories of ions and dust particles at comet 67P close to perihelion (Gunell et al. 2015). The architecture described below enables the model to be run until a stationary state is reached, as in the present study. Moreover, solar wind temporal perturbations and limited temporal changes in ionisation processes can also be investigated.

2.1. Quasi-neutral hybrid architecture

The model describes the cometary plasma using the quasi-neutral hybrid approach (QNH), in which positively charged ions are kinetic (macro) particles whereas electrons act as a charge-neutralising massless fluid moving at bulk velocity U_e. Using a cloud-in-cell implementation, the quasi-charge neutrality is ensured in a grid cell by: (5)where n is the density of electrons (subscript e) and ions (subscript i), e the elementary charge, and q_i the ion charge. Each macroparticle cloud is weighted and represents a large number of real particles, typically of the order of 10²³ ions, which all have the same characteristics such as density, velocity and temperature. On average, there are about 10⁷ macro-ions in the simulation box.

Ions and the electromagnetic fields E and B are self-consistently coupled via Maxwell’s non-radiative equations, Ampère’s and Faraday’s laws: (6)Starting with the magnetic field, the solver determines the total electric current density J and the electron charge density en_e in a cell from the quasi-neutrality assumption (5) and the definition of the electron current density: (7)where U_i is the ion bulk velocity.

The electric field is then determined by the electron fluid momentum equation (Ledvina et al. 2008; Kallio & Jarvinen 2012): (8)where η is the resistivity and p_e the electron thermal pressure at constant temperature T_e, where p_e = n_ek_BT_e, with k_B the Boltzmann constant.

The magnetic field can then be advanced in time by Faraday’s law, and once the updated E(t) and B(t) are known, the ions can be moved and accelerated by the updated Lorentz force: (9)The extra force that transfers momentum from the ions i to the neutrals s is the ion-neutral Langevin frictional term (Cravens 1987; Gombosi et al. 1996; Rubin et al. 2014). The constant k_i,s is the momentum transfer rate between ions, of density n_i and velocity U_i, and neutrals n_s,v_s. Because the neutral momentum is assumed to be much larger than the momentum transfer from ions, the momentum transferred to the neutrals is neglected.

Because of the QNH approach, finite ion gyro-motion effects and the so-called j × B Hall term are automatically included in the model, resulting in kinetic effects and plasma asymmetries. To focus on a specific area of the comet’s environment, grid refinements of different shapes, such as cubic, parabolic or spherical, may be used. The model runs at equilibrium in under t_f = 300 s and has a typical timestep Δt ~ 10^-3 s, small enough to converge towards a quasi-stationary state that satisfies the Courant-Friedrichs-Lewy condition needed for a stable numerical integration of the QNH equations (Birdsall 1991). For a more technical description, the reader is referred to Kallio & Jarvinen (2012) and, for the spherical version of the QNH model, to Dyadechkin et al. (2013).

Because comets are not known to possess an intrinsic magnetic field, the formation of tails and boundaries in the plasma is driven by the solar wind and the state of the comet’s neutral atmosphere. They interact through numerous physico-chemical processes that have to be consistently taken into account in any modelling attempt. Our 3D hybrid model includes photoionisation, charge exchange, electron ionisation, and electron recombination as individual probabilities. This implementation makes it possible, in one run, to follow each newly-born cometary ion particle cloud created by one specific process separately, so that the memory of the origin of each ion is kept throughout the simulation; applications of this technique are discussed in Sect. 4. Particle processes and their implementation are detailed in Appendix A.

2.2. Inputs

Table 1 summarises the input parameters of the model at the heliocentric distance R_h = 1.3 AU, corresponding to perihelion conditions for comet 67P. The hybrid simulation box was cubic and centred on the comet nucleus. The x,y,z dimensions of the box were ± 2.5 × 10⁴ km. Four spherical grid refinements were used, with a total of 92 640 cells in the simulation; the spatial resolution varied from 1600 km close to the external boundary of the simulation to 100 km within the first 800 km from the nucleus. At coordinates (0,0,0) a sphere of r_c = 2 km radius represented the comet’s nucleus, through which the magnetic field diffused freely due to the too coarse resolution around the nucleus. On average, 10 to 200 macroparticles per cell and per species (depending on the thermal velocity and density) were maintained in the interaction region using a split/join algorithm (Dyadechkin et al. 2013). The grid structure in the x − y plane is shown in Fig. 1. The Courant-Friedrichs-Lewy condition for a stable numerical propagation of the fields and particles was ensured everywhere in the simulation box: the signal speed was at worst Δx/ Δt ~ 13 500 km s^-1, with a constant small timestep Δt = 7.4 × 10^-3 s.

Fig. 1 Simulation box and grid refinement. The square simulation box is ± 25 × 103 km large and has an increasing resolution dx when coming closer to the nucleus. dx = 1600 km for x< 25 000 km, dx = 800 km for x< 12 800 km, dx = 400 km for x< 6400 km, dx = 200 km for x< 3200 km and dx = 100 km for x< 800 km. These limits apply to the centre of the cells.

2.2.1. Solar wind upstream parameters

The Sun was situated on the + x axis with the solar wind propagating from + x to − x at a speed of U^sw = 400 km s^-1. Solar wind usptream conditions were taken at 1 AU and then scaled to the perihelion distance of comet 67P. At 1 AU, the modulus of the unperturbed IMF intensity is typically of the order of 5 − 10 nT (Slavin & Holzer 1981; Luhmann et al. 1993). Assuming the higher end of this range, which seems to account better for the observed magnetic field intensities measured by RPC-MAG on Rosetta (Auster et al. 2015; Richter et al. 2015), a B_y component intensity of 7.69 nT was used for which a distance scaling in 1 /R_h was applied following the recommendation of Slavin & Holzer (1981). The Parker spiral angle was thus θ = 90° and the solar wind convection electric field was along the + z axis. At 67P’s perihelion, the IMF Parker spiral angle is about θ = 52°, which is expected to produce additional asymmetries in the B-field draping pattern (see Koenders et al. 2016a; Jarvinen et al. 2013, for Venus) and therefore we chose θ = 90° for simplicity. Ion and electron temperatures of T_i = 10⁵ K and T_e = 1.5 × 10⁵ K were chosen because they are typical values found in the solar wind (Slavin & Holzer 1981). In the hybrid model, by construction, the electron temperature remains constant everywhere in the grid.

2.2.2. Neutral atmosphere

The neutral atmosphere of the comet was taken into account by the neutral spherically-symmetric collisionless expansion model of Haser (1957), containing the photodestruction exponential term for parent molecules (as in Simon Wedlund et al. 2016). The comet outgasses neutral molecules s at a rate Q₀ (s^-1) and at speed v₀ (m s^-1), producing a neutral density n_s(r), (10)which depends on the cometocentric distance r. The total photodestruction rate of neutral molecules s is , with the photodissociation rate (superscript PD) and the photoionisation rate (superscript PI) of neutral species s creating ion i. All photo-rates were taken from Huebner & Mukherjee (2015) for low solar activity and scaled to the heliocentric distance.

Fig. 2 Velocity-dependent charge-exchange cross sections used in the hybrid model for single charge transfer of H+ and He2+ ions with H2O. The points represent the datasets (observations or model), the solid and dashed lines the polynomial fit of degree 4 made in each case. Coefficients for the fits are given in Table 2. The model results of Mada et al. (2007) are adjusted to match the laboratory results of Lindsay et al. (1997) at 1.5 keV impact energy (~ 979 km s-1). The fit is performed after the adjustment.

Daughter and grand-daughter species are ignored because they only start to contribute significantly to the atmospheric content at cometocentric distances of more than 30 × 10³ km. For illustration, at 1000 km cometocentric distance, the main parent molecule is H₂O and is expected to be about 90% of the total mass density of a comet’s atmosphere, whereas CO and H are of the order of 10% and 1% (Tenishev et al. 2008); this tendency remains true at larger cometocentric distances. Therefore, only H₂O is considered in the following. The neutral outgassing velocity in a comet such as 67P/C-G is radial from the comet’s surface and typically in the range 400 − 1000 m s^-1 as measured, for instance, by the ROSINA instrument (Hansen et al. 2016; Galand et al. 2016). It depends on the nucleus’ activity and heliocentric distance. At 67P/C-G’s perihelion at R_h = 1.3 AU, we chose the higher end of the range, v₀ = 1000 m s^-1. For the outgassing rate, we used the recommendation of Snodgrass et al. (2013) obtained from remote-sensing observations and scaled to the heliocentric distance, so that Q₀ = 4.89 × 10²⁷ s^-1. This is in reasonable agreement with the recent measurements of Hansen et al. (2016) made with instruments onboard Rosetta, such as ROSINA, and which have been corrected with an empirical model of the coma for spacecraft motion and position. Hansen et al. (2016) fitted a s^-1 power-law dependence for the inbound leg of the orbit resulting in an outgassing Q₀ = 6.56 × 10²⁷ s^-1 at 1.3 AU, which is reasonably close to the rate of Snodgrass et al. (2013) adopted hereafter. A higher Q₀ value combined with a lower neutral outflow velocity in the Haser model would have resulted in a similar overall neutral density in the simulation box.

2.2.3. Ion environment: production and drag

The equivalent photoionisation rate Q_i, in units of s^-1, in the simulation box is, in the case of an optically thin atmosphere, the integral over r of the neutral density multiplied by the photoionisation rate: (11)where R_box is the radius of the spherical simulation box, assuming it is symmetric with respect to the origin. Even with a cubic simulation box, as it was used in this study, this formula remains accurate because most of the atmospheric density is concentrated in the first few hundred kilometres around the comet. This is a small uncertainty compared with the uncertainties associated with cross sections and reaction rates, as will be discussed below. For the box size and outgassing rates chosen, the numerical application gives Q_i = 3.03 × 10²² s^-1.

Charge-exchange reactions for H⁺ impact on H₂O, noted , were calculated using the total theoretical cross sections of Mada et al. (2007) weighted by the molecular surface and increased by 25% to match the laboratory measurements of Lindsay et al. (1997) at 1.5 keV energy, equivalent to 980 km s^-1 impactor velocity. The calculated results of Mada et al. (2007) were preferred because they give the full energy dependence of the cross section, whereas Lindsay et al. (1997) measured only three energies. Mada et al. (2007)’s results, even increased by 25%, underestimate those of Lindsay et al. (1997) at 0.5 keV by about 15%; consequently, the cross section we adopted gives a lower limit for the charge-exchange efficiency. Experimental uncertainties are of the order of 10% (Lindsay et al. 1996). The corresponding cross sections for He²⁺ on H₂O, , were taken from the measurements of Greenwood et al. (2004), with uncertainties ranging from 25% at low energies to 10% at high energies. These uncertainties in turn result in propagated uncertainties in the hybrid model fields, and especially the density: calculated densities of cometary ions are estimated to be within 25% at most. Charge-exchange cross sections are presented in Fig. 2. Polynomial fitting parameters to both recommended charge-transfer cross sections and depending on the speed of the impactor (in m s^-1) are given in Table 2 for polynomials of order 4.

Electron ionisation and electron recombination use reaction rate constants from Cravens et al. (1987) for a Maxwellian disribution function, and Hollenbach et al. (2012), including all recombination channels for H₂O⁺, respectively, but calculated at a constant electron temperature of K. We note that this temperature is higher than the solar wind electron temperature of the hybrid model by a factor of 2: this value was preferred here in the calculation of reaction rates because it allows us to more accurately model the cometary bow shock region, as emphasised in Cravens et al. (1987). This is consistent with an increase by a factor of 2.4 in the EI reaction rate constant with respect to the nominal electron temperature in the hybrid model; correspondingly, recombination rates are multiplied by a factor of 0.7, thus decreasing the effect of electron recombination in the removal of H₂O⁺ ions. Table 1 gives the electron-related reaction rate constants adopted in the simulations.

For the ion-neutral drag rate, which is essentially the Langevin momentum transfer rate, a conservative value of k_i,s = 1.5 × 10^-15 m³ s^-1 was adopted for H₂O⁺-neutrals reactions following Cravens & Korosmezey (1986), Cravens (1987), Anicich (1993), Gombosi et al. (1996). The net effect of the drag term between cometary ions and neutrals is to slow down the newly born cometary ions to match the speed of the neutrals; this mostly takes place in the denser parts of the comet’s atmosphere. Cravens & Korosmezey (1986) and Körösmezey et al. (1987) also note that, for large relative velocities, momentum transfer is dominated by charge transfer reactions; this occurs for solar wind-neutral collisions upstream of the nucleus where charge transfer reactions are expected to play a major role. That said, the ion-neutral drag term does not result in any scattering of the motion of ions.

To minimise the emergence of large electric field gradients in a low-density medium, as is usually the case for comets, leading to numerical instabilities, a high resistivity of η ~ 34 kΩ m was chosen after several tests and kept constant everywhere in the simulation box. This technique is well-known in hybrid models (Kallio et al. 2006c, for Venus) and was first pointed out by Hewett (1980) to compensate numerically for very small densities in Eq. (8). The trade-off is that the small scales, typically below a few tens of km, may become unreliable. Due to the large mesh size with respect to that of the nucleus, the magnetic field was non-zero inside the nucleus and diffused into the very small obstacle.

2.2.4. Limitations

In their study of the bow shock in perihelion conditions, Koenders et al. (2013) introduced an extended upstream boundary condition on their reduced box size of x_max = 14 000 km to take into account the effect of solar wind mass loading at large distances upstream from the comet nucleus. We performed a test simulation (not shown) in the exact same solar wind upstream and cometary outgassing conditions: B_y = 4.9 nT, U^sw = 400 km s^-1, cm^-3 and Q₀ = 5 × 10²⁷ s^-1. We find that, even without taking into account this extended upstream boundary condition, our hybrid simulation fields, in this case magnetic field and solar wind density, are in good agreement with Koenders et al. (2013): a bow shock standoff distance is found at 2700 ± 100 km in our simulation, to be compared with their finding of ~ 2200 km. The reason for such a difference, without the extended boundary condition, may stem from the fact that in our approach to charge exchange and other processes, a probabilistic method of particle interactions is used that makes full use of energy-dependent cross sections, as described in Appendix A. Moreover, we are mostly interested in the relative effect of each process on the boundary positions and extent, so that the chosen upstream conditions are taken only as a reference. Using a multifluid MHD model, Rubin et al. (2015) find at 1.3 AU a cometary bow shock distance of 3000 km in similar conditions as those of Koenders et al. (2013), using the magnetosonic Mach number as a criterion to identify the boundary; these results are also in good agreement with ours.

The effect of a rotating anisotropic outgassing source of water on the overall plasma environment of the comet is outside the scope and spatial resolution of the present study. For the dwarf planet Ceres, modelling of rotating asymmetric water vapour release is currently in progress by Lindkvist et al., using a hybrid model (Lindkvist 2016). Preliminary results suggest that for a given water production rate, localised sources of water vapour affect the solar wind the most when located near the terminator. Differences in the size and shape of the interaction region may alter the distribution of the ionisation processes of H₂O through photo-, electron, and charge-exchange ionisation. The effect of asymmetric sources on the global plasma interaction is also predicted to be less important for low mass loading, where the interaction depth is much larger than the interaction region. Similar effects are expected at comet 67P. Recently, Huang et al. (2016) used a four-fluid MHD model to investigate the effect of a more realistic nucleus illumination model on comet 67P’s plasma environment at perihelion. They find a bow shock standoff distance increasing from 6000 km to 10 000 km when including an asymmetric gas outflow, suggesting that an equivalent higher outgassing rate in a symmetric Haser neutral model should be used to account for the increase in neutral density on the subsolar side. In both cases, rotation and anisotropic outgassing source, the coarse mesh resolution in our current cometary hybrid simulations prevents such a detailed analysis; the results of Huang et al. (2016) should in particular be kept in mind when interpreting the absolute bow shock standoff distances in our own simulations.

Another limitation, inherent to many hybrid models, including ours, is that we used a constant electron temperature in the box, equal to that of the solar wind electrons. Because electron temperatures are expected to drop rapidly below 10³ K in the inner coma (see Eberhardt & Krankowsky 1995, for 1P/Halley), one would expect an overestimate of the electron ionisation in this region in the model because ionisation frequencies exponentially decrease with T_e in this temperature range (Cravens et al. 1987). Similarly, because the electron thermal velocity is on average about 2500 km s^-1, the effect of electron recombination may be underestimated; however, because the mesh size we used close to the nucleus is 100 km, this effect is expected to be of relatively low importance. Consequently, uncertainties for the calculated densities in the hybrid model are low in the solar wind up to the bow shock-like region, but represent an upper limit inside the cometary bow shock.

3. Methodology

We performed several simulations in the conditions described in Table 1. In each simulation, we used photoionisation as a default process to create a background low-energy cometary ion density which helps to stabilise the simulation numerically and eventually decreases the noise. Simulation runs were performed successively as follows: photoionisation (PI run in the following), photoionisation and electron ionisation (PI+EI run), photoionisation and charge exchange (PI+CX run), and finally, photoionisation, electron ionisation and charge exchange (PI+CX+EI run, all processes at once). In the latter run, each cometary ion macroparticle cloud created by a specific process (photoionisation, electron ionisation or charge exchange) is followed individually, so that the fraction of each process in the creation of slow cometary ions can be evaluated.

In our study, the strategy adopted to compare our hybrid results more quantitatively process-by-process was to characterise the bow shock and cometopause boundaries in terms of standoff distance, curvature and flaring ratio for each simulation. Doing so, we could assess the respective effects of these physico-chemical processes on the formation of these large plasma boundaries on the dayside.

The nature of the cometary bow shock is discussed in Galeev (1987) and can be best approximated by a conic of revolution around the comet-Sun line, or, more precisely, a quasi-paraboloid of revolution close to the subsolar boundary (Farris & Russell 1994). Such a quadric has been successfully used to characterise Earth’s bow shock (Peredo et al. 1995). In the case of comets, a number of authors have remarked that, due to the large mass loading upstream of the nucleus, the cometary bow shock can be physically defined as the region where the magnetosonic Mach number M_ms becomes less than two instead of one (Schmidt & Wegmann 1982; Galeev & Khabibrakhmanov 1990; Coates et al. 1997); as such, the shock is weak and subcritical. Alternatively, we can also empirically define it as the region where the gradient of the solar wind velocity and that of the interplanetary magnetic field reach a maximum due to a large velocity decrease and heavy field compression. These two definitions have been concurrently explored in the following ways.

To quantify how the bow shock-like boundaries are affected by the physico-chemical processes, a quadric surface fit was performed either to the magnetic field intensity and the solar wind velocity boundaries, or to the 3D surface defined by M_ms = 2. For the former technique, we used an edge-detection algorithm based on the Canny method (Canny 1986); following Schug (2012), the boundaries were detected by first smoothing the 3D image and then by selecting the largest gradients of B and U^sw, which showed a steep transition from low to high | B | and from high to low U. A different threshold for each hybrid field was chosen to minimise the number of false positive edge detections. However, after several tests the external boundary could never be fitted with sufficient precision and fell within 30 − 50% of the position of the bow shock nose at best, because the surface of the boundary was never smooth enough to approximate it with a paraboloid function. With this in mind, it was found that the physical criterion of M_ms = 2 yielded consistently better paraboloidal surfaces to fit; the position of the fitted nose of the bow shock could always be found within one grid cell of the M_ms surface at the distance considered. In the following, the bow shock external surfaces are thus extracted from the hybrid model using the magnetosonic Mach number criterion. A discussion on the terminology of the so-called shock surface with respect to the Mach number is proposed in Sect. 4.1.

The general quadratic equation for a paraboloid of revolution k has a main axis along x with the tail of the paraboloid oriented antisunward, in Cartesian coordinates (x,y,z): (12)where A, B, F, G, H, and J are the unknown coefficients of the quadric to determine. Coefficients A, B and F have units of m^-1, J of m, whereas G and H are dimensionless. Using a least-squares numerical algorithm, the best paraboloid fit to the edge points is performed which gives a solution for each coefficient. The study of the residuals at the standoff point of the paraboloid ensures that the convergence and accuracy of the fit is adequate and within a maximum precision of 15%. For example, the bow shock distances, which are calculated either with the fit or directly inferred from the hybrid results, agree within 1% for the run with all processes included.

Using Eq. (12), we can define the standoff distance of the fitted paraboloid k as (13)and the mean curvature , taken in the special case close to the nose of the paraboloid (y = 0,z = 0), as (14)The mean curvature is connected to the mean curvature radius by: (15)The reader is referred to Goldman (2005) and Har’el (1995) for the derivation of these formulae.

Finally, the “flaring ratio” of the paraboloid, noted F_k, is defined for the bow shock as the ratio of the distances from the nucleus to i) the bow shock flank perpendicular to the Sun-comet line, and ii) the subsolar standoff distance (Mendis et al. 1986; Flammer & Mendis 1993, for comets). A measure of the overall bluntness of the paraboloid, this ratio is first evaluated in the x − y and x − z planes and then averaged between them.

This paraboloidal fitting approach is general, and standoff distance, curvature radius and flaring ratio can be applied to any paraboloid surface, such as the bow shock (noted with the underscript k = bs in the following, as in R_bs) and the subsolar-region cometopause (noted with the k = cp underscript, as in R_cp). The cometopause surface is defined in our simulations as the region where the solar wind and cometary ion densities become equal in accordance with Galeev et al. (1988); a paraboloid of revolution is a good approximation only around the nose region. An example of the M_ms = 2 surface for the bow shock, and the corresponding paraboloid fit performed on the run with all processes included overlaid to the magnetosonic Mach number in x − y and x − z planes, are shown in Fig. 3A.

Fig. 3 A) Magnetosonic Mach number Mms and B) plasma β in the x − y (top) and x − z (bottom) planes for the hybrid simulation, with all processes (photoionisation, electron ionisation and charge exchange) included. The 3D bow shock surface fitted for Mms = 2 in panel A is shown as a transparent meshed surface enveloping the inner plasma regions, the outer surface of which is delimited as a black line. The vertical dashed lines mark the subsolar (or standoff) position of the bow shock, Rbs = 6570 km (see text for explanation). The Sun is situated on the + x axis, with the solar wind propagating to − x at velocity . Because the undisturbed interplanetary magnetic field is directed along + y, the undisturbed convection electric field Esw = − Usw × Bsw is in the x − z plane. The position of the comet nucleus at coordinates (x,y,z) = (0,0,0) is shown by a cross.

4. Results and discussion

In the following, we investigate the results of our hybrid simulations and the relative effects of each ionisation process.

Fig. 4 A) Solar wind density; B) bulk velocity; and C) magnetic field intensity in the x − y (top) and x − z (bottom) planes, with all processes included (photoionisation, electron ionisation, charge exchange). Same remaining legend as in Fig. 3.

4.1. General run with all processes included

Figure 4 presents the run including photoionisation, charge exchange and electron ionisation processes at equilibrium for the solar wind H⁺ density (left column), bulk velocity (middle) and magnetic field magnitude | B | (right) in two planes: x − y plane (top, containing the magnetic field vector) and x − z asymmetric plane (bottom, containing the undisturbed convection electric field). The solar wind density changes by a factor of 3–4 across the boundary, with a large increase of proton density just behind it and an extended depletion of protons in the tail. The solar wind speed varies from 400 km s^-1 far upstream to less than 100 km s^-1 close to the nucleus. As for the density, a sharp boundary of parabolic shape resembling a bow wave can be seen where the gradient of the velocity is maximum. The magnetic field interaction region is somewhat similar, with the same sharp boundary separating an outer region of low interplanetary magnetic field intensity of about 7 nT from an inner region of high magnetic fields reaching up to about 25 nT on the dayside. In the tail, a cone of strongly depleted magnetic fields forms at 20 × 10³ km distance from the nucleus and measures 6 − 7 × 10³ km across.

For all parameters, the interaction regions are asymmetric in the x − z plane (y = 0), whereas they are symmetric in the x − y plane (z = 0). Because the magnetic field is oriented along the + y axis, the plane containing the undisturbed convection electric field E^sw = − U_e × B and thus the pickup ion motion is indeed in the x − z plane, giving rise to asymmetries in the distribution, composition and dynamics of the hybrid fields. A different clock angle would result in rotating the results along the x direction. This asymmetry is clearly demonstrated in Fig. 4.

4.1.1. Terminology

Figure 3 shows two important parameters that may help characterise the large shock-like boundary seen upstream of the nucleus: the magnetosonic Mach number (panel A) and the plasma β (panel B). The former controls the size and strength of the shock, the latter the level of turbulence (Schwartz 1998; Balogh & Treumann 2013). The plasma β is defined as the ratio of plasma thermal pressure to magnetic pressure: (16)As seen in Fig. 3B, β varies from less than 1.0 in the solar wind to greater than 2.5 just after the shock-like structure.

Similarly, the magnetosonic Mach number M_ms perpendicular to the IMF is defined as the following ratio: (17)where U^pl is the plasma velocity at the point considered, is the Alfvén velocity with ρ the total mass density, and is the speed of sound in the plasma. The ratio of specific heats, γ, (equal to the polytropic index in the case of an adiabatic process) equals for a mono-atomic gas in thermal equilibrium with three degrees of freedom; this approximation is warranted in the solar wind (Galeev & Khabibrakhmanov 1990), but may be too high in cometary ion-dominated regions that are essentially composed of H₂O⁺. Following the theoretical predictions of Schmidt & Wegmann (1982), Galeev & Khabibrakhmanov (1990), using a 1D hydrodynamic model, and later Khabibrakhmanov & Summers (1997), with an analytical model, also noted that a Mach number value M_ms = 2, due to the deceleration of the solar wind flow, was responsible for the growth of magnetosonic instabilities eventually leading to the formation of a cometary bow shock. A value close to 2 was observationally confirmed by Coates et al. (1997) at comet 1P/Halley and 26P/Grigg-Skjellerup. Using everywhere in the box for simplicity, the magnetosonic Mach number in our simulation with all processes included varies from M_ms ~ 4 in the upstream solar wind to M_ms ≲ 2 downstream of the boundary as seen in Fig. 3A. This indicates a transition from a super-magnetosonic to sub-magnetosonic regimes while transiting across the localised bow wave-like structure. This may justify the designation of a “bow shock-like” plasma boundary in the model.

4.2. Cometary water ions vs solar wind ions

To keep in mind one of the goals of 3D models, for example helping to interpret Rosetta’s complex datasets, a comparison between solar wind ions and cometary ions in the simulations and how they are affected by each physical process is a useful endeavour. Several aspects are of interest: solar wind and cometary ion dynamics, ion production, and ion composition.

Fig. 5 Solar wind proton H+ streamlines in the analysed simulations. A) Photoionisation only. B) Photoionisation + electron ionisation. C) Photoionisation + charge exchange. D) All processes. The colour code shows the bulk velocity of solar wind protons Up ranging between 0 (blue) and 400 km s-1 (red).

Fig. 6 Cometary ions H2O+ streamlines in the analysed simulations. A) Photoionisation only. B) Photoionisation + electron ionisation. C) Photoionisation + charge exchange. D) All processes. The colour code shows the bulk velocity of the cometary ions UH2O+, ranging between 0 (blue) and 500 km s-1 (red). The streamline tracing was performed by choosing equally-spaced seed points on the diagonal line between points of coordinate (x,y,z) = ( − 20,0, − 20) × 103 km and (20,0,20) × 103 km.

Fig. 7 Ionisation frequencies [s-1] of cometary water ions H2O+ due to A) electron ionisation () and to B) charge exchange () in the analysed simulations. The production rates are divided by the neutral cometary density ns to compare each process with respect to the photoionisation frequency s-1, which is shown as contour lines in red with delimited regions having values greater or less than .

4.2.1. Dynamics

Figure 5 shows the solar wind H⁺ streamlines for runs containing the following sets of processes: photoionisation (panel A), photoionisation and electron ionisation (panel B), photoionisation and charge exchange (panel C) and all processes (panel D). The colour code on the streamlines represents the bulk velocity of solar wind protons, ranging from 50–450 km s^-1. Only the x − z plane is considered here: whereas cometary ions follow the direction of the convection electric field, the solar wind ions move towards − E^sw because in their moving reference frame, the Lorentz force points predominantly towards the − E_sw direction (see detailed discussion in Kallio & Jarvinen 2012). The net effect on the proton streamlines of adding each process incrementally is a large solar wind deceleration upstream of the comet and, as for the magnetic field line draping (see section below), an increase in size of the interaction region. The addition of extra sources of cometary ions and the momentum exchange due to pickup becomes expectedly more pronounced when including successively electron ionisation and charge exchange, with charge exchange additionally removing solar wind ions. Both charge exchange and electron ionisation are much more effective than photoionisation at creating ions farther away from the nucleus. Where the solar wind enters the denser parts of the comet’s atmosphere, the efficiency of each of these two processes will increase, until the solar wind flow is sufficiently slowed down and deflected away from the inner regions of the comet’s plasmasphere. For charge exchange, this means that below about 150 km s^-1 impacting velocity, the process becomes exponentially less and less efficient as the cross section decreases (see Fig. 2): it is only barely compensated by the corresponding exponential rise in neutral density. By that point, photoionisation has however become much more efficient than any other process in creating slow-moving cometary ions.

When only photoionisation is taken into account, the solar wind decelerates from 400 km s^-1 along the x axis to less than 100 km s^-1 at a cometocentric distance of about 500 km with a large solar wind deflection within this limit towards − z. The addition of electron ionisation increases the size of the interaction region to about 2000 km and decreases somewhat the deflection because of this increase. Charge exchange, as seen before, expands the interaction region the most. Moreover, one can witness a transition between a Mach-like cone to a more symmetric plasma boundary that is reminiscent of a bow shock; for photoionisation, the (− x, − z) sector is mostly homogeneous, whereas a sharp decrease of U_p in the (+ x, + z) sector can be seen. Adding the other processes make the overall distribution of U_p much more symmetric with respect to the z = 0 line, with the velocity decreasing along the flank of the structures by the same order of magnitude.

Comparatively, Fig. 6 shows the streamlines of cometary H₂O⁺ ions in the x − z plane for the runs containing the different sets of process (same legend as in Fig. 5); this time the colour code represents the ion bulk velocity from a few tens to about 500 km s^-1. As expected, new-born cometary ions are accelerated by the convection electric field and follow gyromotions of typical radius r_g = m_iv₀/q_iB ~ 1.5 × 10⁴ km in the solar wind. Upon arrival closer to the nucleus where magnetic fields increase (see Sect. 4.3), the gyroradius diminishes accordingly, partaking in the piling up of ions in front of the nucleus. Because of this effect, cometary ion fluxes, which are expressed as density multiplied by velocity, increase close to the nucleus and are transported in the cometary tail. When turning on each process separately, fast cometary pickup ions above 400 km s^-1, are progressively replaced by medium-velocity ions, typically 250 km s^-1 and below. This is represented by the green-yellow region taking more and more space progressing from Fig. 6A to Fig. 6D. The closer to the nucleus, the newborn cometary ions have less and less time to accelerate on their gyromotion path. Even though the bulk velocity decreases on average by a factor of about 2, the density in the positive nightside (− x, + z) sector increases by an average factor of 5 − 6 (not shown), resulting in higher cometary ion fluxes that expand farther in the subsolar regions and in the (− x, + z) sector. Lower ion bulk velocities, below 150 km s^-1, reach a maximum regional extent for the PI+CX+EI run, with effects expanding upstream in the + x sector close to the bow shock-like region. Finally, in the − z sector, H₂O⁺ ions appear to be deflected earlier and earlier in their path towards the − x direction as a result of the magnetic field interaction region growing in size.

4.2.2. Production rates

Ion production rates for each process can be defined as Zhang et al. (1993): (18)with the ionisation frequency or ionisation rates in s^-1 of process j, producing ion i from neutral s, with j ∈ { PI,CX,EI }. Photoionisation rates are given by Huebner & Mukherjee (2015) and are tabulated in Table 1. The charge-exchange frequency in s^-1 for solar wind ion species i is simply (19)and thus depends on the solar wind ion flux (see Simon Wedlund et al. 2016). In contrast, the electron ionisation frequency in s^-1 depends on the total electron density, which is theoretically a mixture of cold and warm populations of electrons, but not in the hybrid model where T_e is constant: (20)where is the electron ionisation volumic rate in m³ s^-1.

Fig. 8 Fraction of H2O+ cometary ions (density ratios) created by A) photoionisation, B) electron ionisation and C) charge exchange in the analysed simulations. Density ratios, expressed in %, are calculated as with j ∈ [PI, EI, CX], for photoionisation, electron ionisation and charge exchange, respectively. The position of the comet nucleus at coordinates (x,y,z) = (0,0,0) is shown by a black cross. The white dotted box is the area zoomed-in in Fig. 9.

Calculating the production rates or the ionisation frequencies of cometary ions from each process separately can be done with only one simulation including all processes, as explained in Sect. 3. In this way, Fig. 7 shows the effective ionisation frequency in s^-1 of H₂O⁺ ions, calculated in the model as for electron ionisation (j = EI, panel A) and for charge exchange (j = CX, panel B), so that the production rates are compensated for the neutral atmosphere variations. These frequencies are directly comparable to the photoionisation frequency s^-1, constant everywhere in the box, with the corresponding value shown in the figure. On the dayside of the comet’s plasma environment, charge exchange dominates over all other processes from large cometocentric distances down to about 5000 km from the nucleus, where electron ionisation starts to compete. Regarding electron ionisation productions, a very sharp increase in the efficiency of the production mechanism occurs at around the bow shock distance; at this point, electron densities peak due to slowing down of the solar wind and magnetic field increased pile-up. Because the charge-exchange mechanism depends on velocity-dependent cross sections, this effect is more gradual because it depends on the overall ion flux, as shown in Eq. (19). Deep into the comet’s dayside coma, under ~ 1000 km cometocentric distance, photoionisation becomes dominant with efficiencies larger than charge exchange and electron ionisation by one to two orders of magnitude (regions in blue in Fig. 7).

The total production rates in s^-1 in the simulation box for photoionisation, electron ionisation, and charge exchange are summarised in Table 3. All three processes are of similar global efficiency in creating H₂O⁺ ions, with charge exchange being the most important, followed by photoionisation and finally by electron ionisation. Because electrons are modelled as an isothermal fluid, electron ionisation production rates presented here may give an upper estimate of the overall efficiency of the process in the real environment of the comet. To put these values and ratios in a wider context it is interesting to note that these rates are comparable to oxygen ionisation rates calculated by a gas dynamic model at Mars and Venus (Zhang et al. 1993), which reported values in the range 0.4 − 3.5 × 10²⁴ s^-1 at Mars and 1.0 − 7.0 × 10²⁴ s^-1 at Venus, depending on the process.

Fig. 9 Zoom-in of the fraction of H2O+ cometary ions created by A) photoionisation, B) electron ionisation and C) charge exchange in the analysed simulations. Ratios are in %. The position of the comet nucleus at coordinates (x,y,z) = (0,0,0) is shown by a black cross.

4.2.3. Composition

The production rates drive the density of cometary water ions produced and that of the remaining solar wind, because the only chemical loss considered is electron recombination. Consequently, in this hybrid model, H₂O⁺ ion densities are expected to mostly follow the production rates, especially in regions where pickup mechanisms have not had time to be initiated to change the overall composition significantly. Figure 8, representing the density fraction of H₂O⁺ produced by each individual process, shows exactly this trend: H₂O⁺ from photoionisation are confined to the inner coma (panel A), while electron ionisation (panel B) and charge exchange (panel C) play the major role beyond 1000 km from the nucleus. However, the pickup ion mechanism greatly affects the spatial structures, especially at large distances from the comet and in the tail, giving rise to an asymmetry in the direction of the convection electric field E^sw. At r ≥ 7500 km, charge exchange contributes to about 50 − 70% of the overall H₂O⁺ density, with electron ionisation of the order of 10 − 15%. The gyromotion along E^sw is also seen as streaks of enhanced plasma, though the noise in this area is more pronounced. Photoionisation is responsible for more than 50% of the cometary ion densities close to the nucleus and in the tail, with the pickup ion motion shaping the overall tail structure by carrying photoionised material from the dense dayside inner coma to the relatively depleted nightside outer coma.

Figure 9 displays these density fractions as contour plots zoomed in to the inner coma, within 1000 km from the nucleus. Up to 500 km, charge exchange contributes to about 30% of the cometary ion densities. The contribution of electron ionisation is somewhat similar and at least as important as that of charge exchange, notwithstanding the isothermal electron population considered in the model. The first gyromotion of the cometary ions is also visible in the drift of the created ions towards the positive nightside sector (− x, + z). On average, most of the ions present at these scales come from photoionisation, with the inner coma at less than 200 km distance being produced at 80 − 90% by photoionisation. The maximum dayside cometary H₂O⁺ densities due to photoionisation, charge exchange, and electron ionisation are 152.6, 4.7, and 6.8 cm^-3 respectively and occur in the first few hundred km from the nucleus. Omnidirectional particle flux ratios for H₂O⁺ ions follow a similar overall behaviour as that of densities, with flux maxima of 5.4 × 10⁷, 1.2 × 10⁷, and 1.65 × 10⁷ cm^-2 s^-1, for photoionisation, charge exchange, and electron ionisation, respectively. However, local differences between flux and density as well as their relative process-by-process contribution subsist, because of local changes in the average velocity of the ions.

Fig. 10 Magnetic field draping in the x − y plane in the analysed simulations. A) Photoionisation only. B) Photoionisation+electron ionisation. C) Photoionisation+charge exchange. D) All processes. The colour code shows the intensity of the magnetic field ranging from low (blue) to high intensity (red), with the same scale. The seed points for the field line tracing are shown as a horizontal dotted line along the x axis.

4.3. Magnetic field line draping and compression

Figure 10 presents the draping pattern of the magnetic field lines in the x − y plane for photoionisation (panel A), photoionisation and electron ionisation (panel B), photoionisation and charge exchange (panel C), and finally for all processes combined (panel D).

As expected, the size of the interaction regions increases dramatically when including more ion production processes. Although the magnetic field intensity is somewhat comparable when including each process separately, the typical pile-up region expands in the upstream solar wind in the y − z plane and along + x from about 1000 km (PI run) to 2500 km (PI+EI run) and finally up to 4000 km (PI+CX run). This shows that charge exchange has the most important effect of all in shaping the far upstream boundaries. We will see later in Sect. 4.4 that the effect of each process is somewhat cumulative when it comes to the interaction region size. Regarding the B-field draping, the photoionisation run results are reminiscent of a Mach cone interaction region forming towards a clear boundary. The addition of the other processes creates a more fully-fledged cometary bow shock with less fluctuations and steeper B-field gradients. The size of this interaction region depends on the solar wind and IMF values, the cometary ion source or ion production rate, and how extended it is in the solar wind. This effect can be seen in Fig. 8, as discussed in the previous section, with the main origin of the ions changing from charge exchange and electron ionisation to photoionisation when moving from the solar wind-dominated region to the inner coma.

As the interaction region expands, the magnetic field gradients seen at the point of transition between super-magnetosonic and sub-magnetosonic solar wind become somewhat steeper and smaller such that the transition into the coma is more gradual, especially when considering the PI+CX and PI+CX+EI runs. This effect may be illustrated by studying the width of the transition region, which we defined as the full width at half maximum of a Gaussian distribution fitted to the gradient of | B | along the + x axis. We denote this width Δx_{grad B}. It varies between about 400 km for the PI run, 500 km for the PI+EI run, 1000 km for the PI+CX run, and 1200 km for the PI+CX+EI run. The gradient is larger for photoionisation, and for photoionisation and electron ionisation by a factor of 3 on average. To complement this picture, the maximum intensity of the magnetic field on the x axis changes between 27.7 nT (PI run), 29.2 nT (PI+EI run), 25.7 nT (PI+CX run), and 25.9 nT (PI+CX+EI run), and is located at about 500 km cometocentric distance for photoionisation and ~ 800 km for all other runs. All of these values are within the first grid refinement of the simulation box, with a mesh resolution of 100 km, and therefore numerical noise should be minimal. The maximum intensity can be compared to that of the undisturbed IMF of | B^SW | ~ 7.7 nT; the pickup process and subsequent piling-up is thus responsible for a factor of 3.3 − 3.8 increase of B-field magnitudes on average, with photoionisation and electron ionisation being responsible for the largest increase overall. Though the size of the interaction region and the characteristic width of the transition from super-magnetosonic to sub-magnetosonic appear to be nearly cumulative when adding all processes, this is certainly not true for the maximum B-field intensity. All values pertaining to the magnetic field are given in Table 5.

Fig. 11 Bow shock 3D paraboloid quadric fits projected on the x − y (top) and x − z (bottom) planes for photoionisation (PI run), photoionisation + electron ionisation (PI+EI run), photoionisation + charge exchange (PI+CX run), and for all processes (PI+CX+EI run). The position of the nucleus is marked as a black cross. The bow shock surface is asymmetric with respect to the z = 0 line in the x − z plane.

4.4. Bowshock

As a summary of the discussions in the preceding paragraphs, we can quantify in more detail the evolution of the geometry of the bow shock boundary. Figure 11 presents the results of the fits to the bow shock in planes x − y (top) and x − z (bottom) for each run corresponding to the ionisation processes (photoionisation, photoionisation+electron ionisation, photoionisation+charge exchange and all processes combined). These fits are performed following the methodology described in Sect. 3: the shock-like boundary surface of each run is extracted for M_ms = 2 and then fitted in 3D with a paraboloid of revolution. Parameters of the fits are given in Table 4, and standoff distance and mean curvature radius, defined in Eqs. (14) and (15) from the fits, are given in Table 5. The average flaring ratio of the bow shock is also given; it has often been used to extrapolate satellite flyby measurements to the position of the bow shock’s subsolar point (Mendis et al. 1986; Flammer & Mendis 1993).

In Fig. 11, the shock surface appears noticably asymmetric in the x − z plane as expected from the solar wind flow direction around the obstacle. Two global effects are seen: one that increases the lateral size of the interaction region, represented by the change in curvature radius, the other that pushes the subsolar boundary towards + x, represented by the change in standoff distance.

As discussed in the previous sections, the effect of adding processes, from photoionisation to charge exchange, significantly increases the size and shape of the interaction region, which is better quantified by the value of the curvature radius around the nose of the bow shock. Whereas the shock with photoionisation only has a very narrow parabolic opening with a curvature radius km, the addition of electron ionisation and charge exchange shows shocks with much larger mean curvature radii, about 10 000 km and 18 000 km, respectively. The PI+CX+EI run has a curvature radius of km. Similarly, the incremental addition of all processes significantly increases the bow shock standoff distance, from about 830 km (PI run) to 2320 km (PI+EI run), and from 3840 km (PI+CX run) to 6570 km (PI+CX+EI run).

The effect of all processes together further amplifies both parameters, in a nearly incremental way, but not quite linearly. This is to be expected because the ions produced far away upstream, mainly by charge exchange and electron ionisation, will most influence the final shape and distance of the bow shock.

Our simulations show that the enhancement of the bow shock standoff distance and the curvature radius is maximum when charge exchange with solar wind H⁺ and He²⁺ ions is turned on. Second only to charge exchange, electron ionisation also plays a very important role, but has the most uncertainty in the model. In contrast, photoionisation is controlling the inner coma plasma, where cometary ions dominate (see next section).

The combination of these two effects, standoff distance and curvature, is further summarised by the average flaring ratio of the paraboloids . A value of 2.6 − 3.1 is derived from our simulations, depending on the processes taken into account. This value may be compared to the flaring ratio of 2.15 calculated by MHD models in the case of comet 1P/Halley, and subsequently applied to the Giotto measurements of comets 1P and 26P/Grigg-Skjellerup (Flammer & Mendis 1993). We note that our simulation’s interaction region is somewhat blunter, such that the bow shock flares out on the flanks more than in the historical studies. This is possibly due to the intrinsic characteristics of the coma of comet 67P, which is less active comet than 1P/Halley, and the use of a different model, namely MHD vs hybrid.

4.5. Formation of cometopause

We hereafter define the cometopause as the transition between a solar wind ion-dominated region and a region dominated by ions of cometary origin, assumed here to be H₂O⁺ (Gombosi 1987). Formally, (21)Figure 12 presents the ratio of cometary ions to solar wind densities in the asymmetric x − z plane for photoionisation (panel A), photoionisation and electron ionisation (panel B), photoionisation and charge exchange (panel C), and all processes (panel D). In this figure, the cometopause is seen as a white transition region, where the ratio R= 1, which is narrow at the subsolar point and more extended in the tail along the x axis. We performed paraboloid fits on the subsolar-region cometopause surface where the density ratio R= 1 to characterise this boundary with more precision around the subsolar point. Parameters of these fits are presented in Table 4. Because the cometopause surface is only of parabolic-like shape close to the nose of the cometopause, these fits are only valid for certain values of x so that x>x_lim. Analysis of the fits in terms of standoff distance, curvature radius and flaring ratio is presented in Table 5. Uncertainties in the determination of the standoff distance are higher than for the bow shock calculations. This is because the cometopause’s shape is more difficult to approximate with a paraboloid of revolution due to ion dynamics close to the nucleus, especially in the case of the PI run. The parabolas are almost identical between the PI+CX and PI+EI runs, showing that charge exchange and electron ionisation play a similar role. Although nearly constant for all the other runs, the curvature radius increases to about 900 km for the run with all processes included. The flaring ratio is below 2 in all cases, indicating a somewhat narrow and pointy cometopause envelope, with a more pronounced tendency for the PI+EI and PI+CX runs.

The position of the subsolar point, or standoff distance, from the nucleus varies approximately from 210 km (PI run) to 550 km (PI+EI and PI+CX runs) and 650 km (PI+CX+EI run). As explained in Gombosi (1987), the cometopause forms as a consequence of the rapid deceleration of the solar wind flow in this region; we can see this effect in our simulation in Fig. 4B for the PI+CX+EI run, or comparatively for all the runs in Fig. 5.

Asymmetry in the ion content is also seen in Fig. 12. The photoionisation results display a plume of cometary ions dominating towards the − z direction, an effect marginally seen in the other processes within the first few hundred kilometres from the nucleus on the nightside. Farther away from the nucleus on the nightside, the cometary ion tail globally moves towards the + z direction and is most extended for the PI+CX and PI+CX+EI runs. This is due more precisely to the cometary ion motion as explained in Sect. 4.2 (see also Fig. 6), because ions created very close to the nucleus, which is most relevant for the PI run, do not have the time to be accelerated by the convection electric field and will thus follow the flow of the solar wind towards − z. This is reminiscent of the hybrid simulation results of Koenders et al. (2016a) for a comet of intermediate activity. We also note that upstream of the cometopause the proportion of solar wind ions increases locally, seen as a redder region on Fig. 12, corresponding to the respective bow shock position of each run. Upstream and downstream of the bow shock, cometary ions densities increase with respect to solar wind ion densities. This is seen most clearly in Figs. 12B and C, and is linked to the increased efficiency of electron ionisation and charge exchange at these cometocentric distances.

As mentioned in the introduction, the original detection of a cometopause at comet 1P/Halley has been the subject of much debate (see Coates & Jones 2009, and references therein). The transition from solar wind-dominated to cometary ion-dominated regimes in our simulations is sharp, as predicted by other models (Gombosi 1987). It occurs in a region where the magnetosonic Mach number becomes much less than one, where the solar wind flow nearly stagnates and reaches velocities of 50 km s^-1 and below.

Fig. 12 Ratio R of cometary ion density to solar wind density in the x − z plane, where the undisturbed convection electric field lies. A) Photoionisation only. B) Photoionisation + electron ionisation. C) Photoionisation + charge exchange. D) All processes. The cometopause location is seen as the thin white region between the blue (cometary ion-dominated) and red (solar wind-dominated) regions, defined here as the positions where R= 1. The bow shock can be seen as a darker red paraboloid region that corresponds to an increased solar wind density upstream of the cometopause.

4.6. Discussion in the context of Rosetta

The absence of detection with Rosetta, or of positive indication throughout the mission, of a bow shock boundary at comet 67P has been the subject of much thought in the Rosetta community as part of an effort to bridge the gap between observations and modelling results. In late September 2015 at about 1.4 AU from the Sun, Rosetta was set to explore the dayside sector up to 1500 km in an attempt to probe the bow shock region (Edberg et al. 2016a). The dayside excursion lasted about two weeks with a heliocentric distance spanning 1.34 − 1.41 AU. Though no indication of a bow shock boundary was found, a coronal mass ejection (CME) on 5 − 6 October 2015 significantly compressed the comet’s plasma environment while Rosetta was at about 800 km from the nucleus and returning back to the inner coma. A sudden reappearance of a weak solar wind ion signal was observed in the field of view of the RPC-ICA ion spectrometer, whereas unprecedentedly high magnetic field amplitudes were observed, up to ~ 200 nT, in RPC-MAG. We note that the proton signal was short-lived, lasting a few minutes, and that this constitutes the only detectable signal of solar wind protons during that time. This is in contrast to previous model predictions, such as that of Koenders et al. (2013), which showed a shocked proton signal very close to the nucleus. This indicates that Rosetta was likely inside the cometopause boundary, not just inside a possible bow shock. For this observational period, the reader is referred to the case study of Edberg et al. (2016a).

Our simulations suggest that, with the solar wind upstream inputs we used in the nominal run, or even in the more quiet solar wind conditions of Koenders et al. (2013) that we used to validate our work, no bow shock could have been detected by Rosetta at about 1500 km cometocentric distance. The effect of charge-exchange reactions with both H⁺ and He²⁺ solar wind ions increases the interaction regions to several thousands of kilometres upstream of the comet. Also, the shock being weak, Rosetta’s transition into such a post-shock area may have been too gradual to detect any sharp changes in the plasma parameters. As a comparison, in the quiet solar wind case used in the present study, our simulation with all ionisation processes included predicts very little to no change in the magnetic field, 20 nT or so, for a virtual spacecraft moving out and in again in the simulation box corresponding to the dayside excursion of Rosetta through the comet’s environment. The cometary ion fluxes and ion densities in the simulation decrease by a factor of 5 or so while moving outwards, a result qualitatively seen in the Rosetta data (Edberg et al. 2016a). Because this virtual spacecraft’s orbit would effectively cross the cometopause boundary, the simulated solar wind fluxes increased by a factor of 30 or so; however, no solar wind signature was observed in the Rosetta ion spectrometer data (RPC-ICA or -IES). It is important to note here that Hansen et al. (2016) reported with ROSINA a neutral outgassing rate of Q = 1.58 × 10²⁹R^-7.15 s^-1 for the outbound leg of the orbit, which was larger around perihelion than that of the inbound leg. Thus, at 1.4 AU, Q ~ 1.4 × 10²⁸ s^-1, which is about a factor of 3 greater than the value adopted in our simulations. This would substantially increase the neutral atmosphere and plasma interaction region, likely expanding the cometopause in the direction of the subsolar point, with cometary ions becoming the major species at larger cometocentric distances and possibly at the maximum distance of Rosetta’s excursion. The effect of an illumination-driven model of the cometary neutral outgassing should also increase the bow shock distance substantially, as shown by Huang et al. (2016) with their multi-fluid MHD model.

Because the cometary plasma environment has been shown to be much more dynamic than expected for such a weakly outgassing comet, a correspondingly dynamic model should also be attempted. This is the purpose of a forthcoming study by our team, investigating the effect of solar-induced sudden variations, reminiscent of CMEs (solar wind) or solar flares (photoionisation), on the overall comet’s environment.

5. Conclusions

We present a new hybrid model of the cometary plasma environment, capable of simulating the macroscopic boundaries expected at comets and seen with pre-Rosetta space missions. The model is also capable of accounting for rapid temporal changes in the solar wind upstream conditions and activity of the comet, which is the subject of a companion paper in preparation.

In this study, we investigate in detail the stationary effect of photoionisation, charge exchange, and electron ionisation processes on the extent of dayside macroscopic boundaries such as the bow shock or the cometopause. We perform simulations in perihelion conditions (1.3 AU) for a virtual comet resembling comet 67P/Churyumov-Gerasimenko, the target of the Rosetta mission. A simulation containing all processes constitutes our base simulation from which all other simulations are run. Because we use a probabilistic algorithm for the creation of heavy cometary ions, these ions keep the memory of their original ionisation process, so that the effect of each ionisation process can be spatially monitored in the box.

The ionisation processes were also studied nearly separately. Charge exchange contributes to about 42% of the total production rate in the simulation box, with production rates from electron ionisation and photoionisation reaching 33% and 25%, respectively. The total H₂O⁺ ion production rates are calculated at ~ 9 × 10²⁵ s^-1 and were, interestingly, comparable to oxygen ionisation rates found at Mars and Venus (Zhang et al. 1993).

To quantify further the evolution of each ionisation process, we performed quadratic fits to the bow shock and dayside cometopause surfaces for hybrid simulations including, successively, photoionisation only, photoionisation and electron ionisation, photoionisation and charge exchange, and finally all processes combined. We explored three parameters of the fits, namely the standoff distance, the curvature radius, and the flaring ratio. Starting with photoionisation only and subsequently adding electron ionisation, and charge exchange processes, the bow shock region expands laterally in the y − z plane, while the standoff distance increases by a factor of 3 to 5 towards the subsolar direction. The run using solely photoionisation shows the least developed bow shock of all runs performed. Including the individual ionisation processes does not result in a linear combination of individual processes for these two parameters, especially the standoff distance, showing that a simple additive model of ionisation would be unable to approximate the more realistic case with all ionisation processes combined.

Charge-exchange reactions are found to be the dominant process triggering the formation of the fully-fledged bow shock boundary, followed by electron ionisation. Both charge exchange and electron ionisation are responsible for the bulk of the cometary ion densities far away from the nucleus for cometocentric distances exceeding 2000 km. Photoionisation, as expected, plays the major role around and inside the cometopause transition, below about 1000 km from the nucleus.

When Rosetta performed its month-long excursion to explore the dayside sector of the plasma environment of comet 67P (up to 1500 km from the nucleus) in September to October 2015, the RPC instruments detected no apparent changes in the magnetic field or the ion content on the first leg of the orbit outwards from the nucleus, which would be expected for a possible transition into the shocked region downstream of the bow shock-like boundary. This would imply that the shock structure, if it existed, must have been located at much larger cometocentric distances. We show in this study that this interpretation could be supported by our hybrid model results, with a bow shock position found at about 6500 km for neutral atmosphere and upstream solar wind conditions that are not dissimilar to the ones expected at the comet during the dayside excursion.

Acknowledgments

Rosetta is a European Space Agency (ESA) mission with contributions from its member states and the National Aeronautics and Space Administration (NASA). The work at University of Oslo was supported by the Research Council of Norway grant No. 240000. Work at Aalto University was supported partly by Academy of Finland, grant No. 251573. Work at NASA/SSAI was supported by NASA Astrobiology Institute grant NNX15AE05G and by the NASA HIDEE Program. Work at IRF was funded by the Swedish National Space Board under contracts 108/12, 112/13 and 96/15. Work at the Royal Belgian Institute for Space Aeronomy was supported by the Belgian Science Policy Office through the Solar-Terrestrial Centre of Excellence and by PRODEX/ROSETTA/ROSINA PEA 4000107705. Work at Umeå University was funded by the Swedish National Space Board (SNSB project 201/15). We wish to acknowledge the AMDA science analysis system provided by the Centre de Données de la Physique des Plasmas (CESR, Université Paul Sabatier, Toulouse) supported by CNRS and CNES. C.S.W. thanks Dr. Mea Simon Wedlund for helpful comments and suggestions.

References

Appendix A: Probability of interaction

Following Birdsall (1991), physico-chemical processes are implemented in the hybrid model as an individual probability of interaction in a manner similar to Monte Carlo particle-in-cell simulations (Vahedi & Surendra 1995).

• Charge exchange:X+Y⁺ → X⁺ + Y. It is taken into account by calculating the probability of interaction between a fast ion i of absolute velocity v_i with a neutral target s of density n_s, creating the new slow ion j: (A.1)using the energy-dependent charge-exchange cross section (in m²).

• Electron ionisation: X+e⁻ → X⁺+2e⁻. The probability of ionisation by impact of a suprathermal electron on a neutral species s is, creating in the process the slow ion species j: (A.2)where is the reaction rate constant for electron ionisation of neutral s (in m³ s^-1), depending on the electron temperature T_e.

• Electron recombination: X⁺+e⁻ → X. The associated probability of recombination of ion species j with an ambient electron impactor of density n_e is: (A.3)where is the electron recombination rate constant for ion species j (in m³ s^-1), depending on T_e.

For each species and interaction, a separate probability is evaluated so that each cloud of particle can be followed separately in the simulation. This provides the opportunity, in one run, to calculate how cometary ions created by a specific process affect the total cometary ion content. Cometary photo-ions due to photoionisation of the neutrals are produced analytically in the model from the 3D neutral profiles in Eq. (11), resulting in an ever-present background photo-ion population.

All Tables

All Figures

	Fig. 1 Simulation box and grid refinement. The square simulation box is ± 25 × 103 km large and has an increasing resolution dx when coming closer to the nucleus. dx = 1600 km for x< 25 000 km, dx = 800 km for x< 12 800 km, dx = 400 km for x< 6400 km, dx = 200 km for x< 3200 km and dx = 100 km for x< 800 km. These limits apply to the centre of the cells.
In the text

Fig. 2 Velocity-dependent charge-exchange cross sections used in the hybrid model for single charge transfer of H+ and He2+ ions with H2O. The points represent the datasets (observations or model), the solid and dashed lines the polynomial fit of degree 4 made in each case. Coefficients for the fits are given in Table 2. The model results of Mada et al. (2007) are adjusted to match the laboratory results of Lindsay et al. (1997) at 1.5 keV impact energy (~ 979 km s-1). The fit is performed after the adjustment.

In the text

Fig. 3 A) Magnetosonic Mach number Mms and B) plasma β in the x − y (top) and x − z (bottom) planes for the hybrid simulation, with all processes (photoionisation, electron ionisation and charge exchange) included. The 3D bow shock surface fitted for Mms = 2 in panel A is shown as a transparent meshed surface enveloping the inner plasma regions, the outer surface of which is delimited as a black line. The vertical dashed lines mark the subsolar (or standoff) position of the bow shock, Rbs = 6570 km (see text for explanation). The Sun is situated on the + x axis, with the solar wind propagating to − x at velocity . Because the undisturbed interplanetary magnetic field is directed along + y, the undisturbed convection electric field Esw = − Usw × Bsw is in the x − z plane. The position of the comet nucleus at coordinates (x,y,z) = (0,0,0) is shown by a cross.

In the text

	Fig. 4 A) Solar wind density; B) bulk velocity; and C) magnetic field intensity in the x − y (top) and x − z (bottom) planes, with all processes included (photoionisation, electron ionisation, charge exchange). Same remaining legend as in Fig. 3.
In the text

	Fig. 5 Solar wind proton H+ streamlines in the analysed simulations. A) Photoionisation only. B) Photoionisation + electron ionisation. C) Photoionisation + charge exchange. D) All processes. The colour code shows the bulk velocity of solar wind protons Up ranging between 0 (blue) and 400 km s-1 (red).
In the text

Fig. 6 Cometary ions H2O+ streamlines in the analysed simulations. A) Photoionisation only. B) Photoionisation + electron ionisation. C) Photoionisation + charge exchange. D) All processes. The colour code shows the bulk velocity of the cometary ions UH2O+, ranging between 0 (blue) and 500 km s-1 (red). The streamline tracing was performed by choosing equally-spaced seed points on the diagonal line between points of coordinate (x,y,z) = ( − 20,0, − 20) × 103 km and (20,0,20) × 103 km.

In the text

	Fig. 7 Ionisation frequencies [s-1] of cometary water ions H2O+ due to A) electron ionisation () and to B) charge exchange () in the analysed simulations. The production rates are divided by the neutral cometary density ns to compare each process with respect to the photoionisation frequency s-1, which is shown as contour lines in red with delimited regions having values greater or less than .
In the text

Fig. 8 Fraction of H2O+ cometary ions (density ratios) created by A) photoionisation, B) electron ionisation and C) charge exchange in the analysed simulations. Density ratios, expressed in %, are calculated as with j ∈ [PI, EI, CX], for photoionisation, electron ionisation and charge exchange, respectively. The position of the comet nucleus at coordinates (x,y,z) = (0,0,0) is shown by a black cross. The white dotted box is the area zoomed-in in Fig. 9.

In the text

	Fig. 9 Zoom-in of the fraction of H2O+ cometary ions created by A) photoionisation, B) electron ionisation and C) charge exchange in the analysed simulations. Ratios are in %. The position of the comet nucleus at coordinates (x,y,z) = (0,0,0) is shown by a black cross.
In the text

Fig. 10 Magnetic field draping in the x − y plane in the analysed simulations. A) Photoionisation only. B) Photoionisation+electron ionisation. C) Photoionisation+charge exchange. D) All processes. The colour code shows the intensity of the magnetic field ranging from low (blue) to high intensity (red), with the same scale. The seed points for the field line tracing are shown as a horizontal dotted line along the x axis.

In the text

Fig. 11 Bow shock 3D paraboloid quadric fits projected on the x − y (top) and x − z (bottom) planes for photoionisation (PI run), photoionisation + electron ionisation (PI+EI run), photoionisation + charge exchange (PI+CX run), and for all processes (PI+CX+EI run). The position of the nucleus is marked as a black cross. The bow shock surface is asymmetric with respect to the z = 0 line in the x − z plane.

In the text

Fig. 12 Ratio R of cometary ion density to solar wind density in the x − z plane, where the undisturbed convection electric field lies. A) Photoionisation only. B) Photoionisation + electron ionisation. C) Photoionisation + charge exchange. D) All processes. The cometopause location is seen as the thin white region between the blue (cometary ion-dominated) and red (solar wind-dominated) regions, defined here as the positions where R= 1. The bow shock can be seen as a darker red paraboloid region that corresponds to an increased solar wind density upstream of the cometopause.

In the text