RPC-MIP observations at comet 67P/Churyumov-Gerasimenko explained by a model including a sheath and two populations of electrons

The response of the mutual impedance probe RPC-MIP on board Rosetta orbiter electrostatically modeled considering an unmagnetized and collisionless plasma with two Maxwellian electron populations. A vacuum sheath surrounding the probe was considered in our model in order to take the ion sheath into account that is located around the probe, which is immersed in the cometary plasma. For the first time, the simulated results are consistent with the data collected around comet 67P/Churyumov-Gerasimenko (67P), but strong discrepancies were identified with the previous simulations that neglected the plasma sheath around the probe. We studied the influence of the sheath thickness and of the electron populations. This work helps to better understand the initially unexpected responses of the mutual impedance probe that were acquired during the Rosetta mission. It suggests that two electron populations exist in the cometary plasma of 67P.


Introduction
Originally intended for geological and archaeological prospection, mutual impedance probes have been adapted (Storey et al. 1969) and successfully used in space plasma investigations since the 1970s (Beghin & Debrie 1972;Décréau et al. 1978;Beghin et al. 1982;Bahnsen et al. 1986;Trotignon et al. 2007), and it is planned to use them for future space missions such as BepiColombo to the planet Mercury (Trotignon et al. 2006) or the Jupiter icy moons explorer (JUICE). This active electrostatic probe introduces a sine wave electrical current with a constant amplitude in the probe surroundings, and it measures the induced electrical potential drop between two receiving electrodes (Storey 1998). The potential drop is proportional to the inverse of the relative permittivity in a conventional dielectric medium that shows a scalar permittivity. The situation can be much more complex in a plasma because the apparent relative permittivity is no longer a scalar, but is a function of the frequency, the distance relative to the emitting electrodes, the plasma electron and ion velocity distribution functions, the collision rates, and the magnetic field.
In a plasma that surrounds a comet, the electron collision frequencies as well as the electron cyclotron frequency are negligible compared to the plasma frequency. Thus the cometary plasma in the vicinity of the plasma frequency can be considered collisionless and unmagnetized. Based on the same approach as in Geiswiller et al. (2001) and Béghin et al. (2005), electrostatic simulations have been performed about 20 yr ago in order to predict the response of the mutual impedance probe in such a plasma (Geiswiller 2001). The probe was modeled as a mesh immersed in a Maxwellian plasma, and the Debye length spanned from several millimeters up to a few centimeters. These conditions correspond to what was expected by way of a surrounding around comet 46P/Wirtanen, the original destination of the European Space Agency's Rosetta spacecraft. However, a large discrepancy was observed between simulations and experimental data after the first measurements had been acquired around comet 67P/Churuymov-Gerasimenko (hereafter 67P). In particular, the Debye length was observed to be much larger, from a few tens of centimeters to a few meters. Thus, the first modification of the simulation was to extend the Debye length toward higher values. Then, a second Maxwellian electron population was introduced in the plasma because suprathermal electrons have been commonly observed in space plasmas, superposed on a Maxwellian-like core (Pilipp et al. 1987;Pierrard & Lazar 2010). These modifications failed to bring the simulated results closer to experimental data (Béghin et al. 2005;Gilet et al. 2017). Finally, it was suspected that the problem came from the ion sheath between the probe and the plasma, which is caused by the mobility difference between electrons and positive ions in the cometary plasma; this had been neglected in the modeling. This omission was probably reasonable in the terrestrial ionosphere where the electrostatic sheath that surrounds the mutual impedance probes used in the different experiments launched in the 1970s was far smaller than the dimensions of the probes. However, the expectations of the Debye length at comet 67P show that the Debye length should be comparable to the dimensions of the mutual impedance probe. Because the sheath thickness is usually in the same order of magnitude as the Debye length, a study of the electrostatic influence of a large sheath Table 1. rms value and phase shift of V A = |V A | × exp ( j (ωt + φ A )) and V B = |V B | × exp ( j (ωt + φ B )) for the different operating modes of the probe. Mode 1V π 0 Notes. The rms values are given with respect to the spacecraft potential, and they correspond to the so-called half-level of emission of the probe (see text for details). n.a. = not applicable. ω = 2π f . around the probe was necessary. Consequently, the main aim of this work is to introduce a large sheath to modeling the response of the mutual impedance probe on board the Rosetta spacecraft. The probe is presented in Sect. 2, and the modeling including a vacuum sheath around the probe is fully described in Sect. 3. The influences of the probe surroundings and of the input parameters of the model are studied in Sect. 4. A comparison between simulated results and experimental data is performed in Sect. 5, followed by an estimation of the electron velocity distribution function (evdf) at comet 67P in Sect. 6. Finally, we discuss the different hypotheses in our modeling in Sect. 7.

