Plasma source and loss at comet 67P during the Rosetta mission

Context. The Rosetta spacecraft provided us with a unique opportunity to study comet 67P / Churyumov-Gerasimenko from a close perspective and over a two-year time period. Comet 67P is a weakly active comet. It was therefore unexpected to ﬁnd an active and dynamic ionosphere where the cometary ions were largely dominant over the solar wind ions, even at large heliocentric distances. Aims. Our goal is to understand the di ﬀ erent drivers of the cometary ionosphere and assess their variability over time and over the di ﬀ erent conditions encountered by the comet during the Rosetta mission. Methods. We used a multi-instrument data-based ionospheric model to compute the total ion number density at the position of Rosetta. In-situ measurements from the Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA) and the Rosetta Plasma Consortium (RPC)–Ion and Electron Sensor (IES), together with the RPC–LAngmuir Probe instrument (LAP) were used to compute the local ion total number density. The results are compared to the electron densities measured by RPC–Mutual Impedance Probe (MIP) and RPC–LAP. Results. We were able to disentangle the physical processes responsible for the formation of the cometary ions throughout the two-year escort phase and we evaluated their respective magnitudes. The main processes are photo-ionization and electron-impact ionization. The latter is a signiﬁcant source of ionization at large heliocentric distance ( > 2 au) and was predominant during the last four months of the mission. The ionosphere was occasionally subject to singular solar events, temporarily increasing the ambient energetic electron population. Solar photons were the main ionizer near perihelion at 1.3 au from the Sun, during summer 2015.


