Plasma distribution around Comet 67P in the last month of the Rosetta mission

Abstract After accompanying comet 67P/Churyumov–Gerasimenko on its journey around the Sun and observing the evolution of its induced magnetosphere throughout the comet’s life-cycle, the Rosetta operations concluded at the end of September 2016 with a controlled impact on the cometary nucleus. At that time, the comet was located more than 3.7 AU from the Sun, but the data still show clear indications of a weak but well developed plasma environment around the nucleus. Rosetta observed this fading cometary magnetosphere along multiple recurring elliptical orbits, which allow us to investigate its properties and spatial structure. We examined the measured electron and neutral densities along these consecutive orbits, from which we were able to determine the structure of the spatial plasma distribution using a simple latitude and longitude dependent model.


Introduction
At 3.6 AU from the Sun, on 6 August 2014, the Rosetta spacecraft  has rendezvoused with comet 67P/ Churyumov-Gerasimenko (67P) (Churyumov and Gerasimenko, 1972) and began to monitor its nascent atmosphere as the comet travelled towards its perihelion. The Jupiter-family comet 67P currently has a 6.44 years long orbit around the Sun with an aphelion distance of 5.68 AU and a perihelion distance of 1.24 AU. After accompanying 67P on its journey and observing the evolution of its plasma environment throughout the comet's life-cycle for more than two years, the operations of the Rosetta orbiter concluded on 30 September 2016, at 3.8 AU from the Sun, with a controlled impact on the cometary nucleus. Throughout these two years, the ESA Rosetta mission collected a variety of measurements that provide an immense insight into cometary physics.
Nearing perihelion, the activity of comets rises, and the neutral coma expands Biver et al., 2019). The large number of neutral particles are continuously ionized by photoionization, electron impact ionization and charge exchange with solar wind ions (Mendis et al., 1985;Cravens, 1991;Vigren et al., 2015;Galand et al., 2016;Madanian et al., 2016;Wedlund et al., 2017;Heritier et al., 2018). During the evolution of the cometary coma of 67P, photoionization and electron impact ionization were both shown to be necessary * Corresponding author.
to explain the observed electron densities over the southern, winter hemisphere while over the illuminated, northern hemisphere photoionization alone was reported to dominate the ionization processes Vigren et al., 2016). After perihelion, at large heliospheric distances (2 AU), electron impact ionization dominated over photoionization and was predominant during the last 4 months of the mission on both the southern and the northern hemispheres (Heritier et al., 2018).
An early sign of the cometary plasma environment around comet 67P was observed by  through the detection of water ions in the coma on 7 August 2014. At this time, the comet was located 3.6 AU from the Sun and the comet-spacecraft distance was approximately 100 km. The newly created heavy cometary ions are accelerated by the solar wind convective electric field and are picked up by the solar wind flow. As a result of the mass loading of the solar wind with cometary ions, the solar wind suffers an energy loss and is slowed down, piled up and deflected upstream of the comet (Coates, 1997;Szegö et al., 2000) although this close to the nucleus the spacecraft detected only the beginning of the mass loading process, apparent in the deflection of the solar wind ions .
During early activity, the high density plasma in the inner coma was investigated by Yang et al. (2016) who found that comet 67P's https://doi.org/10.1016/j.icarus.2020.113924 Received 4 April 2020; Received in revised form 27 May 2020; Accepted 10 June 2020 early plasma environment at a heliocentric distance of 3.4 AU consisted of two regions: an outer part mostly dominated by the solar wind convection electric field and an inner region of enhanced plasma density.
The evolution of the cometary ion environment was described during early activity in 2014 as the heliocentric distance decreased from 3.6 to 2.0 AU  as well as throughout the entirety of the mission . As the activity of the comet increased, the accelerated cometary ions became more common and reached higher energies. In April 2015, the solar wind disappeared from the vicinity of Rosetta -a solar wind cavity formed around the cometary nucleus . Inside the boundary called cometopause, the ion composition changes from a mixture of cometary and solar wind ions to picked-up cometary ions .
In the coma of comet 67P, at relatively large heliocentric distances (2.5 AU), the ion densities fall off with the radial distance from the comet with approximately r −1 based on both photochemical equilibrium and transport dominant models Vigren et al., 2016). In a model, presented by Nemeth (2020), in addition to the transport, production and loss, the effects of the magnetic field gradients were also taken into account. It was shown, that even in the presence of strong magnetic field gradients the plasma density features a r −1 radial dependence, except in the immediate vicinity of the diamagnetic cavity boundary. Edberg et al. (2015) reported a r −1 dependence of the electron densities based on Rosetta measurements performed in early 2015 within 260 km from the nucleus. These results also agree with the observations made at comet 1P/Halley during the Giotto mission (Cravens, 1987).
This observed vertical cometary density profile has been confirmed down to about 3 km from the nucleus surface with the observations made on the last day of operations (30 September 2016), during the controlled descend of the Rosetta orbiter (Heritier et al., 2018), using the combined measurements of the Mutual Impedance Probe (RPC MIP)  and the Langmuir Probe (RPC LAP)  instruments of the Rosetta Plasma Consortium (Carr et al., 2007). The findings were in a close agreement with cometary vertical ionosphere models predicting a maximum in the ionospheric densities close to the surface (Vigren and Galand, 2013) and a sharp decrease below this ionospheric peak .
Rosetta offers the unique opportunity to observe the fading cometary plasma environment in September 2016 through several similar, consecutive orbits. Our aim in this paper is to map the plasma environment around the nucleus of comet 67P through the electron densities measured by the RPC MIP experiment during the last month of the Rosetta mission. Our findings are explained and summarized by a distance, latitude and longitude dependent model of the plasma density of comet 67P.