Mutual impedance probe
The mutual impedance probe of the Rosetta mission has been accurately described by Trotignon et al. (2007). We hereafter concentrate on the so-called SDL operational mode of the RPC-MIP instrument, which consists of four cylindrical conducting rods (T 1 , T 2 , R 1 , R 2 ) supported by a conducting carbonreinforced plastic frame (see Fig. 1). The active rods T 1 and T 2 (transmitters) are connected through the 0.12 pF series capacitor C 1 and C 2 to an individual voltage sine generator operating in the range [27 kHz, 3.5 MHz] (see Fig. 2). The rms value of the signal delivered by the generators (|V A | for T 1 and |V B | for T 2 ) can be set to 0V, 0.25V, 0.5V, 1V, or 2V relative to the spacecraft potential V SC . The rms value and the phase shift φ A,B of these voltages versus the different operating SDL modes of the probe are detailed in Table 1. For a large part of the RPC-MIP operations, the voltage levels were 1V (half-level) from about the beginning of the cometary operations until September 2015, and then 2V (full level) until the end of the cometary operations in September 2016. The passive rods R 1 and R 2 (receivers) were connected to the spacecraft through the 0.12 pF series capacitor C 3 and C 4 (see Fig. 2). The frame of the probe as well as the boom supporting the probe were at the spacecraft potential V SC . Because of the voltage drop induced by the series capacitors C 1 , C 2 , C 3 , and C 4 , the potentials V T 1 and V T 2 of the transmitters and V R 1 and V R 2 of the receivers are at the angular frequency ω with the impedances where q T 1,2 is the overall charge carried by T 1 or T 2 , and q R 1,2 is the overall charge carried by R 1 or R 2 . The value ∆V R of the induced voltage between the receivers of the probe is given by After it was acquired, |∆V R | was normalized by 0.05 mV, and the spectral power of the signal P = 10 log 10 (20 |∆V R |) 2 is referred to as the response of the probe (in dB units). Thus, 0 dB corresponds to |∆V R | = 0.05 mV. It has been experimentally checked with the engineering model of the probe that a 71 mV rms sine voltage operating at 28 kHz that is directly applied between the receivers provides a power response of 63 dB. The response of the probe is increasingly attenuated with the frequency by a low-pass anti-aliasing filter, and the attenuation at a frequency of 1 MHz was experimentally estimated between 2 and 3 dB with respect to the level provided at 28 kHz. The geometry and the electrical bias of the considered parts of the probe are listed in Table 2.

Electrostatic modeling
The discrete surface charge distribution (DSCD) approach (Béghin & Kolesnikova 1998) is considered in this work. It consists of distributing point charges on the surface of the boundaries in order to determine the potential or the electric field within the electrostatic approximation, that is, if the size of the considered domain is much smaller than the electromagnetic wavelength associated with the frequencies under consideration.
The electrostatic modeling depending on the probe surroundings is detailed for three different situations. The first is the probe in vacuum (Sect. 3.1), the second is the probe immersed in a plasma without an isolating sheath around the probe (Sect. 3.2), and the third situation is the probe isolated from the plasma by a vacuum sheath (Sect. 3.3). The probe, the boom, and the sheath (if any) are modeled as grids, and their respective discrete charges are located on the nodes of the grids in question. The dielectric parts isolating the transmitters or receivers from their holders have not been considered. When the charge distribution at the probe surface was determined, the overall charge carried out by each receiver was computed. These charges lead to the potentials of the receivers through Eq. (2). Then the spectral power of the voltage difference ∆V R between the receivers (or the response of the probe) is finally computed through P = 10 log 10 (20 |∆V R |) 2 .