Introduction
The Rosetta spacecraft first encountered comet 67P/Churyumov-Gerasimenko in August 2014. While other cometary missions to 21P/Giacobini-Zinner, 1P/Halley, 26P/Grigg-Skjellerup, and 19P/Borelly were fly-bys, Rosetta was the first mission to escort a comet and orbit around it to study the evolution of its coma and its interaction with the space environment (Glassmeier et al. 2007). As a comet gets close to the Sun, the volatiles in the nucleus sublimate, forming an extended coma, which is partially ionized by photo-ionization, electron-impact ionization, and charge-exchange with the solar wind (Cravens et al. 1987), creating the cometary ionosphere. The escort of comet 67P by Rosetta lasted for two years and ended on 2016 September 30. During this period, the comet evolved from a low outgassing object at 3.6 au from the Sun, where Rosetta encountered it, to an active comet near perihelion (1.24 AU) on 2015 August 13 and went back to its quiet state at 3.8 au from the Sun in September 2016. The data collected during this two-year escort period provided us with a unique opportunity to assess the evolution and variability over seasons, spatial locations, and heliocentric distances.
The Rosetta Plasma Consortium (RPC) instruments monitored the plasma environment during the full mission. Among them, the RPC-Mutual Impedance Probe (MIP) (Trotignon et al. 2007) and the RPC-LAngmuir Probe (LAP) (Eriksson et al. 2007) were used to derive electron densities and electron temperatures (Galand et al. 2016;Eriksson et al. 2017;Gilet et al. 2017;Heritier et al. 2017a). The bulk of the electron population comprises two distinct components. A warm component, with energies lying between 5 and 10 eV, is observed at all heliocentric distances by RPC-LAP . A cold component, with energies typically lower than 0.1 eV, was principally observed during high neutral activity Odelstad et al. 2018), at short heliocentric distances (<2 au). Such a cold temperature is most likely explained by electron collisions with neutrals, which are enhanced during high activity conditions. These main electron populations coexist with a less dense but more energetic component with energies typically between 10 and 200 eV (Clark et al. 2015). These energies are high enough to ionize the cometary neutrals. The ionizing electron fluxes are thought to be accelerated through wave-particle interactions (Broiles et al. 2016) or through the effect of the cometary ambipolar electric field (Madanian et al. 2016) accelerating solar wind electrons towards the comet (Deca et al. 2017). They were continuously measured by the RPC-Ion and Electron Sensor (IES) (Burch et al. 2007). When RPC-IES is corrected for the spacecraft potential bias (generally negative) measured by RPC-LAP (Odelstad et al. 2015;Odelstad et al. 2017), one can compute the local electron-impact ionization frequency (Galand et al. 2016;Heritier et al. 2017a). Apart from electron-impact ionization, photo-ionization through extreme ultraviolet (EUV) radiation is an active source of plasma around comet 67P, especially near perihelion where the solar flux is at its peak of intensity. There is no instrument on-board Rosetta to measure the EUV spectral solar flux. We therefore use the Thermosphere Ionosphere Mesophere Energetics and Dynamics (TIMED) -Solar EUV Experiment (SEE) (Woods et al. 2005) to measure the daily EUV spectral solar flux from the orbit of Earth and extrapolate it to the position of comet 67P, taking into account heliocentric distance and phase shift. Once the ionization frequencies are derived, it is essential to understand the neutral environment surrounding the cometary nucleus. The Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA) -COmet Pressure Sensor (COPS) (Balsiger et al. 2007) provided us with a continuous measurement of the total neutral number density whereas the ROSINA -Double Focusing Mass Spectrometer (DFMS, Balsiger et al. (2007)) supplied the neutral composition of the cometary coma. The main neutrals present are H 2 O, CO 2 , and CO (Fougere et al. 2016a;Le Roy et al. 2015), and these are the only chemical species considered here for our assessment of the ion total number density.
The electron density profiles of comet 67P were first assessed by Edberg et al. (2015) and Vigren et al. (2016) who used constant ionization frequencies that did not capture the high variability of the electron-impact ionization frequencies, which were later shown to vary by an order of magnitude within hours (Galand et al. 2016). The full multi-instrument approach, featuring time-dependent electron-impact ionization frequency, was first introduced by Galand et al. (2016) in order to study the ionospheric population at large heliocentric distances (3.1 au) and small cometocentric distances (10-20 km) in October 2014. The data-based ionospheric model used RPC-IES based electron-impact ionization frequencies, TIMED-SEE based photo-ionization frequencies, ROSINA-COPS total neutral densities, and neutral composition based on the work by Le Roy et al. (2015) using the ROSINA-DFMS instrument. The modelled total ion number density unambiguously matched with the electron number densities measured by RPC-MIP and RPC-LAP. Photo-ionization and electron-impact ionization were found to be the two main and sufficient sources of ions and their respective contributions were highly variable. Indeed, over the summer (northern, in the case of October 2014) hemisphere, where the neutral densities were the highest, photo-ionization was found to be the dominant source of ions whereas over the winter (southern) hemisphere, where the neutral densities were the lowest, the ionizing electron flux was more intense and both electron-impact ionization and photo-ionization were necessary to explain the plasma densities observed by RPC-MIP and RPC-LAP. This anticorrelation between neutral density and ionizing electron flux was not an isolated case and is illustrated in many of the cases studied in this paper. As this model is able to disentangle the sources contributing to the RPC-MIP and RPC-LAP plasma density observations, it was particularly useful to investigate unusual events, such as the impact of a cometary outburst on the ionosphere . The same data-based ionospheric model was also used to assess the generation of the near-surface cometary ionosphere probed during the end of mission, in September 2016 (Heritier et al. 2017a). At that time, electron-impact was the overwhelming source of ionization and the variations of the local ionizing electron fluxes strongly influenced the ionospheric densities along the spacecraft trajectory. While the ionospheric model of Galand et al. (2016) used constant neutral and ion velocities, the model was upgraded to be applied to the end of mission study (Heritier et al. 2017a) to feature neutral and ion velocity profiles (calibrated to ROSINA-COPS) reproducing the expansion and acceleration of the cometary gas as it moves away from the nucleus. This new feature was essential to explain the near-surface plasma densities observed by RPC-MIP and RPC-LAP as most of the acceleration of the outgassing flow takes place in the first few kilometres above the surface (Tenishev et al. 2008;Fougere 2014). This allowed us to constrain the cometary gas expansion velocity, which is a parameter that has not yet been directly reported by Rosetta instruments even though some estimations were made by the Microwave Instruments for the Rosetta Orbiter (MIRO, Marshall et al. (2017); Gulkis et al. (2015); Lee et al. (2015)). As the trajectories studied in the present paper do not vary in a similar way to the end of mission descent, we use a constant neutral velocity profile (while still allowing margins of error in these velocities) in our ionospheric model. Heritier et al. (2017a) also computed the charge-exchange frequencies between solar wind protons and cometary neutrals. As this process was found to not heavily influence the absolute ion number density at the location of the spacecraft, it is neglected in the present study.
In this paper we provide an overall view of the ionospheric drivers of comet 67P by covering a wide range of phases throughout the mission, in addition to the cases presented by Galand et al. (2016), Hajra et al. (2017), andHeritier et al. (2017a). We identify the ionospheric drivers throughout the two-year escort phase and their variability over time and heliocentric distance. As Galand et al. (2016) studied pre-perihelion cases, we focus in this paper on complementary post-perihelion case studies where we apply the same data-based ionospheric model in order to validate the findings of Galand et al. (2016) and Heritier et al. (2017a), and test the robustness of this model to singular events, such as the impacts of co-rotating interacting regions (CIRs) onto the cometary atmosphere or the ability of the model to reproduce measurements obtained during and after abrupt changes in the cometocentric distances covered by Rosetta. This paper mostly assesses the total ion number density without giving further details on the ion chemical composition. The reader is invited to consult Heritier et al. (2017b), Beth et al. (2016), Fuselier et al. (2016 for more information on the ion composition of comet 67P. Section 2 presents an overview of the external conditions encountered by the comet in terms of neutral environment, solar flux, and loss processes influencing the ion number density. In Section 3 we present the data-based ionospheric model and its results in five different post-perihelion cases. Finally, Section 4 presents a two-year comparison and interpretation of the different ionization sources in the vicinity of comet 67P, incorporating our earlier works (Heritier et al. 2017a,b;Galand et al. 2016).

Neutral conditions
The Rosetta spacecraft escorted 67P over an elliptic arc along its orbit, witnessing the comet approach and retreat from perihelion (August 2015) until the end of mission (September 2016), at 3.8 au from the Sun. ROSINA-COPS therefore measured an increase and a decrease in the outgassing rate by about a factor 1000 (Hansen et al. 2016). Figure 1 shows the time evolution of an approximated measurement-based global outgassing rate Q (black dots, bottom panel) as seen by the Rosetta spacecraft. It is an estimate assuming a Haser model (Haser 1957) and a constant outflow velocity v of 600 m.s −1 such that Q = 4πr 2 nv, where n stands for the total neutral number density measured by ROSINA-COPS. The red line shows the cometocentric distance, r, of the Rosetta spacecraft. The cometocentric distance is highly correlated with the outgassing rate as the Rosetta spacecraft had to moderate its exposure to dust in order to safely operate using its star tracker. Consequently, the measurements of the total neutral gas density by ROSINA-COPS were not directly correlated to the activity of the comet. The highest number densities measured (apart from the end of mission near-surface measurements) did not occur near perihelion but in May 2016, when the spacecraft was very close to the nucleus (∼7 km) although the outgassing rate was rather low (Q ∼ 10 26 s −1 ).
The peak in the local outgassing probed by ROSINA-COPS occurred slightly after perihelion, as pointed out by Hansen et al. (2016). This is due to latitude effects. At that time, spring (almost summer), was in the southern hemisphere and autumn in the northern hemisphere (see Figure 1 and Table 1). However, the Rosetta spacecraft was covering the northern hemisphere during perihelion (2015 August 13, Figure 1 top panel) and the whole first part of August until the 23, before visiting the southern hemisphere from August 23 to September 4. After September 4, Rosetta returned to the northern hemisphere and we can see the local outgassing rate decreasing accordingly in Figure 1, shortly before the excursion. More details on the latitude effects on the near-perihelion outgassing rates can be found in Gasc et al. (2017a).
In general, seasons introduce a significant difference in the local outgassing rate over the different hemispheres of the nucleus. The mission started with a northern hemisphere summer and a southern hemisphere winter where seasonal variations in the local outgassing were very intense (Hässig et al. 2015;Galand et al. 2016) until the equinox on 2015 May 10 where the seasonal effects reversed (autumn in the northern hemisphere and spring in the southern hemisphere). These effects gained in intensity after the solstice in September 2015. Interestingly, after the outbound equinox on 2016 March 21, even if the northern hemisphere went from winter to spring and the southern hemisphere went accordingly from summer to autumn, the local outgassing remained intense in the southern hemisphere (Gasc et al. 2017a). While the sub-solar latitude was in the northern hemisphere, the southern hemisphere outgassing was dominated by CO 2 that has a sublimation temperature lower than of H 2 O, of which the outgassing pattern followed the sub-solar latitude. A summary of these conditions can be found in Table 1.
An overview of the variation of the neutral composition over the mission can be found in Fougere et al. (2016a), Luspay-Kuti et al. (2015), and Gasc et al. (2017a). While the variations of the neutral densities have an important impact on the ionospheric densities, the neutral composition does not alter the total number of ions in a significant manner (though it does alter ion composition). Local CO 2 number densities can be occasionally higher than H 2 O. Galand et al. (2016) showed that, even though the CO 2 cross sections for photo-ionization and electron-impact ionization are overall higher than those of H 2 O (Itikawa 2002;Itikawa & Mason 2005), the effect of these variations on the total ion density is rather small. Indeed, the photo-ionization frequency is only increased by a factor 1.18 at most from a pure H 2 O to a half CO 2 half H 2 O mixture, while the electron-impact ionization frequency is increased by about 1.14 or less. [NH spring / SH autumn] Outgassing in the southern hemisphere is still higher due to the CO 2 in the southern hemisphere, which has a lower sublimation temperature than H 2 O. The two hemispheres have about the same magnitude of outgassing.

Solar flux and photo-absorption
The solar flux is a key driver in the ionospheric population. Extreme ultraviolet [EUV] radiation is able to efficiently ionize the main neutral species of the coma of 67P.

Solar flux
The number of sunspots has an eleven-year periodicity and the EUV solar flux intensity varies accordingly. Some smaller and less predictable timescale variations also occur. Figure 2 shows the time evolution of the number of sunspots as computed by the Sunspot Index and Long-term Solar Observations (SILSO World Data Center 2000-2017, gathering data from several solar observers (30% observatories and 70% amateur astronomers). The blue dots correspond to the monthly averaged sunspot number and the black line is a yearly moving average applied to these monthly sunspots to emphasise the 11-year period. The two-year escort phase of the Rosetta mission, shown in red, took place during a decreasing phase of the solar cycle. The photo-ionization frequencies at a fixed heliocentric distance are therefore expected to decrease during this period.
To assess the photo-ionization of the coma, we need a daily observation of the spectral solar flux. It is continuously mea-sured by the Thermosphere Ionosphere Mesosphere Energetics and Dynamics (TIMED)-Solar EUV Experiment (SEE) (Woods et al. 2005), from Earth's orbit, adjusted to 1 au from the Sun. Most of the recent literature uses the V11.02 of the Level 3 (L3) TIMED-SEE solar irradiance datasets and complementary solar models. A new version, V12.0, was released on 2017 October 31, featuring solar fluxes with overall higher intensity. The effect of the update is assessed in this section. Alternatively, the Flare Irradiance Spectral Model (FISM, Woods et al. (2006); Peterson et al. (2008)) also provides solar irradiance datasets in the EUV range. FISM uses observations from TIMED-SEE, the Solar Stellar Irradiance Comparison Experiment (SOLSTICE) on the Upper Atmospheric Research Satellite (UARS) and, the Solar Dynamics Observatory (SDO) EUV Variability Experiment (EVE). However, SOLSTICE-UARS only measures far ultraviolet (FUV) and ultraviolet (UV), from 120 nm to 420 nm, which are not energetic enough to contribute to photo-ionization of the main neutral species at the comet. Due to a power anomaly that occurred on 2014 May 26 (Milligan & Chamberlin 2016), SDO-EVE no longer gathers data within the 6 to 37 nm wavelengths. Therefore, most of the FISM data relevant to the photo-ionization of cometary neutral species originate from TIMED-SEE.
The spectral daily solar flux is first corrected for phase shift, as the part of the Sun facing the Earth is not necessarily the same that faces comet 67P. Therefore, a time shift, varying between 0 and 12 days depending on the day of the escort phase, was applied. Furthermore, as the solar flux is inversely proportional to the heliocentric distance, the values of the solar spectral flux at 1 au must be divided by the square of the heliocentric distance (in au) of comet 67P to obtain the local solar flux near the comet.
The solar spectral flux F is related to the photo-ionization frequency through the photo-ionization cross section of the neutral species considered. For instance, the photo-ionization frequency of H 2 O is obtained through where σ hν,ioni H 2 O (λ) is the photo-ionization cross section of H 2 O, taken from Vigren & Galand (2013), λ th corresponds to the ionization threshold wavelength below which absorption of solar photons possibly leads to ionization (about 100 nm for H 2 O), and λ min (here 0.1 nm) is a wavelength below which σ hν,ioni n (λ) and F(λ) drop drastically. Figure 3 shows the time evolution of the H 2 O photoionization frequency at Earth's orbit (dashed curve), at 1 au from the Sun, and at the location of comet 67P (solid curve). The ionization frequencies are computed using the three different datasets introduced earlier: TIMED-SEE L3 V11.02, its recently-updated version TIMED-SEE L3 V12.0, and FISM. The photo-ionization at 1 au is subject to a decrease over the solar cycle (see Figure 2) while the photo-ionization at the location of the comet is mostly sensitive to the heliocentric distance and peaks close to perihelion (2015 August 13), unlike the outgassing rate that peaks slightly after (see Figure 1). For a given heliocentric distance, the photo-ionization frequency at the comet is found to be higher pre-perihelion than post-perihelion. The effects of the time shift applied to extrapolate the solar flux from Earth's orbit to the position of comet 67P are visible in Figure 3 through the phase differences in the oscillations of the ionization frequencies. The phase shift is maximum (180 • , or 12 days) in February 2015, when comet 67P and the Earth are located at exactly opposite sides of the Sun, and minimum (0 • , or 0 days) in March 2016, when the two bodies are aligned on the same side of the Sun. There is a significant increase of about 30% for the photo-ionization frequencies computed with TIMED-SEE between V11.02 and V12.0. FISM lies in between these two versions in terms of magnitude. In Section 3, we chose to use TIMED-SEE V12.0 as it is the most recent released version of the EUV solar flux. Nevertheless, differences between the different solar irradiance datasets, in terms of ionospheric densities, are not significant at the large heliocentric distances cases treated in this paper. However, they should be important near perihelion.

Photo-absorption and Beer-Lambert law
The solar flux is attenuated by the cometary neutrals along the photons' trajectory. This phenomenon is well described by the Beer-Lambert law (Rees 1989). If we call s the coordinate along the trajectory of the solar flux and s 0 the coordinate of interest, the spectral solar flux at coordinate s 0 and for a given wavelength λ can be written as where n stands for any neutral species present along the photons' trajectory, σ hν,abs n (λ) is the cross section of photoabsorption at wavelength λ for species n, and n n is the number density of n. The quantity within the exponential function is referred as the optical depth τ(λ, s 0 ). With spherical symmetry, the equation can be solved analytically in the reference frame of the comet to get F(r, χ, λ), where r is the cometocentric distance and χ is the solar zenith angle at the spacecraft location (Beth et al. 2016). Fig. 4. Total ion density profiles as a function of cometocentric distances computed using an ionospheric model (Heritier et al. 2017b). The blue line corresponds to the output of the model neglecting both photo-absorption and electron-ion dissociative recombination (DR), the orange is considering only photo-absorption assuming a 90 • incidence angle, the yellow is corresponding only DR (assuming T e =200 K), and the purple is considering both effects.
We have assessed the effect of the solar absorption on the ionospheric densities as a function of the cometocentric distance using the ionospheric model presented in Heritier et al. (2017b). However, unlike this previous study in which adiabatic expansion was included, the runs here presented assume a constant neutral outflow velocity, as the impact of adiabatic expansion onto the total ionospheric density is not critical, except for the final plunge onto the comet. Figure 4 shows the ion total density profile as a function of the cometocentric distance considering (orange curve) and neglecting (blue curve) photo-absorption through the Beer-Lambert law. We have assumed an incidence angle, χ, of 90 • as the Rosetta spacecraft was mostly located in the terminator plane. The effects of electron-ion dissociative recombination (DR), which results in the formation of dissociated neutrals through the recombination of an ion and an electron (discussed in Section 2.3), are also displayed. The total ion density profile was computed for a typical day near perihelion (2015 August 31) where the spacecraft was located 380 km from the nucleus surface. The assumed neutral conditions (neutral composition and total neutral density) are those assessed by the ROSINA-COPS and ROSINA-DFMS (Balsiger et al. 2007) instruments, respectively (Heritier et al. 2017b). For this case, photo-absorption has an effect on the total ion density at small cometocentric distances (typically 20 km for this case) but has practically no effect at cometocentric distances larger than 100 km (comparing the purple and yellow curves). In particular, the effect of photo-absorption on the total ionospheric density at the location of Rosetta (380 km), for this typical near-perihelion day, is negligible.
During the whole escort phase, photo-absorption does not affect the plasma densities significantly at the location of Rosetta. Indeed, even when the spacecraft was located at lower cometocentric distances, the outgassing rate was accordingly lower (see Figure 1), therefore the effects of photo-absorption were negligible, even at 10 km. This is illustrated in Figure 5. The purple dots Article number, page 5 of 18 A&A proofs: manuscript no. plasmaSource represent the relative errors between the total ion densities at the location of Rosetta considering and neglecting photo-absorption, while the white circles represent the daily averages of these errors. We can see that neglecting photo-absorption yields an error of 10% at most. Consequently, it is safe to assume no photoabsorption of the neutral molecules when assessing the ionospheric densities at the location of the spacecraft. However, it would be wrong to say that photo-absorption has no appreciable effect at all on the ionospheric density profiles as we clearly see an induced error at small cometocentric distances (see Figure 4), when the outgassing rate is sufficiently high. Furthermore, we only considered volatile photo-absorption. There are currently some new investigations revealing that dust photo-absorption could affect the solar flux, and thus the ionospheric densities, in particular around perihelion (Johansson et al. 2017).

Electron-ion dissociative recombination
Electron-ion dissociative recombination (DR) is an ionospheric loss process that occurs when ions and electrons recombine to form dissociated neutrals: It takes place especially when the ionospheric densities are high and the electrons are efficiently cooled through collisions with cold neutrals. Figure 4 shows the impact of DR on the ionospheric densities (with DR in yellow, without DR in blue) as a function of the cometocentric distance for near-perihelion neutral conditions (31 August 2015). The DR kinetic coefficients were taken from Heritier et al. (2017b). The electron temperature is assumed to be 200 K. While the electron temperature is generally higher for the different electron populations in the coma Odelstad et al. 2018), it also strongly depends on the heliocentric and cometocentric distances. The use of 200 K (the neutral temperature) gives a lower bound that defines an upper bound for the electron-ion dissociative recombination effects. Indeed, the ion-electron dissociative recombination kinetic coefficients are proportional to the inverse square root of the electron temperature (McElroy et al. 2013;Heritier et al. 2017b). On 2015 August 31, the Rosetta spacecraft was located at 380 km from the nucleus. We can see that, at this cometocentric distance, neglecting DR yields a total ion density over-estimation of up to 50%. This relative over-estimation is increased at lower cometocentric distances. It is therefore not a valid assumption to neglect DR when assessing the cometary ionosphere subject to the high outgassing conditions encountered around perihelion.
When the outgassing rate is low, the number density of ions and electrons is reduced and reaction (3) occurs less frequently. We have assessed the periods when the outgassing and cometocentric conditions were such that we could neglect DR when computing the ionospheric densities at the position of the spacecraft during the escort phase of Rosetta. Figure 6 shows the relative error (blue) between considering and neglecting dissociative recombination in the local (at the position of Rosetta) ionospheric densities. The ion densities are computed according to the analytical model presented in Beth et al. (2018) (taking into account photo-ionization, ion transport, and DR, while assuming a constant outflow velocity of 600 m.s −1 and a constant electron temperature of 200 K). As in Figure 5, the white circles correspond to the daily average of these errors. Unlike photo-absorption, the impact of DR on the ionospheric densities is much more pronounced, and not only near perihelion. The dashed red line represents the 20% error threshold. From June 2015 (1.52 au) to December 2015 (1.78 au), the impact of DR cannot be neglected when assessing the ionospheric densities. Outside these periods, the error is less than 20%. Furthermore, the computed error is an upper bound. We expect the actual error to be even less than that computed for low outgassing conditions as we have under-estimated T e and considered all electrons to be cold. Based on this analysis, we have neglected DR when we computed the ionospheric density profiles post-perihelion, from April 2016 to September 2016 (see Section 3.3).

Post-perihelion ionizing sources
Since pre-perihelion cases were studied by Galand et al. (2016), in this section we focus on the ionospheric population postperihelion measured by the Rosetta instruments at large heliocentric distances (> 2.5 au). The ionospheric population at small heliocentric distances (near perihelion) will be the object of a separate study as not only dissociative recombination needs to be taken into account (with an accurate electron temperature, see Section 2.3) but also because solar flux becomes the main ionizer (see Section 4) and there is uncertainty on its local magnitude, due to possible attenuation through dust near perihelion (Johansson et al. 2017). We first give an introduction to the data-based model (Section 3.1) that is used. We then give a brief description of the instruments used to measure the total plasma densities (RPC-MIP and RPC-LAP) and the way the data has been processed (Section 3.2). Finally, we present five post-perihelion case studies that were chosen to assess the sources of variability of the ionospheric population as well as the robustness of the model (Section 3.3).

A data-based ionospheric model
We derived a model to assess the total ion density at the location of the Rosetta spacecraft using multi-instrument measurements. Post-perihelion, at large heliocentric distances (r h > 2.5 au) and small cometocentric distances (r < 20 km except for the April 2016 case study of Section 3.3.4, during which 20 km < r < 70 km), assumptions can be made to reduce the complexity of the ion continuity equation. Section 3.1.1 focuses on the continuity equation applied to the total number of ions. Section 3.1.2 justifies the assumptions made for the photo-ionization and electronimpact ionization frequency profiles. Finally, Section 3.1.3 depicts an analytical solution to the continuity equation.

Continuity equation
In Section 2.3 we showed that dissociative recombination, as an ion loss process, was negligible at large heliocentric distance. Under this assumption and assuming steady-state conditions and spherical symmetry, the total ion continuity equation can be written as (Galand et al. 2016;Heritier et al. 2017a) n i (r) [m −3 ] is the total ion number density at the location of the spacecraft; u i (r) [m.s −1 ] is the ion bulk velocity. In this study we assume that the ions keep the inertial velocity of the parent neutrals from which they originate. At very small cometocentric distances, in the densest part of the cometary atmosphere, there is evidence for outwards acceleration (Heritier et al. 2017a). However, this acceleration only takes place in the near-surface region. Assuming a constant velocity is therefore justified for the Rosetta cometocentric distances considered here (r > 7 km) as ions are not only produced near the surface but throughout the full [r 0 , r] range, where r 0 is the effective comet radius (2 km). At similar heliocentric distances but large cometocentric distances (> 100 km), ion transport is influenced by the solar wind magnetic field, enabling the ion pick-up process (Glassmeier 2017;Nilsson et al. 2017). The electrons originating from photo-ionization or electron-impact ionization typically have higher energies than the new-born ions and would tend to leave the ambient medium. Instead, they are retained by the ambipolar electric field through small-scale departure from charge neutrality, in order to ensure large-scale quasi-neutrality (Cravens 1997). At the same time, ions are expected to be accelerated outwards through the effect of this radial electric field. However, evidence for ion acceleration has as yet only been presented at relatively large cometocentric distances (> 100 km, see Vigren & Eriksson (2017); ; Odelstad et al. (2018)). We have therefore chosen to take a constant ion velocity (though within a given uncertainty range) in our study. This assumption was proved to be justified in Galand et al. (2016). The neutral (thus ion) outflow velocity is estimated from previous ionospheric studies (Heritier et al. 2017a;Galand et al. 2016) and measurements from MIRO (Gulkis et al. 2007), which give values lying within the 300 to 800 m.s −1 range Lee et al. 2015;Marshall et al. 2017) for large heliocentric distance (low activity) conditions. -ν e− (r) [s −1 ] is the electron-impact ionization frequency. It is estimated from in-situ measurements of RPC-IES (Burch et al. 2007). RPC-IES electron flux densities are retrieved from counts by using the in-flight calibrated geometric factor and the MicroChannel Plate (MCP) efficiencies applied to each energy bin. The fluxes are then integrated along azimuthal and elevation angles. We assume isotropy of the electron population to estimate the fluxes in the missing field of view. The negatively charged spacecraft affects the electron flux density measured by IES. It is, however, possible to partially retrieve the original electron flux density from the spacecraft potential measured by RPC-LAP using Liouville's theorem (Galand et al. 2016). This not only corrects the energy repartition of the ionizing electron population but also restores consistency in the dataset as the corrected frequencies are generally less noisy (Heritier et al. 2017a). , assuming a constant outflow velocity for each neutral species. It is therefore useful to apply a correction factor depending on the neutral composition. In the following model, this correction factor is concealed within the ionization frequencies (also affected by the neutral composition), as in Galand et al. (2016).