Data
We investigated the spatial distribution of the cometary plasma around comet 67P in September 2016, more than one year after perihelion. At that time the comet was located at 3.7-3.8 AU, with sub-solar latitudes around 18-20 • on the northern hemisphere (Preusker et al., 2017). The Rosetta spacecraft had a highly elliptical orbit at 4-17 km from the nucleus with periods of approximately 3 days (Fig. 1). The nucleus had a rotation rate of 12.4 h. The top panel of the figure shows the trajectory in comet-centred solar equatorial (CSEQ) coordinates (the +X axis points from centre of mass of the nucleus towards the Sun, the +Z axis is the component of the Sun's north pole of date orthogonal to the +X axis, the +Y axis completes the right-handed reference frame). The bottom panel uses the body-fixed 67P/C-G_CK coordinate frame. (The origin of the frame is located in the centre of the comet, the +X axis points towards the prime meridian, the +Z axis towards the north pole while the +Y axis completes the right handed frame.) During this month, Rosetta performed eight very similar, consecutive orbits around comet 67P, suitable to perform a comprehensive 3D mapping of the cometary ionosphere.
We show the total electron density measured by RPC MIP and the total neutral density as measured by the Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA) on Fig. 2. The ROSINA instrument contains two sensors to determine the composition of the comet's atmosphere and ionosphere, the velocities of electrified gas particles, and reactions in which they take part (Balsiger et al., 2007). The main objective of the MIP experiment is to provide in situ electron density and temperature in the inner coma of 67P through the measurement of the mutual impedance between two electric dipoles embedded within the plasma to be investigated . The MIP sensor is made of two receiving and two transmitting electrodes, mounted on a 1 m long bar, itself mounted on a boom on the Rosetta orbiter. The instrument is capable of measuring plasma properties in two different operational modes. First, the so-called ''Short Debye Length'' mode (SDL), uses different combinations of a single or of the two MIP transmitters to access dense enough plasmas. Second, the so-called ''Long Debye length mode'' (LDL) uses the spherical probe of the LAP experiment, mounted on another boom and located 4 m from the MIP Z. Nemeth et al. antenna, as a monopolar transmitter. This LDL mode has been used to access lower electron densities than those accessible with the SDL mode, down to a few tens of cm −3 . During September 2016, MIP operated essentially in SDL mode, measuring densities up to thousands of cm −3 . The uncertainty of the measured electron density is estimated to be around 10%. The uncertainty is obtained from (i) the frequency discretization of the mutual impedance spectra and (ii) the stiffness of the spectral resonant signal in the mutual impedance spectra used to retrieve the electron density. Detailed information on the computation of the RPC-MIP plasma density uncertainty -as well as explanations on possible data gaps in the electron density due to the used MIP operation mode -is given in the RPC-MIP User guide  and reference therein, available in the Planetary Science Archive RPC-MIP archive .
Since the length scales over which we study the density distribution is much larger than the Debye length, we assume quasi-neutrality and take the MIP electron density results as a measure of the overall plasma density. Second, since the solar wind density at 3.8 AU is much smaller than the plasma density measured by MIP around 67P, we assume these measurements correspond to the overall cometary plasma density.
The plasma and neutral density curves in Fig. 2 feature clear periodicity corresponding to the orbital period of the spacecraft, but the signal is complex, not at all symmetric around the position of the closest approach to the nucleus. In addition to the main recurring peak, the data also show recurring fine structure. On the top and middle panels of Fig. 2 we also show the spacecraft's radial distance from the nucleus and its latitude and longitude in the body-fixed 67P/C-G_CK coordinate frame.
On the first days of September 2016 a corotating interaction region (CIR) impacted on the comet and disrupted the measured electron densities . In order to concentrate on the unperturbed cometary plasma, the present investigation focuses on the measurements from 4 September 2016 to 24 September 2016 (Fig. 2), before the spacecraft manoeuvred itself to collision course with the cometary nucleus.
By investigating the position of the measurements with respect to the surface of the nucleus, we can conclude that both the measured electron and neutral densities show a maximum at the southern hemisphere, after that the density falls off rapidly shortly before the spacecraft enters the northern hemisphere. On the top panel of Fig. 3, we show a projection of the trajectories onto the terminator plane in comet-centred solar equatorial (CSEQ) coordinates. We observe that the higher density measurements occur when the spacecraft is close to the nucleus, but the high density region is offset towards the negative z region. The bottom panel of Fig. 3 shows the data in comet fixed coordinates; the southern hemisphere clearly dominates. Although at this time the subsolar point is located at the northern hemisphere, the active regions (for cometary neutral production) were reported to be above the southern hemisphere during this period. Hansen et al. (2016) presented water distribution around the nucleus at 1.5 AU after perihelion with a maximum above the southern hemisphere, around latitudes −30 • . Kramer et al. (2017) showed how the highest neutral density regions 100 km above the nucleus shift from the northern to the southern hemisphere between April 2015 and May 2016. In May 2016, the highest density regions were above latitudes around −60 • and longitudes of −10 • . Combi et al. (2020) investigated H 2 O, CO 2 , CO and O 2 production rates throughout the Rosetta mission. They confirmed the H 2 O production rates to be dominant during the mission, except from mid-2016 where CO 2 gradually became dominant over all other species, its activity peaking at the southern hemisphere. As the main source of cometary plasma is the neutral outgassing of the nucleus, a strong correlation between the neutral and electron densities is expected.