Probe in vacuum
According to the previous section, the potential of every element i of the probe or of the boom reads where Z S is the impedance crossed by the charges q S of a surface S (q i belonging to q S ) to reach the spacecraft, and V 0i is the potential applied by a generator (if any) on S . Consequently, for the subscripts i belonging to T 1,2 , we have V 0i = V A,B , Z S = Z C 1,2 . For the subscripts i corresponding to R 1,2 , we have V 0i = 0, Z S = Z C 3,4 . For the subscripts i corresponding to the rest of the probe and for the boom, we have V 0i = 0, Z S = 0.    In the following, the charges located on the probe and on its boom are referred to as q i or q j with i, j ∈ [1; N] (N charges) and j i. Assuming the superposition principle, the potential V i induced on an elementary surface S i of the probe or of the boom can also be expressed according to the charge distribution at the probe surface, where r i j is the distance between the ith point charge q i located at the center of S i and the N−1 other charges q j at the surface of the probe and of the boom. q i /4πε 0 α i is the potential induced on S i by its own charge q i , which is supposed to be uniformly distributed on this surface, where d is the distance from the center of S i . The details computing α i for the facets of cylindrical or spherical objects are included in Appendix A. According to Eqs. (6) and (7), the potential distribution at the surface of the probe leads to N equations with N + 1 unknowns (q i and V SC ), This set of equations is completed by a neutrality condition at the considered frequency over the conducting surfaces, Thus the determination of the charge distribution over the probe surface in vacuum requires solving a linear system of N + 1 equations with real coefficients.

Probe in a plasma without a sheath
In order to take the dielectric characteristics of the plasma into account, Eq. (9) has to be modified when the probe is immersed in an infinite homogeneous kinetic plasma without any sheath between the probe and the plasma, Here, φ/φ 0 is the radiated potential of a pulsating point charge in the plasma normalized by the potential radiated by this charge in vacuum. It can be seen as the inverse of the plasma relative permittivity function, which depends on the frequency, the distance r i j from the corresponding point charge j, the electron and ion velocity distribution functions (evdf and ivdf), the electron/neutral collisional rate, and the magnetic field. This function was computed for a Maxwellian evdf in collisionless, non-magnetized plasma and for a frequency significantly higher than the ion plasma frequency (Beghin 1995). More recently, φ/φ 0 was determined in the same conditions for a twoelectron temperature plasma modeled by a double-Maxwellian evdf (Gilet et al. 2017) and for a kappa evdf (Gilet et al. 2019). β i is analogous to α i , but for an object in a homogeneous infinite kinetic plasma. Its expressions for the facets of cylindrical or spherical objects are included in Appendix B.
As in vacuum, the charge neutrality condition on conducting surfaces (Eq. (10)) applies. Therefore, determining the potential distribution over the probe surface in a plasma without a sheath requires solving a linear system of N + 1 equations with complex coefficients and requires knowing the function φ/φ 0 . This model has been considered in the electrostatic simulations of the response of the mutual impedance probes embedded in the terrestrial ionosphere or during space plasma investigations before this work (Storey et al. 1969;Grard 1969Grard , 1997Beghin & Debrie 1972;Chasseriaux et al. 1972;Rooy et al. 1972;Pottelette et al. 1975;Décréau et al. 1978;Pottelette & Storey 1981;Beghin et al. 1982;Storey 1998;Geiswiller et al. 2001;Béghin et al. 2005;Trotignon et al. 2007;Geiswiller 2001), implicitly meaning that the probe and the boom have an average electrical potential that corresponds to the plasma potential or that the electron photoemission at the probe surface is very high. Neglecting the ion sheath surrounding the probe is sometimes applicable if the sheath thickness is much smaller than the distance between the transmitters T 1 and T 2 of the probe (Geiswiller 2001). This can be argued to be the case for the Rosetta mission based on the expected Debye length around comet 67P.