Ionization frequencies independent of r
Photo-ionization frequencies:~4 km~1 .6km dV z x y Fig. 8. In spherical coordinates, the ions produced in a volume element dV 1 (solid red sector) travel outwards, and the volume containing them expands with r 2 to reach a magnitude dV at r (dashed red sector). The neutral number density n decreases with r −2 (Haser 1957) so the ion production rate (per cubic meter) decreases but the number of ions produced in each element of constant solid angle is conserved.
In Section 2.2.2, we have assessed the effect of photoabsorption. Its effect on the total ion density at the location of the spacecraft was found to be insignificant throughout the two-year escort of comet 67P. The photo-ionization frequency ν hν (r) should therefore be independent of the cometocentric distance r. However, the comet nucleus is not a fully convex shape and some regions can end up in the shadow, out of direct reach of the solar EUV radiation (or, only reachable through reflection). In particular, the larger lobe can obstruct the near-surface neck region (see Figure 7).
During most of the escort phase, the Rosetta spacecraft was located in the terminator plane (as depicted in Figure 7). The ion density, measured at the location of the spacecraft, is driven by the ionization taking place along the path-line from the nucleus surface to the spacecraft (Equation 4).
With no obstruction, the ions reaching the spacecraft at a cometocentric distance r are equally likely to be produced close or far from the nucleus. Indeed, if we call r 1 and r 2 such that r 0 < r 1 < r 2 , at r 1 , the number of ions produced per unit time, within dV 1 (see Figure 8), is ν hν n(r 1 ) dV 1 . These ions are transported to r with constant velocity and expand to a volume dV (with dV > dV 1 , see Figure 8) located at r. The number of ions produced at r 1 and arriving at r, per unit volume and unit time is ν hν n(r 1 ) dV 1 /dV = ν hν n(r 1 )(r 1 /r) 2 . Similarly, the number of ions produced at r 2 and arriving at r, per unit volume and unit time is ν hν n(r 2 )(r 2 /r) 2 . Assuming Haser (1957) model, we have n(r 2 )(r 2 /r) 2 = n(r 1 )(r 1 /r) 2 , so their respective contributions are identical. In short, every cometocentric distance contributes uniformly to the number density of ions at a given range r, since the loss by transport affects ions and neutrals identically.
Consequently, if obstruction occurs along the particle trajectory, the maximum relative decrease of the observed n i is of the same order as the fraction of the trajectory that is shadowed. In a situation similar to Figure 7, the length of the obscured section is about 1.6 km. If the spacecraft is located at 10 km from the surface, n i at Rosetta is reduced by about 16%. That is only if we assume that all ions are produced by photo-ionization. As electron-impact ionization is also a significant contributor (Section 3.3), and is thought to be driven by ionizing electron fluxes coming from the space environment (Deca et al. 2017) -possibly undisturbed by the obstruction -the relative contribution of photo-ionization is reduced and the overall effect of obstruction on the ion density near the spacecraft is diminished.
This result also explains why photo-absorption has so little effect on the total ionospheric densities near perihelion. Indeed, despite the reduction of the photo-ionization frequency by a whole order of magnitude at a cometocentric distance of 10 km (Heritier et al. 2017b), the relative error resulting from neglecting this effect is less than 4% (see Figure 5). That is because the spacecraft was located at ranges of 170 to 400 km from the nucleus during this period. As the ions produced in the first 10 km do not contribute more than the rest of the ions produced along the full nucleus-spacecraft line, omitting these ions should only cause a reduction of about 2.5 to 5.9% on the total ion density at the location of Rosetta.
Furthermore, this effect is partially negated by the assumption of a spherical nucleus. In the model (Section 3.1.1), the cometary nucleus is assimilated to a sphere of 2 km radius. However, shadowing effects occur in priority around concave regions, where the distance between the nucleus centroid and the surface is actually much smaller than 2 km. By assuming that the comet is a sphere, we have already removed the contribution of any ion created below 2 km anyway. We therefore neglect any obstruction effect in the case studies covered in Section 3.3 and it is justified to assume that the photo-ionization frequency is independent of r.