Model and discussion
Figs. 2 and 3 show that although the radial distance plays an important role in the plasma density variation, it cannot be the sole player responsible for the observed structures. It is a reasonable hypothesis that the plasma density depends on the latitude and longitude coordinates in comet fixed frame. This hypothesis is supported by earlier results, e.g. Hansen et al. (2016) has shown that the neutral density features such angle dependence. The strongly non-spherical shape of the comet nucleus (Preusker et al., 2015;Jorda et al., 2016) and the solar-wind comet interactions (Deca et al., 2017(Deca et al., , 2019Koenders et al., 2016;Huang et al., 2016Huang et al., , 2018 can also influence the density distribution. In this section, we aim at providing a distance, latitude and longitude dependent model of the plasma density of comet 67P, which is able to reproduce the observed cometary data. Since for these highly excentric trajectories the vicinity of closest approach is associated to a fast latitude scan, it is possible that the rapid change in latitude is responsible for the drastic variation (strongest peaks followed by very low densities in Fig. 2) found close to the nucleus. Fig. 3 qualitatively supports this hypothesis. In addition to the highly apparent slow periodicity, the data in Fig. 2 also shows fine structures (secondary and sometimes higher order peaks before the main peaks for each orbit, see e.g. Sept. 8, 11, 14 and 17 in Fig. 2). These seem to follow the rotation period of the nucleus, which suggests that the plasma distribution may be best modelled in a comet fixed coordinate system. Thus, we modelled the 3D spatial distribution of cometary electrons and plasma around comet 67P in September 2016 in comet fixed spherical polar coordinates. We fitted by least squares method the parameters of the following test function ( , , ) to the in situ measured electron densities: where is the distance from the comet, is a constant corresponding to the angle averaged mean electron density on a hypothetical spherical source surface one kilometre over the centre of the comet. The angles and are the latitude and longitude of the spacecraft in the comet-fixed 67P/C-G_CK frame. This is the simplest possible expression, which describes a smooth partial angle dependence for both angle coordinates together with a −1 radial decay. The function describes the 3D cometary plasma distribution surprisingly well. The expression in the first parenthesis determines the latitudinal behaviour of the electron density. Here measures the relative weight of the latitude dependent part, 0 is the latitude where the electron density has a maximum. The expression in the second parenthesis determines the longitudinal behaviour of the density, where gives the relative weight of the longitudinal variations and 0 is the longitude where the electron density has a maximum.
If we carefully examine the density curve shown in the bottom panel of Fig. 2, we find that there is a decreasing trend, the recurring structures have generally diminishing magnitudes. This feature can be easily understood by taking into account the diminishing activity of the comet. The simple Ansatz presented in Eq. (1) cannot capture this feature, and thus the model in its simplest form strongly overestimates the last two structures (19-23 Sep) of the density curve. We overcome this problem by making our k parameter dependent on the distance from the Sun (primarily the Sun-comet distance ( ) determines the cometary activity). According to Hansen et al. (2016), the production rate suffers approximately a tenfold decrease in every 0.58-0.6 AU travelled away from the Sun. In the last month the comet moved from 3.68 to 3.84 AU, which suggests that the activity was almost halved during this period. We can take into account this factor by defining ( ) as where 0 = 3.68 AU, = 0.6 AU and 0 = ( 0 ) is the new constant parameter to fit. The quality of the fit depends only slightly on the value of D, similar results can be achieved if we choose the parameter anywhere from the 0.5 AU < < 1.2 AU range. We fitted the density measurements by inserting the time variation of the ( , , ) coordinates of the spacecraft into the simple function presented in Eq. (1), and used Eq. (2) to take into account the influence of the changing Sun-comet distance ( ) on the value of = ( ). Fig. 4 shows the very good agreement between the model (red curve) and the MIP cometary plasma density in situ measurements (black). After combining Eqs. (1) and (2), the final form of the model can be written as (3) and we used 0 = 3.68 AU and = 0.6 AU as above. We do not expect such a simple model to account for all the short scale features observed in the measurements, which can be associated with local plasma dynamics and/or variations in solar wind forcing. However, the model reflects the large-scale behaviour very well, in particular the main periodicity, the abrupt drops after the main density peaks, and also the presence of secondary peaks next to the main peaks. Moreover, it fits well both the peak widths and amplitudes. The amplitudes and sometimes the positions of the third and fourth peaks show significant deviations, which are probably due to a more complex source structure than the simple first order angle dependence we used. The amplitude of the main peaks are usually somewhat underestimated by this first-order model. This means that the angular distribution features a sharper (higher-order) peak over the highest activity source region.
We assume a single smoothly varying source region in our model, from which the majority of the ionized particles originate. The fact that this simple assumption describes the density distribution so well probably means that most of the small scale density variations are smoothed out before the gas and the plasma reach the sampled altitudes. This does not require a collisional process, since the measured density is the sum of the contributions of all the individual sources. If the measurement is performed far enough from the sources (the distance from the surface is much larger than the source separation) then all the sources are summed up with similar geometric attenuation factors, and the result will be a smooth function reflecting the average source strength. (In contrast, close to the surface, material sources closest to the spacecraft would dominate the measurements, but the 4 km minimum altitude of our orbits ensure an already significant averaging. As the main peaks in the data occur close to the surface, some traces of the low altitude inhomogeneities show up as deviations from the model near the main peaks.) Since our model captures the main features of the plasma density structure the observed deviations can be used to investigate the fine structure caused by transient or local effects. Such effects can be for example the local spatial variations due to the fine structure of the source or the temporal variation of the ionization rates or even solar wind transients. The event on 17 September is the most significant example of such deviations, we are currently investigating its possible cause. The only peculiarity of the 17 September event revealed so far is an excess of suprathermal electrons. All the other plasma density peaks are accompanied by a depletion and cooling of the suprathermal electrons, which can be expected as these events take place in the densest region of the neutral atmosphere and the electrons are cooled by neutral collisions. The excess on 17 September is indicative of a singular process replenishing the suprathermal electron component. The anomalous increase in plasma density observed on Sept 17th is therefore likely caused by an increase of electron impact ionization associated with this excess of suprathermal electrons, as reported for different periods in previous studies Heritier et al., 2018) In agreement with previous results based on measurements from earlier phases of the comet's lifetime Edberg et al., 2015;Vigren et al., 2016;Nemeth, 2020), the electron density falls off with approximately r −1 in the fading coma of comet 67P. This r −1 dependence of the electron density is a remarkably persistent feature of the cometary environment everywhere, where the energy density of the cometary plasma dominates over that of the solar wind. Further away form the nucleus, where the effects of the solar wind dominate, this rule is not expected to hold. Behar et al. (2018) created a semianalytical model of this region, in which the transition from new born ions into pick-up ions is treated as a loss term for the newborn ion population. Nilsson et al. (2018) interpreted the energy spectra of the pick-up ions in terms of their source region ion density, which appeared to fall off as r −2 in accordance with the expected production rate.
It is important to note that the validity of using separate radial and angular variables in our analytical model is equivalent with the angular structure being independent of the radial distance. This means that the plasma motion is essentially radial at the distances considered here, irrespectively of the location in latitude and longitude.
The electron density features a maximum in the southern hemisphere, the best fit to the measured MIP data is achieved when we set the location of the maximum of the density around the south pole. This result agrees well with the findings of investigations of the neutral density after perihelion Kramer et al., 2017;Combi et al., 2020) that found an active southern hemisphere and Z. Nemeth et al. showed the separation of the sub-solar point and the highest density areas above the comet. According to the findings of Combi et al. (2020) during the last months of the Rosetta mission the dominant CO 2 surface activity distribution showed strong latitudinal dependence with a maximum at latitudes around 90 • in the southern hemisphere and a longitudinal dependence with a faint maximum at longitudes around 0 • . Kramer et al. (2017) reported that in May 2016 the neutral densities had a maximum above longitudes around −10 • . In agreement with this we have found an electron density maximum at 0 = −15 • in our model. Values between −30 • and 0 • give similar results.
This study shows that the latitude plays a very important role in the density distribution: the high = 0.76 latitudinal modulation amplitude means that the density over the north pole is only 14% of the density over the south pole, the ratio of the two values is (1 − 0.76)∕(1 + 0.76) ≈ 0.14. In contrast, the longitudinal position influences the density only slightly, with an = 0.13 modulation amplitude. Thus the minimum in longitude is 77% of the maximum since (1 − 0.13)∕(1 + 0.13) ≈ 0.77.
A radial distance -latitude map of the model density distribution is shown in Fig. 5, to be compared to Fig. 3. The model explains the cometary plasma densities measured along the Rosetta orbiter trajectories very well. Fig. 6 is a longitude-latitude map of the electron density 4 km over the centre of the nucleus (top panel). This 4 km altitude is the minimum altitude sampled in this time period, but according to Heritier et al. (2017), this altitude also coincide with the height featuring the peak ionospheric density. The bottom panel projects the density contours onto a map (El-Maarry et al., 2016) showing surface features and regions of 67P. The highest densities were measured over the Bes region, while the lowest activity corresponds to Seth.
These maps show the plasma distribution in comet fixed coordinates. Since at this time of the mission both the neutral flow and the plasma is tenuous, the bulk motion of plasma particles points radially outwards from the cometary nucleus in inertial frame. This means that in comet fixed coordinates they move along slightly bent trajectories. Since close to the nucleus the radial flow speed is much larger (∼500-1000 m/s, Hansen et al., 2016) than the apparent tangential speed (∼2 m/s at 15 km from the comet) in the comet fixed frame, this effect does not change the picture described above; close to the nucleus the plasma motion can be assumed to be approximately radial in comet fixed frame as well. In the 4-15 km radial range of our study we see a plasma cloud radially expanding with respect to the comet and preserving the original latitude-longitude distribution of the source surface.