Probe isolated from the plasma by a vacuum sheath
The occurrence of the sheath is due to the mobility difference between the electrons and the positive ions in the plasma. An electron-depleted region builds up around any object that is exposed to a plasma associated with a negative electric field that balances the electron and positive ion fluxes at the surface of the object when this last is electrically floating. If the object is biased to a voltage value lower than the plasma potential (which is the case for the RPC-MIP instrument), the (positive ion) sheath remains electron-depleted. In a double-Maxwellian evdf plasma, the sheath is governed by the hot electrons, and the sheath thickness is expected to be comparable to the Debye length of this hot electron population. This positive ion sheath is considered as a vacuum region in our modeling, which is justified because the ions are unable to follow the oscillations of the electric field at frequencies that are close to the plasma frequency. Thus, the surroundings of the probe become inhomogeneous because of the interface between plasma and sheath.
Originally developed for homogeneous environment, the DSCD method applies to inhomogeneous dielectrics by introducing fictitious charges on both sides of the interfaces. As an extension of the method of the electrostatic images, these fictitious charges are intended to provide the contribution of the (real) charges embedded in a dielectric homogeneous medium in a second homogeneous dielectric medium interfacing with the first, and vice versa. Béghin & Kolesnikova (1998) have demonstrated that a kinetic plasma is equivalent to the unusual situation of a dielectric having a free charge density implanted in its entire volume (Béghin & Kolesnikova 1998). Thus, the occurrence of a vacuum sheath around the probe is modeled by introducing fictitious charges on both sides of the plasma-sheath interface. The fictitious charges on the internal side of the plasma-sheath interface are called q 1k or q 1m . They correspond to the plasma contribution on the electrostatic potential in the vacuum sheath. The fictitious charges on the external side of the plasma-sheath interface are called q 2k or q 2m . They correspond to the contribution of the (real) charges located inside the sheath (with the exception of q 1 charges) on the electrostatic potential induced in the plasma region. On both sides of the interface, we have k, m ∈ [1; N 1 ] (N 1 charges on both sides of the interface) and m k. The distance between the internal side holding q 1k and the external side holding q 2k of the plasma-sheath interface is referred as 2 × dr. The sheath region is assumed to be electronfree, and the plasma region is assumed to be homogeneous. Thus the potential inside the sheath is computed by considering q i and q 1k , whereas the potential in the plasma region is computed according to q 2k . Consequently, Eq. (9) has to be modified in order to consider the influence of the fictitious charges q 1k at the sheath edge on the distribution of the potential at the surface of the probe, with r ik the distance between the ith point charge q i and the kth point charge q 1k located at the plasma-sheath interface.
The continuity of the potential at the plasma-sheath interface reads (see Appendix C for the details of the computation) The mutual impedance probe is located inside a cylindrical electron-free sheath ending by a hemisphere. The spacecraft body is modeled by the square plate. Plasma is only present at the external side of the sheath, which is presented in sectional view in order to show the probe.
r ki is the distance between the kth point charge q 1k and the ith point charge q i . r km is the distance between the kth point charge q 1k and the N 1 − 1 other charges q 1m of the internal side of the plasma-sheath interface. The Maxwell-Gauss equation or the conservation of the perpendicular component of the displacement at the plasmasheath interface reads (see Appendix C for the details of the computation) .
with r ki1 and r ki2 the distances between the kth point charge q 1k (r ki1 ) or q 2k (r ki2 ) with the ith point charge q i . r km1 and r km2 are the distances between the kth point charge q 1k or q 2k with the mth point charge on the inner (r km1 ) or the outer (k km2 ) interface. Finally, Eq. (10), which describes the charge neutrality over the conducting surfaces, is included to build a linear system of N + 2N 1 + 1 equations with complex coefficients that allow us to determine the spacecraft potential and the charge distributions on the probe, the boom, and the internal and external sides of the plasma-sheath interface according to Eqs. (10), (12), (13), and (14). The N + 2N 1 + 1 complex unknowns are q i , q 1k , q 2k , and V SC , respectively, this last being the spacecraft potential at the angular frequency ω = 2π f .
The geometry of the sheath and of the spacecraft were simplified in our simulation. The mutual impedance probe is located inside a cylindrical sheath with a radius R and a constant length L of 3 m ending by a hemisphere having the same radius (see Fig. 3). In order to compute the perpendicular component of the electric field at the interface (Eq. (14)), the q 1k charges are located on the internal side of the sheath with a radius R − dr, while the q 2k charges are located on the external side of the sheath with a radius R + dr with dr = λ Dh /10 4 . We note that λ Dh corresponds to the Debye length of the hot electron population in the case of a double-Maxwellian evdf. In oder to mimic the spacecraft body, a square conductive plate was also introduced, as displayed in Fig. 3. This sheath geometry is obviously a crude modeling of the real sheath around the probe. This last was investigated using the SPIS software (Johansson et al. 2016).