Electron-impact ionization frequencies:
The ionizing electron flux can also interact with neutrals and lose significant energy as it penetrates into the coma, a process commonly referred to as energy degradation. However, the mechanisms governing the electron energy balance are much more complex than for photo-ionization. Unlike photons, electrons are not absorbed by neutrals during ionization (or other inelastic interactions). They give some of their energy to trigger the ionization process and to the newborn, secondary electron. After ionization, the primary electron can possibly ionize subsequent molecules, and so can the secondary electron, if their energies allow it.
As a first approximation, we assume that the energetic electrons detected by RPC-IES can only perform a single ionization and that secondary electrons do not ionize subsequent neutrals.
The energy lost from ionization is at least equal to the ionization energy threshold (12.62 eV for H 2 O, Itikawa & Mason (2005)). The ionization threshold is the energy needed to take a charge from the molecule without dissociating it (i.e. H 2 O → H 2 O + ). Ionization reactions that yield dissociated products (i.e. H 2 O → OH + +H) take even more energy from the primary electron as more energy is required to break molecular bonds. We further assume that during the process, the ionizing electron gives half of its energy to the newborn electron. That is somewhat arbitrary as in reality, the relative amount of energy given to the newborn electron is a stochastic random variable, which can vary between 0% to 100% of that remaining energy pool but the shape of the corresponding probability distribution function is poorly known (e.g. Lummerzheim & Lilensten (1994)). Under such a set of approximations, for a typical energy of 40 eV, as witnessed in the RPC-IES electron flux (Galand et al. 2016), the primary and secondary electrons would have an energy of 13.7 eV after the ionization reaction. That is almost below the H 2 O ionization energy threshold and the corresponding cross section is very low (σ e−,ioni H 2 O (13.7 eV) = 2.5 × 10 −18 cm 2 , Itikawa & Mason (2005)). Now, even if one of these two electrons has slightly more energy, say 20 eV, and the other 6.5 eV, the (total) ionization cross section at 20 eV is still significantly smaller than at 40 eV, as σ e−,ioni H 2 O (40 eV) > 3 × σ e−,ioni H 2 O (20 eV). So the single ionization approximation is not too unrealistic as, in the energy ranges encountered in the plasma surrounding 67P, the ionization capability of the electrons is greatly decreased after one ionization. We have neglected other inelastic electron-neutral interactions, such as dissociative excitation, which can also drain energy out of the ionizing electron fluxes. The dissociation energy is about 5 eV for H 2 O and the corresponding cross sections are of the same order as the ionization ones (Itikawa & Mason 2005). We can expect the ionization capability of an electron to be not too affected by such interaction, keeping in mind the modest loss in energy through dissociation, though a rigorous calculation would be required to validate this statement.
In that vein, we make the ionizing electron flux behave similarly to the energetic photon flux discussed in Section 2.2.2. In particular, the optical depth τ, defined in Equation 2, can be generalised to the electron depth τ e− (E) such that τ e− (E) = n r r 0 n n (s)σ e−,ioni n (E) ds.
We are integrating between r 0 and r as the ionizing electron flux is expected to be directed inward from the space environment (Deca et al. 2017). For a typical energy of 40 eV for the ionizing electron flux (Galand et al. 2016) and neutral conditions encountered at large heliocentric distances such as 3.1 au in October 2014 (Q = 10 26 s −1 at 500 m.s −1 ), we get τ e− (40 eV) = 0.12, from the surface to 10 km (typical position of Rosetta during large heliocentric distance escort). This is lower than 1, but that could still be temporarily increased by an order of magnitude during outgassing peaks. Near-perihelion, for extreme outgassing conditions of Q = 3 × 10 28 s −1 at 1000 m.s of the diamagnetic cavity, which was also found to greatly affect the ionizing electron distribution function (Madanian et al. 2017).
The case studies covered in Section 3.3 are located at large heliocentric distances. We assume no energy degradation and an electron-impact ionization frequency independent of r while keeping in mind that occasional energy degradation is not inconceivable during outgassing peaks. This assumption turns out to be valid for most cases, when we compare the results of the model to the plasma densities dataset (Section 3.3).