Conclusions
Near the end of the Rosetta orbiter operations, although comet 67P was more than 3.7 AU from the Sun, in situ measurements still show clear signs of a weak but well defined cometary plasma environment. During the last month of the Rosetta operations, in September 2016, the spacecraft moved along a periodic, recurrent orbit that made it possible to study the 3D spatial distribution of the plasma density near the nucleus. In this paper, we derived a simple and useful model to explain the plasma density distribution in the coma of comet 67P in September 2016.
Based on in situ MIP electron density measurements we defined a simple distance, latitude and longitude dependent first order cosine function to model the 3D spatial distribution of the cometary plasma. The model features a −1 dependence on the distance from the centre of the nucleus. It slightly depends on the Sun comet distance as well, because the cometary activity diminishes as the comet moves away from the Sun. A remarkable advantage of this model is that the four variables of interest are separated, thus showing the role of each independent variable in the 3D mapping of the cometary ionosphere.
This 3D cometary plasma density distribution model reproduced the Rosetta MIP observations remarkably well. The model reflects the observed structures in the plasma density distribution, in particular the main periodicity, the abrupt drops after the main peaks, even the presence of secondary peaks next to the main peaks; it fits well the peak widths as well as the amplitudes. We trust that this first-order 3D model of the cometary ionosphere of 67P will also make it possible to better understand the local plasma dynamics identified as local discrepancies between the Rosetta plasma observation and the model described in this work.
The plasma density distribution shows a strong latitudinal dependence: the plasma density is highest above the southern hemisphere. This is consistent with the neutral density observations after the comet's perihelion passage Kramer et al., 2017;Combi et al., 2020). Indeed, the southern, nightside hemisphere produces more plasma than the sunlit northern hemisphere -mostly due to the higher neutral outgassing rates. Our model shows that the plasma density can be described well by assuming only a single plasma source in longitudes around −15 • . This also correlates with the findings of Kramer et al. (2017) who found that in May 2016 the neutral densities had a maximum above longitudes around −10 • .