Simulation setup
The response of the probe was modeled in vacuum and in a vacuum sheath surrounded by a kinetic homogeneous infinite plasma (see Fig. 3). Single-and double-Maxwellian evdf were considered, but we only report double-Maxwellian evdf simulations because they agree much better with experimental data than the single-Maxwellian evdf simulations. The double-Maxwellian plasma is composed of a hot and a cold Maxwellian electron population, whose densities are n h and n c and whose temperatures are T h and T c , respectively. λ Dh stands for the Debye length of the hot electron population, and n tot = n c + n h is the overall electron density. Ω = f / f p is the normalized frequency, where f is the frequency of the excitation and f p is the plasma frequency. The considered electron density in f p is n tot where q e = −1.6 × 10 −19 C is the charge of an electron, m e = 9.1 × 10 −31 kg is the mass of an electron, k B = 1.38 × 10 −31 J K −1 is the Boltzmann constant, and ε 0 = 8.85 × 10 −12 F m −1 the vacuum permittivity.

Simulation in vacuum
The N point charges q i = |q i | × exp ( j (ωt + ϕ i )) and the oscillating spacecraft potential V SC were determined by solving the linear system of Sect. 3.1 considering N = 3348 (2360 discrete charges for the probe, 888 for the boom, and 100 for the square plate). Similar results were obtained for N = 928 (600 discrete charges for the probe, 228 for the boom, and 100 for the square plate). Except for the so-called antiphase mode, where T 1 and T 2 operate in counter-phase (push-pull) and for which the charges are mainly distributed on those transmitters, |∆V R | is highly dependent on the position of the probe with regard to the end of the boom. Moreover, the presence of a Langmuir probe at the end of the boom might influence the distribution of the charges because it is grounded when the impedance probe operates. Therefore we set the end of the impedance probe 15 cm before the end of the boom in the simulations in order to match the experimental responses in vacuum. Table 3 summarizes the simulated response of the probe in vacuum for the different operating modes with regard to experimental responses acquired at about 3.4 AU from the Sun in September 2014, which are expected to be comparable to the vacuum signals. Finally, the simulated distribution of q i at the surface of the probe in vacuum is shown in Fig. 4 for the antiphase mode.

Sensitivity of the model with a vacuum sheath
The simulation of the probe surrounded by a vacuum sheath in a kinetic plasma was carried out with 600 discrete charges on the probe, 228 on the boom, and 100 on the square plate (N = 928). The step between two point charges located on the same side of the plasma-sheath interface was set between λ Dh /2 and λ Dh /4, leading, for example, to N 1 = 410 charges q 1k and q 2k on both sides of the interface for λ Dh = 50 cm, R = 2λ Dh and a step of λ Dh /2. It was checked for every run that increasing the number of point charges did not significantly change the simulated results.
Most of the experimental data were acquired in phase mode because the highest signal is observed experimentally at around the plasma frequency. We therefore now focus on this operational mode.
According to the simulations, the spectral power is particularly sensitive to the sheath radius. This is illustrated in Fig. 5 (top left) for the so-called phase mode where the spectral power tends toward its value in vacuum when the sheath radius R is increased, as expected. The influence of the Debye length on the hot electrons, the ratio n h /n tot , and T h /T c are presented in Fig. 5, where we observe that the spectral power converges at high frequencies toward almost the same value regardless of these parameters. It is noteworthy that the sheath radius acts as a baseline on the response of the probe, whereas λ Dh , n h /n tot , and T h /T c affect the shape of the response of the probe in the vicinity of the plasma frequency. Figure 6 confirms that the sheath radius R is the most influential parameter on the converging level of the response of the probe above the plasma frequency. This finding was verified for λ Dh and R ∈ [0.4; 2] m, T h /T c ∈ [5; 250], and n h /n tot ∈ [0.1; 0.9] and remains the same regardless of the operating mode of the probe. The modification of the evdf does not affect the response level of the probe for Ω 1 because of the function φ/φ 0 , which is always found in a range of values close to the inverse of the cold plasma permittivity ε c = 1 − Ω −2 when Ω 1, regardless of n h /n tot , T h /T c , and the ratio between the distance and λ Dh . This is illustrated in Fig. 7, where |φ/φ 0 | covers a wide range of values for Ω = 0.6, whereas this range is significantly narrower and close to 1/ε c for Ω = 2.
Moreover, as for the model in which the sheath located around the probe was neglected, a large resonance peak (maximum level of the response) is observed at a frequency f 0 , which can be up to 20% lower than the plasma frequency f p (or Ω = 1) in the case of a double-Maxwellian evdf plasma. As for the experimental RPC-MIP responses, there is a large antiresonance peak (minimum level of the response) located at a frequency f 1 that is lower than f 0 when the sheath is considered in the model. This peak vanishes when the sheath is neglected in the model (see Fig. 8).  Finally, similar simulated responses corresponding to different input parameter values (R, λ Dh , n h /n tot , T h /T c ) have not been observed in the frequency domain f p /2 ; 2 f p .