Analytical solution
Assuming that both photo-ionization and electron-impact ionization frequencies are constant with cometocentric distance, Equation 4 can be solved analytically (Galand et al. 2016): Despite its apparent simplicity, Equation 6 is extremely powerful. Under the Haser (1957) model and uniform ionization frequencies, it describes that the local ion density only depends on the local neutral density even though these ions are produced along the entire path from the nucleus surface to the point of interest. It also perfectly suits Rosetta multi-instrument in-situ measurements as depicted in the schematic of Figure 9. It is our main tool to understand the cometary ionosphere at large heliocentric distances across the many post-perihelion case studies covered in Section 3.3.

Instrumental observations from RPC-MIP and RPC-LAP
The modelled ion density is compared to RPC-MIP (Trotignon et al. 2007) and RPC-LAP (Eriksson et al. 2007) observations (see Figure 9). The RPC-Mutual Impedance Probe (MIP) allows us to derive the local plasma frequency from the mutual impedance spectrum between two electric antennae. From the plasma frequency, the in-situ electron number density is obtained and used for our comparison. The RPC-MIP instrument has two operating modes. The Short-Debye-Length (SDL) mode is designed to probe plasma frequencies within the 7 kHz to 3.5 MHz range and uses the RPC-MIP transmitting antennae separated by 40 and 60 cm from the RPC-MIP receiving antennae. It is used to probe the electron density in plasmas where the Debye Length (λ D ) is expected to be lower than half the distance to the closest receiving antenna (λ D <20 cm). When the electron density is too low to be detected by MIP in SDL mode, and this condition is no longer respected, RPC-MIP can be used in Long-Debye-Length (LDL) mode. This mode takes advantage of the RPC-LAP2 probe, used as a transmitter and located 4 m away from the RPC-MIP receiving antennae. Just like SDL mode, the range of plasma density that can be retrieved in LDL mode is limited. On one hand, if the electron density is too low, the Debye Length is too high (λ D > about 2 m), even for the LDL mode, and the instrument becomes blind to very low plasma densities. This typically happens for electron densities below 50 cm −3 . On the other hand, if the electron densities are too high, as this mode only covers plasma frequencies lying within the 7-168 kHz range, the plasma frequency is too high to be detected. This typically happens for electron densities above 300 cm −3 .
This last point is illustrated in Figure 10 where we plot the time evolution of the electron densities as probed by RPC-MIP from 2016 July 9 to 2016 July 11. The pink dots correspond to the instantaneous RPC-MIP measured electron densities. The instantaneous RPC-MIP electron densities exhibit significant variations on very short timescales because of high frequency plasma dynamics. Plasma dynamics are not included in the model presented in Section 3.1 so we are mostly interested in low frequency variations of RPC-MIP for comparison purposes. It is therefore usually practical to use moving averages when assessing the ionospheric densities as done previously in Galand et al. (2016) and Heritier et al. (2017a). The moving averages are performed during a ten-cycle RPC Plasma Interface Unit (PIU) (Carr et al. 2007) interval, which corresponds to roughly five minutes. On 2016 July 10 14:00 UT, as the local electron densities decreased, RPC-MIP was switched from SDL to LDL mode. However, the plasma frequency was higher than the 168 kHz threshold. This can be observed from the RPC-MIP mutual impedance spectra (not shown here) and can also be seen in the saturation of the measured electron densities at about 270 cm −3 (Figure 10, red dashed line). Due to the presence of an upper boundary, similar to what is experienced on 2016 July 10 14:00 UT, a standard moving average would be biased towards lower densities. When computing them at time t, we make sure that more than 70% of the data points used for the computation (i.e. the data points within the time range [t − 2.5 min : t + 2.5 min]) lie within the range of detectability. When these conditions are not fulfilled, the average is not computed and no value for the ionospheric density is provided. The moving averages are plotted in Figure 10 (purple dots). As expected, most of the moving averages between 2016 July 10 14:00 UT and 2016 July 11 00:00 UT cannot be computed due to these blindness effects.
Such blindness effects on RPC-MIP are common and it is useful to use RPC-LAP electron densities to fill the missing data gaps and build consistency within the existing RPC-MIP data. We use Langmuir probe bias voltage sweeps to determine dI/dU (Eriksson et al. 2007). Assuming a drifting Maxwellian distribution for the ion population, dI/dU is proportional to n i /u i (Jacobsen et al. 2009), where n i is the ion number density and u i is the effective ion velocity (combining thermal and drift velocities). The latter being unknown, we calibrate n i to n e as measured by RPC-MIP (in the regions where the moving averages are well defined) and apply a constant factor for the full period considered (a few days). While this factor is still quite uncertain and could vary within shorter timescales than those we have considered, this method still gives an idea of what should be the ionospheric densities in regions where RPC-MIP moving averages are uncertain. In particular, we can see in Figure 10 that the ionospheric densities estimated by RPC-LAP (green circles) are higher than the saturated RPC-MIP electron densities between 2016 July 10 14:00 UT and 2016 July 11 00:00 UT. This also strengthens the idea that performing a standard average of the RPC-MIP instantaneous electron densities would have resulted in an under-estimation of the ionospheric densities due to the boundaries of RPC-MIP operational modes.