Simulation results
In order to determine whether our modeling is valid, the simulated results were compared with experimental data that were acquired during the Rosetta mission by the mutual impedance probe MIP of the Rosetta Plasma Consortium (Carr et al. 2007) on board the spacecraft. Figure 9 (top panel) displays a set of spectral responses collected on 20 March 2016 between 4:10 and 5:10 (UTC), when a sequence of operations that successively involved the four different SDL modes of RPC-MIP had been played on Rosetta. The reproducibility of the experimental responses allows us to consider that the plasma was stable. This enables us to efficiently compare the different RPC-MIP SDL operational modes by considering the same plasma conditions during the whole interval. The accordance with the simulated spectral responses can be evaluated in the bottom panel of the same figure, where the four operating modes were simulated in the same plasma conditions. It is noteworthy that the simulation is always consistent with the experimental data regardless of the operating mode of the probe, as soon as the plasma sheath is taken into account in the modeling (plain lines). In contrast, this is clearly not the case when the sheath is not taken into account in the model (dashed lines). Simulations predict the plasma signature to be more marked when the probe operates in phase mode and to be particularly singular in the singleE1 operating mode, as observed during RPC-MIP operations. After matching experimental and simulated probe responses many times, we conclude that the best accordance was always reached when the Debye length of the hot electron population λ Dh was similar to the sheath thickness R, which was chosen according to Fig. 6. Then the frequency range located between 0.5 × f p and f p can be more or less well fit by adjusting the density and temperature values of the two electron populations, with most of the time n h /n tot = 0.3 ± 0.1 and T h /T c = 30 ± 20. Some experimental   responses can be very well fit by simulated responses (Fig. 10, 11, 13), but a qualitative agreement can only be reached over this frequency range for a significant amount of experimental data (Fig. 12). This is probably due to the hypotheses considered in our modeling, which are discussed in Sect. 7. The results are gathered in Fig. 14 and show the electron energy distribution functions (eedf) that provide the best agreement between simulations and experimental responses acquired by the probe; they are presented in Figs. 10-13. The eedf is the sum of two Maxwellian eedfs and reads where the energy E and the temperatures T c,h are in eV.

Estimating the evdf at comet 67P
The RPC-MIP acquisitions in SDL mode at comet 67P are consistent throughout the entire frequency range of RPC-MIP with the electrostatic DSCD simulations based on a model including a large ion sheath around the probe. The agreement was much better with a double-Maxwellian evdf plasma than with a single-Maxwellian evdf plasma, suggesting that the evdf at comet 67P presented an evdf that was broader than a single-Maxwellian evdf. Keeping this consideration in mind, we can estimate the plasma evdf during the Rosetta mission at comet 67P. Using the model described above, it is possible to compute the response of the probe for a variety of values of n h /n tot and T h /T c , λ Dh and R. However, fitting the huge amount of mutual impedance spectra acquired during the comet phase of the Rosetta mission is beyond the scope of this work. Thus we chose another way to process the experimental data. The simulation results and the experimental data often have only a qualitative agreement around the plasma frequency, where the mutual impedance responses are more sensitive to the details of the electron distribution function. Instead, the agreement is most of the time very good above the   plasma frequency, leading to the determination of the sheath radius R (see Fig. 6). After we realized that the best accordance between the simulated and the acquired data is reached for a sheath radius similar to the Debye length of the hot electron population (λ Dh ≈ R), when the ratios n h /n tot and T h /T c were in the range 0.3 ± 0.1 and 30 ± 20 (see Sect. 5), respectively, the experimental data could be processed rapidly according to Eqs. (15) and (16).
According to the multiple comparisons between simulated and experimental probe responses (Sect. 5), the following ranges were considered in order to estimate the plasma parameters at comet 67P: These considerations lead to the relative uncertainties we list in Table 4. We estimated the plasma evdf for the measurements performed on 20 March 2016 and on 24 August 2016 by the mutual impedance probe of Rosetta. The comet was located at 2.62 AU from the Sun, and the spacecraft altitude was about 12 km from the comet center on 20 March 2016. These distances became 3.63 AU and varied from 5 to 14 km during the observations reported on 24 August 2016. Figure 15    It is noteworthy that a hot and a cold electron population emerge for these days, but this finding only points out an excess of hot electrons in the evdf compared to a single-Maxwellian distribution function. At this stage, it is not possible to assert from the RPC-MIP diagnostics that the evdf is made of two different electron populations, that is, which are created by different mechanisms. However, the occurrence of a hot and cold electron population has been identified at comet 67P (Eriksson et al. 2017;Gilet et al. 2017;Engelhardt, I. A. D. et al. 2018). A hot electron population is observed in the range [5-15 eV], while a cold electron population is observed in the range [0.1-1 eV]. Interestingly, the temperatures of the two electron populations are consistent with those observed independently during another time interval by the Langmuir Probe RPC-LAP (Eriksson et al. 2017), which is another instrument on board the Rosetta orbiter. This further supports the model described in this work. The hot electron population (also referred to as a "warm" electron population in the literature to avoid confusion with another suprathermal electron population that was observed in the cometary plasma and is referred to as a different "hot" component in other studies) is likely to be associated with the cometary electrons that are freshly ionized in the cometary ionosphere, either by photoionization (Vigren et al. 2016), electron impact ionization, or a mix of the two (Galand et al. 2016;Heritier, K. L. et al. 2018). Instead, the observed cold electron population is most probably associated with electrons that would have cooled down because of collisions on cometary neutrals (Eriksson et al. 2017). It is interesting and somewhat surprising to observe signatures of Table 4. Relative uncertainties on the determination of the evdf parameters.

Relative uncertainty
cold electrons in the data even when 67P was far from the Sun (3.6 AU) and the neutral outgassing rate was consequently low. During this period, the so-called electron exobase (Mandt et al. 2016), which represents the region around comet 67P were the electron dynamics is expected to be dominated by collisions on neutrals, was not expected to be located above the comet surface, so that one could have expected to observe no cold electrons. The evdf might have been different from a pure Maxwellian evdf on this day, and the consideration of a second electron population in the model helped fit the experimental data better. Finally, it may be possible that the approximation of the sheath geometry and the absence of electrons inside it influences the simulation results. This point is discussed in Sect. 7.

Discussion
Considering a vacuum sheath around the mutual impedance probe RPC-MIP exposed to the cometary plasma enabled us to make the simulations consistent with the experimental data that were collected during the Rosetta mission. This highlights that the sheath surrounding the mutual impedance probe needs to be considered for electrostatic modeling. As observed during the Rosetta mission, the plasma signature on the simulated response of the probe is much more marked in phase mode than in the other operating modes (see Fig. 9). The singular behavior of the response in singleE1 mode compared to the behavior of the response in the other operating modes shown in Fig. 9 is also observed experimentally. This mode was suspected to operate incorrectly before this work because of this singularity. The discrepancy between the response in singleE1 and singleE2 modes is due to the asymmetry that is induced by the boom that supports the probe, by the spacecraft body, and by the sheath end.
The remaining discrepancies between the simulations and the experimental data are attributed to the assumptions made by the model, such as the shape of the sheath and the absence of electrons inside it. A vacuum sheath around the mutual impedance probe leading to a sharp transition of the electron density at the interface with the plasma remains a first approximation of the real sheath, where the electron density is expected to increase with the distance from the probe. Considering an electron density profile in the sheath should be feasible with the DSCD approach by considering multi-layer plasma regions stacked around the probe with different evdf and fictitious charge distributions at the interfaces. Despite a strong increase in random access memory needs and computation time, such a modeling would require to know the spatial distribution of the evdf in the sheath. This is expected to be reachable with a future release of the modeling software SPIS. However, if the evdf in the sheath does not correspond to a Maxwellian or a to double-Maxwellian evdf, which is likely to append, the function φ/φ 0 will have to be determined as it is an input of the simulation. The evdf inside the sheath is expected to be a shifted and/or truncated (double-) Maxwellian because of the local electric field. If confirmed, the incomplete plasma dispersion function (Baalrud 2013) will arise when the function φ/φ 0 in the sheath is computed.
In addition, it is possible that the double-Maxwellian evdf is sometimes a crude approximation of the real plasma evdf at the comet. This could be the reason why a significant amount of experimental data are only in a partial agreement with the simulations, such as in Fig. 12. Thus the simulation with a sheath should also be carried out with another evdf, such as the Kappa evdf, which was observed to correspond to the evdf in the solar wind (Maksimovic et al. 1997) and in the cometary plasma (Clark et al. 2015;Broiles et al. 2016). A mix of a Maxwellian evdf for the cold electrons and a Kappa evdf for the hot electrons might likewise be imagined (Gilet et al. 2019).
Finally, our work points out an interesting and important result regarding the ability of the RPC-MIP experiment to efficiently operate plasmas characterized by a large Debye length. Based on the results of numerical simulations that considered a single-Maxwellian electron population and that neglected the plasma sheath around the probe, the maximum Debye length value that could be reached by the RPC-MIP probe operating in SDL mode was strongly underestimated (20 cm) before our modeling (Trotignon et al. 2007). However, neglecting the sheath is appropriate if the sheath thickness is negligible compared to the size of the probe. A 20 cm long Debye length associated with a sheath thickness of only few centimeters is unlikely to occur, however, because the Debye length is usually comparable to the sheath thickness unless the probe is biased close to the plasma potential or above it, or if the electron emission at the probe is high. These two situations are very unlikely in the case of the Rosetta mission. According to our investigations, which suggest (as expected) that λ Dh ≈ R, cometary plasmas with a Debye length up to 1.5 m (note that this is the Debye length of the hot electron population, the cold electron population is much smaller) would be analyzed by the probe in SDL mode, as shown in Fig. 6, where the response of the probe differs from the vacuum response when λ Dh < 1.5 m. Then the response of the probe for higher values of λ Dh are indistinguishable from the response of the probe in vacuum.