Case studies
In this work, five different case studies have been performed. Each case study was selected for a specific purpose. First, we considered a post-perihelion case (May 2016), symmetric -in terms of heliocentric distance -to a case previously analysed by Galand et al. (2016) in October 2014, in order to illustrate the influence of seasonal effects on the cometary ionosphere. For simplicity, sub-solar latitudes are not discussed but should be included for a deeper analysis of the illumination conditions throughout the seasons. Second, we considered another post-perihelion case (March 2016) symmetric to October 2014, but this time in terms of local photo-ionization rates. Then, two cases (July and August 2016) were studied to show the effect of high variability and intensity in terms of local ionizing fluxes. They illustrate the effects of solar wind transients, such as co-rotating interacting regions (CIRs), onto cometary ionospheric densities. Finally, the return from the night-side excursion of April 2016 was selected as a final case, as it features strong variations in cometocentric distances and demonstrates the robustness of the ionospheric model. These different case studies and their respective purposes are summarised in Table 2. A last case (April 2016) was selected on the return from a night side excursion, at higher cometocentric distances, and demonstrates the robustness of the ionospheric model to significant variations in terms of cometocentric distances. These different case studies and their respective purposes are summarised in Table 2. Case r h Sec. Purpose May 2016 3.1 3.3.1 Symmetric to pre-perihelion October 2014 study performed by Galand et al. (2016) in terms of heliocentric distance Mar 2016 2.6 3.3.2 Symmetric to pre-perihelion October 2014 study performed by Galand et al. (2016) in terms of local photo-ionization frequency Jul 2016 3.4 3.3.3 Effect of a CIR on the electronimpact ionization frequency and the induced ionospheric population Aug 2016 3.5 3.3.3 "Return" of the July 2016 CIR, after one solar rotation, and illustration of a case of high plasma densities purely dependent on the electron-impact ionization. Apr 2016 2.8 3.3.4 Illustration of the robustness of the model to high variations in cometocentric distance undergone during the return from the cometary tail excursion. Table 2. Summary of cases studied in Section 3.3 and their respective purposes.

2016 May 25-29
The first case study is a post-perihelion symmetric case with respect to the pre-perihelion October 2014 period analysed by Galand et al. (2016). For the latter, the heliocentric distance was about 3.1 au and the cometocentric distance was rather small (about 10 km). Galand et al. (2016) found strong seasonal variations within the ionospheric population and their different sources. Over the summer (northern) hemisphere, where the outgassing rate was the highest, photo-ionization was found to be the dominant source of ions. However, over the winter (southern) hemisphere, electron-impact ionization was found to be the dominant source. These seasonal effects are generally manifested by the anti-correlation of the ionizing electron flux with respect to the neutral densities. One of the reasons for this can be energy degradation through electron-neutral collisions. Our calculation performed in Section 3.1.2 provided an electron depth of only 0.12 for such conditions but the somewhat simplistic assumptions behind it leave that possibility open, especially during outgassing peaks. Through electron-neutral inelastic collisions, the ionizing electrons could be subject to energy loss, driving their ionizing capability down, therefore letting unattenuated EUV radiation be the dominant ionizer. We wished to investigate if the results and seasonal effects described Average heliocentric distance, in au, for the period studied.
For the post-perihelion case analysed here, the spacecraft was back at 3.1 au from the Sun, at a similar cometocentric distance of about 8 km. The results of the ionospheric model are plotted in Figure 11. The top panel corresponds to the time evolution of the total neutral number density (solid, black curve) plotted together with the cometocentric distance (dashed, black curve). The seasons are colour coded in the top panel, with spring in the northern hemisphere (positive latitudes) and autumn and the southern hemisphere (negative latitudes), happening shortly after equinox (see Table 1). The neutral densities (and the neutral outgassing, since we are at constant cometocentric distances and assume constant neutral outflow velocity) are lower in spring than in autumn. This counter-intuitive result is partially due to the heterogeneity of the nucleus sub-surface (Gasc et al. 2017a;Fougere et al. 2016a). While the southern hemisphere (autumn) is cooler, the mixing ratio of CO2 in the south is likely higher (in the sublimating layer) than in the north, and it possesses a lower sublimation temperature than H 2 O.
The middle panel shows the photo-ionization (blue) and electron-impact ionization (red dots) frequencies. Despite the seasonal singularity explained above, the results are consistent with Galand et al. (2016). High neutral densities contribute to lowering the ionizing flux and we observe an anti-correlation between the neutral densities and the electron-impact ionization frequencies.
The bottom panel displays the results of the data-based ionospheric model presented in Section 3.1 compared with RPC-MIP measurements. We plotted the estimated ionospheric population using the data-based model through photo-ionization only (blue ribbon), and using both photo-ionization and electron-impact ionization (red ribbon). The width of the ribbons illustrates the uncertainty in the outflow velocity. The upper bound of each ribbon corresponds to the slow flow case of 300 m.s −1 and the lower bound corresponds to the fast flow case of 600 m.s −1 along the whole column from the surface to the spacecraft. These assumed outflow velocities are consistent with the observations of MIRO neutral outflow velocities (Marshall et al. 2017;Gulkis et al. 2015;Lee et al. 2015). The results of the ionospheric model are compared with RPC-MIP instantaneous electron density measurements (pink). These measurements are smoothed using a five-minute moving average excluding blindness effects (purple), as explained in Section 3.2. We have also plotted RPC-LAP ionospheric densities, cross calibrated with RPC-MIP as described in Section 3.2 (green).
There is an excellent agreement between the modelled densities and the plasma density measurements. We can see that the ionospheric population is overall higher in autumn where the neutral densities are higher. While a positive correlation between neutral and ionospheric populations can appear trivial as more neutrals mean more sources, it was not the case in October 2014 where low neutral number densities were associated with high energetic electron fluxes yielding high ion number densities.
Photo-ionization alone (blue ribbon) is not able to fully explain the observed ionospheric densities, even when the neutral densities are high. Electron-impact ionization seems to be the dominant source (or, at least, of similar importance) over both hemispheres. The local maxima and minima in plasma density observed by RPC-MIP and RPC-LAP are only captured if electron-impact is taken into account. Before perihelion, at the same heliocentric distance, we could explain the ionospheric densities in summer using photo-ionization only, so why would the situation be different after perihelion? The main difference between this case and its pre-perihelion counterpart is that the solar flux has decreased (see Figure 2) and therefore the photoionization frequency (see Figure 3) at a given heliocentric distance has decreased with time. Another effect could be the different illumination conditions throughout the regions of the comet (Fougere et al. 2016b). To assess a true symmetric case with respect to 2014 October 17-20, we need to find a period where the photo-ionization frequency was similar at the location of the comet. To do this, we need to go back in time to March 2016 (see Section 3.3.2) as a reduced heliocentric distance can compensate for the decrease in the solar flux, the solar flux being inversely proportional to the square of the heliocentric distance.

2016 March 18-23
In March 2016, the heliocentric distance was 2.6 au and the cometocentric distance (about 12 km) was close and comparable to October 2014. The photo-ionization rate at 1 au was lower than in 2014 due to the variations of the solar cycle (see Figure 2). If we correct for the heliocentric distance, the local photo-ionization rate in March 2016 was similar to the local photo-ionization rate in October 2014 (on average ν hν 7 × 10 −8 s −1 , see Figure 12 and Galand et al. (2016)). Figure 12 shows the results of the data-based model for the 2016 March 18-23 period. The colour coding is similar to Figure 12 but it must be highlighted that this case begins just before equinox: the northern (winter) hemisphere in cyan and the southern (summer) hemisphere in red are shown in the top panel. After the equinox on 2016 March 21, the northern (spring) hemisphere is shown in pink and the southern (autumn) hemisphere is shown in yellow. The seasonal variations are apparent: the summer hemisphere is associated with higher neutral number densities than any other season. An interruption of measurements in the neutral densities is visible at the beginning of 2016 March 19. This is due to a spacecraft manoeuvre.
The middle panel compares the different ionization frequencies while the bottom panel compares the ionospheric densities observed and modelled. There is still an excellent agreement between the modelled densities and the RPC-MIP/LAP observations. Photo-ionization is the dominant source during the summer (March 19 and 20) and is similar or weaker than electronimpact ionization during the rest of the period. In that vein, the results are more similar to the case study of 2014 October 17-20 without being identical. Indeed, just like for 2016 May 25-28, the highest ionospheric densities occur when the neutral densities are the highest (here summer). This was not the case in October 2014 (especially October 18 and 19, 2014, see Galand et al. (2016)). Even if the energetic electron flux seems to be stronger with low neutral densities, it is not strong enough to yield ionospheric densities higher than that observed over the summer hemisphere.

Summer 2016: Electron-impact ionization and CIRs
As the mission progressed, EUV radiation became less and less intense due to the decrease in the solar flux (Section 2.2.1). In the last three months of the mission, the relative importance of electron-impact ionization with respect to photo-ionization was high and very variable (see Section 4). Occasionally, it was so high that photo-ionization can be omitted when modelling the total ionospheric densities. Figures 13 and 14 show the results of the data-based model for the 2016 July 03-17 and 2016 August 03-10 periods, respectively. The colour coding is similar to Figure 11 and 12 with the exception of the middle panel, on which we have added the estimated solar wind dynamic pressure (green) based on the Tao et al. (2005) model. During these two time windows, the spacecraft was performing consistent elliptic trajectories with an apoapsis (respectively, periapsis) of 28 km (resp. 10 km) for the July period and an apoapsis (resp. periapsis) of 12 km (resp. 8 km) for the August period of study. The total neutral densities monitored by ROSINA-COPS are consistently anti-correlated to the cometocentric distances covered, while featuring typical (12 hour) diurnal variations (Hässig et al. 2015).
We can see that the ion number densities assuming photo-ionization only (in blue, bottom panel) are nowhere near the measured RPC-MIP/LAP plasma densities. Indeed, electron-impact ionization frequencies (Figures 13 and 14, middle panels) are often an order of magnitude higher than photo-ionization frequencies, though they are very variable over time. This variability cannot be explained either by the neutral densities, or the geometrical conditions. The orbits covered the same range of cometocentric distances (Figures 13  and 14, top panels) twice in a row in a periodic manner, along with longitudinal and latitudinal conditions (not shown here). However, the energetic electron fluxes were not periodic at all. For instance, the electron-impact ionization frequency encountered at Rosetta 10 km periapses on July 4 and July 10 were about 6.44 (±2.22) × 10 −9 s −1 and 1.52 (±1.07) × 10 −6 s −1 (see Figure 13), respectively. Similarly, between 2016 August 3 and 2016 August 4, the electron-impact ionization increased from about 5.43 (±6.27) × 10 −8 s −1 to 9.04 (±9.71) × 10 −7 s −1 (see Figure 14) whereas the neutral and geometrical conditions barely changed over these two days.
These irregularities can be explained by the presence of a corotating interaction region (CIR) in the solar wind (Smith & Wolfe 1976;Pizzo 1985;Balogh et al. 1999). This CIR was observed four times in the vicinity of comet 67P: June 13-15, July 6-11, August 3-6, and September 1-3. The average 28-day period corresponds to the sidereal rotation of the Sun. This CIR is analysed in detail in Hajra et al. (2018). We can point out that its signature is observable within the solar wind dynamic pressure estimates (Figures 13 and 14, middle panels, in green). The increase in the solar dynamic pressure is correlated with an increase of the energetic electron flux. The peak in solar wind dynamic pressure on July 10, 00:00 UT is correlated with a peak in RPC-MIP, completed by RPC-LAP due to the upper boundary of RPC-MIP SDL mode (see Section 3.2), but lagging with respect to the peak in the electron-impact ionization frequency (July 10, ∼12:00 UT). This can be related to the neutral environment. The peak in neutral density occurs at 00:00 UT (Figure 13, top panel) and is essential for the plasma density peak. However, high neutral densities are often Seasonal conditions are presented with the following colour code: cyan area for winter (northern hemisphere before equinox), red area for summer (southern hemisphere before equinox), yellow for autumn (southern hemisphere after equinox), and pink for spring (northern hemisphere after equinox).  Figure 11 applied to 2016 July 3-17. The neutral outflow velocities are assumed to lie within 400 to 700 m.s −1 in the ionospheric model. We have also plotted the solar wind dynamic pressure (green curve, middle panel) as predicted by the Tao et al. (2005) model to illustrate the effect of the CIR impact. associated with weaker ionizing electron fluxes (see Section 4 and Galand et al. (2016)). We note that when the peak in solar wind dynamic pressure occurs at low neutral densities, such as on July 6, there is no lag between the increase of the modelled solar wind dynamic pressure and the increase in electron-impact ionization frequency.
These case studies from summer 2016 show that the databased ionospheric model (Section 3.1) is robust to solar events due to the continuous stream of data coming from RPC-IES and provide further evidence of the strong anti-correlation between the ionizing electron flux and the neutral environment.

2016 April 17-27
In addition to being robust to singular events such as CIRs increasing the local energetic electron flux (see Section 3.3.3), the model is robust to significant variations in cometocentric distances. Figure 15 shows the time evolution of the ionospheric densities (modelled and measured, bottom panel), ionization frequencies (middle panel), and neutral conditions (top panel) from 2016 April 17 to 2016 April 27. These measurements occurred on the return from the tail excursion between 2016 March 23 and 2016 April 09. They spread from a large cometocentric distance (70 km) down to a small cometocentric distance of 20 km, with ionospheric densities increasing from about 100 cm −3 to 1000 cm −3 .
During the descent, the local total neutral number density measured by ROSINA-COPS increased roughly with the inverse square root of the decreasing cometocentric distance (Figure 15, top panel) while undergoing local diurnal variations. The photo-ionization frequency stayed constant throughout the descent but the electron-impact ionization frequency underwent significant variations within 10 −8 -10 −6 s −1 . Despite these variations, the modelled densities remain consistent with RPC-MIP and RPC-LAP measurements and demonstrate the robustness of the data-based model. Once again an anti-correlation of the electron-impact ionization frequency with respect to the neutral density can be observed: on 2016 April 26, diurnal variations of the electron-impact ionization frequencies are visible, but in opposite phase with respect to the neutral diurnal variations.