Conclusion
With the insight brought by our electrostatic modeling of the mutual impedance probe RPC-MIP on board the Rosetta orbiter, the initially unexpected experimental responses observed around comet 67P can henceforth be related to the plasma evdf and to the sheath thickness surrounding the probe. We confirmed the occurrence of a large ion sheath around the probe in almost all the Rosetta mutual impedance probe responses that show a plasma signature. The evdf in the cometary plasma is found to present an excess of electrons at high energy levels compared to a single-Maxwellian evdf, which means that a second electron population needs to be introduced in the modeling of the response of the probe to match the experimental responses acquired at comet 67P. It is important to underline that the RPC-MIP experiment is unable to determine whether the cometary plasma at comet 67P was effectively made of two Maxwellian electron populations or if it consisted of a single population with more electrons at high energies than for a single-Maxwellian population. However, the electron populations we discussed here are very likely independent of the spacecraft itself (i.e., they not associated with the spacecraft photoelectrons, whose density is far lower than that of the cometary electrons). Instead, the presence of the two electron populations would correspond to (i) cometary electrons with an energy of about 5-15 eV coming from the ionization of the expanding cometary atmosphere through photoionization and electron impact ionization, mixed with solar wind electrons, and (ii) cometary electrons that have cooled down to lower energies (0.1-1 eV) through collisions with the cometary neutral molecules (mostly H 2 O and CO 2 ).
Even if the modeling that includes the sheath is unable to always exactly reproduce experimental data, preventing an accurate plasma diagnostic to be performed, with the exception of the determination of the plasma frequency (and of the overall electron density through it), which always arises close to the maximum level of the response of the probe according to the simulations. However, after our study, several hypotheses on the density and temperature ratio between the cold and the hot electron populations can be made to estimate the plasma evdf at comet 67P despite the challenging operating conditions of the probe, which has never been intended to operate in such a large sheath.
Finally, an important result of this work is that although the RPC-MIP instrument is completely embedded in the plasma sheath that surrounds the spacecraft and its booms, the resulting observations still enable us to retrieve the plasma parameter away from that same sheath. This points out the non-locality of plasma measurements made by the mutual impedance method, in contrast to the locality of plasma diagnostics performed by Langmuir probes. This emphasizes the high complementarity of these two measurement techniques.   In order to verify the validity of our determination of α i as well as the DSCD method in vacuum itself, simulations were carried out for classical electrostatic configurations for which analytic solutions are available using the method of electrostatic images, for instance. An example is given in Fig A.2 with a comparison between the DSCD result and the analytical result in Fig. A.3.