Overview of the ionization sources
In this section, we compare photo-ionization and electronimpact ionization throughout the full escort phase. Figure 16 shows the time evolution of the H 2 O photo-ionization frequency at the location of comet 67P (blue curve, bottom panel). The photo-ionization frequency has been computed using TIMED-SEE (L3 V12.0) measured EUV solar fluxes (see Section 2.2 for details). We have neglected photo-absorption as it does not have a significant effect on the ionospheric population at the range of cometocentric distances covered by Rosetta (see Section 2.2.2). Electron-impact ionization frequencies for H 2 O have been computed during 82 days of the escort phases using RPC-IES local ionizing electron fluxes (see Galand et al. (2016);Heritier et al. (2017a) for details). They are represented by the red dots in Figure 16; the error bars show the standard deviations of the variations for each day.
We only show the ionization frequencies for H 2 O as it was the most abundant neutral species during most of the mission but the same trends are observed for any other neutral species. In the case of comet 67P, the main factors influencing the ionization frequencies are rather more the solar flux, ionizing electron flux intensities, and energy distributions (visibly affected by the neutral outgassing rate) than the local differences in the neutral composition of the gas mixture.

Ionization sources at large heliocentric distances
At large heliocentric distances, electron-impact ionization frequencies are generally equal to or larger than the photoionization frequencies. Galand et al. (2016) pointed out that in October 2014, electron-impact ionization was larger than photo-ionization over the southern (winter) hemisphere but was of the same order of magnitude (or even lower on the 2014 October 17) over the northern (summer) hemisphere. The same seasonal effects were observed in May 2016 (see Section 3.3.1) and March 2016 (see Section 3.3.2) where the electron-impact ionization frequency was anti-correlated with the neutral density. As the comet retreated to even larger cometocentric distances, electron-impact ionization became the dominant source on both hemispheres. Photoionization became negligible as both heliocentric distance and decreasing solar activity (see Figure 2) contributed to its decline.
The standard deviation of the daily variations of electronimpact ionization is quite significant at large heliocentric distances. When the spacecraft covers both hemispheres within the day, seasonal variations contribute to these variations. For instance, on 2016 May 26 (see Figure 11), the electron-impact ionization frequency decreased by about an order of magnitude (from ∼ 4 × 10 −7 s −1 to ∼ 4 × 10 −8 s −1 ) as the spacecraft moved from the northern (spring) hemisphere to the southern (autumn) hemisphere. These variations are also sensitive to singular events like the CIR encountered on 2016 July 9-10 and 2016 August 05-06 (see Section 3.3.3). In the case of July 2016 (see Figure 13), the electron-impact ionization frequency increases by an order of magnitude (from ∼ 2 × 10 −8 s −1 to ∼ 2 × 10 −7 s −1 ) at the impact of the CIR on July 9. Overall, it seems that at large heliocentric distances, and therefore low outgassing activity, the comet is less shielded from the external space environment and solar events can strongly affect its existing ionosphere.

Ionization sources near perihelion
As the solar flux intensity is proportional to the inverse square of the heliocentric distance, the photo-ionization frequency peaks near perihelion (2015 August 13) to reach values of the order of 4 × 10 −7 s −1 . Near perihelion, photo-ionization dominates over electron-impact ionization. While photo-ionization is high due to the proximity of the Sun, we note that electronimpact ionization frequencies measured in the vicinity of the spacecraft are rather low (around 5 × 10 −8 s −1 ) when compared to the values reached during the rest of the two-year escort phase.
It is not surprising to witness low electron-impact ionization frequencies near perihelion. Indeed, in the previous case studies (see Section 3.3), as well as in Galand et al. (2016), we saw anti-correlation between the electron-impact ionization frequency and the neutral activity, which is at its peak near perihelion. The electron depth computed in Section 3.1.2 for near-perihelion conditions is found to be three, which indicates significant energy degradation of the ionizing electron flux between Rosetta and the nucleus surface. Depending on where the ionizing electron flux is actually coming from, this could potentially attenuate the flux measured at the spacecraft location.
The standard deviations of the daily variations of electronimpact ionization frequencies are very low near perihelion but this may be due to orbital geometry. We note that daily variations at large heliocentric distances were mostly attributed to seasonal changes (i.e. changes of hemisphere). While at large heliocentric distances, the spacecraft followed bound orbits that covered both hemispheres in one day, this was not the case near perihelion where the spacecraft remained in each hemisphere for about 15 consecutive days (see Figure 1) and therefore these seasonal variations are not discernible at the daily timescale.

Conclusions
The Rosetta mission provided us with a unique opportunity to characterise the variability of a cometary ionosphere as it approached and retreated from the Sun. The main and sufficient sources that explain the observed plasma densities in the vicinity of comet 67P are photo-ionization and electron-impact K. L. Heritier et al.: Plasma source and loss at comet 67P during the Rosetta mission ionization.
Photo-ionization can be attenuated through photo-absorption by neutral molecules (Section 2.2.2) but the effect on the electron density at the location of Rosetta is minor. Indeed, at large heliocentric distances, the flux of ionizing photons remains practically unattenuated. At small heliocentric distances, the photo-ionization frequency can be decreased by an order of magnitude at small cometocentric distances (<10-20 km, see Figure 4). However, as ions are produced along the entire path from the comet to the spacecraft, this does not affect the ionospheric density significantly at the location of Rosetta. Moreover, the effect of photo-absorption on the ionospheric densities was always less than 10% throughout the two-year escort phase (see Figure 5).
Dissociative recombination is another phenomenon that affects the ionospheric densities. It is found to have a significant effect at small heliocentric distances, when the comet is near its peak of activity. Consequently, the error made by neglecting this phenomenon is higher than 20% (and can reach up to 150%) at the location of the spacecraft between June and December, 2015 (see Figure 6); that is, assuming an electron population that is cooled by the neutral population (which gives an upper bound for the error). With this effect in mind, we performed quantitative studies of the ionospheric population at large heliocentric distances where dissociative recombination is negligible and where we can use the ionospheric data-based model (Section 3.1) previously introduced in Galand et al. (2016) and Heritier et al. (2017a).
The post-perihelion counterpart, in terms of heliocentric ditances, of the case studied by Galand et al. (2016) (2014 is 2016 May 25-29 (Section 3.3.1). However, as the solar activity decreased throughout the mission (see Figure 2), 2016 March 18-13 (Section 3.3.2) turns out to be subject to the same conditions in terms of EUV radiation intensity and photo-ionization frequency. Similarly to October 2014, these cases feature seasonal variations that affect both neutral and ion populations. Over the colder, less illuminated hemisphere, electron-impact ionization is dominant with respect to photo-ionization whereas over the warmer hemisphere, the two ionization processes are almost equally important. As the solar activity continued to decrease during the post-perihelion phase, electron-impact ionization became increasingly significant (Heritier et al. 2017a). In summer 2016, a CIR hitting the cometary ionosphere was observable every solar rotation (∼28 days). In 2016 July 10 and August 4 (Section 3.3.3), it is found to strongly increase the ionizing electron flux. Consequently, the electron-impact ionization frequency increased by about an order of magnitude and drove the ionospheric population to increase accordingly (see Figures 13 and 14). As the ionospheric model is data-dependent, it is very sensitive to singular events, such as CIRs, or changes in cometocentric distances (Section 3.3.4) and is perfectly able to capture the observed ionospheric density, building consistency between instruments and giving us an unique insight into the ionospheric drivers. Since electron-impact ionization has been shown to be a major source for the cometary ionosphere, our results emphasise the need to better understand the physical mechanism responsible for the generation of suprathermal electrons observed at comet 67P/CG. So far, two mechanisms have been suggested to explain the observed ionizing suprathermal electrons: heating by lower hybrid waves through wave-particle interactions (Broiles et al. 2016) and parallel acceleration by the cometary ambipolar electric field (Madanian et al. 2016). Full kinetic numerical simulations (Deca et al. 2017) should help to disentangle the relative contribution of these two mechanisms.
An in-depth, near-perihelion, quantitative study of the ionospheric population as seen by the Rosetta instruments would complete the full picture of this two-year cometary ionosphere story. As mentioned in Section 4, photo-ionization is known to be the main source of these ions. Furthermore, their chemical