© G. Wattieaux et al. 2019

1 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 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.

2 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₁, T₂, R₁, R₂) supported by a conducting carbon-reinforced plastic frame (see Fig. 1). The active rods T₁ and T₂ (transmitters) are connected through the 0.12 pF series capacitor C₁ and C₂ 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 ( for T₁ and for T₂) 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₁ and R₂ (receivers) were connected to the spacecraft through the 0.12 pF series capacitor C₃ and C₄ (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₁, C₂, C₃, and C₄, the potentials and of the transmitters and and of the receivers are at the angular frequency ω

with the impedances

where is the overall charge carried by T₁ or T₂, and is the overall charge carried by R₁ or R₂.

The value ΔV_R of the induced voltage between the receivers of the probe is given by (5)

After it was acquired, was normalized by 0.05 mV, and the spectral power of the signal is referred to as the response of the probe (in dB units). Thus, 0 dB corresponds to 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.

Fig. 1 Schematic of the mutual impedance probe mounted on its boom.

Fig. 2 Equivalent electrical arrangement of the emitters and receivers. This schematic does not correspond to the geometry considered in the electrostatic simulation (see Fig. 3).

3 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 .

3.1 Probe in vacuum

According to the previous section, the potential of every element i of the probe or of the boom reads (6)

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 toT_1,2, we have V_0i = V_A,B, . For the subscripts i corresponding to R_1,2, we have V_0i = 0, . 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 probeor of the boom can also be expressed according to the charge distribution at the probe surface, (7)

where r_ij 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πε₀α_i is the potential induced on S_i by its own charge q_i, which is supposed to be uniformly distributed on this surface, (8)

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), (9)

This set of equations is completed by a neutrality condition at the considered frequency over the conducting surfaces, (10)

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.

3.2 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, (11)

Here, ϕ∕ϕ₀ 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_ij 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, ϕ∕ϕ₀ was determined in the same conditions for a two-electron 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 ϕ∕ϕ₀.

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 1969, 1997; Beghin & 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₁ and T₂ 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.

3.3 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 bythe 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₁ charges) on the electrostatic potential induced in the plasma region. On both sides of the interface, we have k, m ∈ [1;N₁] (N₁ 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 electron-free, 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, (12)

with r_ik the distance between the ith point chargeq_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) (13)

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 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 plasma-sheath interface reads (see Appendix C for the details of the computation) (14)

with r_ki1 and r_ki2 the distances between the kth point chargeq_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 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 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⁴. 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).

Fig. 3 Spatial configuration considered in the electrostatic simulations. 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 presentat the external side of the sheath, which is presented in sectional view in order to show the probe.

4 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⁻¹⁹ C is the charge of an electron, m_e = 9.1 × 10⁻³¹ kg is the mass of an electron, k_B = 1.38 × 10⁻³¹ J K⁻¹ is the Boltzmann constant, and ε₀ = 8.85 × 10⁻¹² F m⁻¹ the vacuum permittivity.

4.1 Simulation in vacuum

The N point charges 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₁ and T₂ operate in counter-phase (push-pull) and for which the charges are mainly distributed on those transmitters, 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.

Fig. 4 Vacuum charge distribution on the probe operating in antiphase mode with a half-level of emission according to the DSCD model. Left: . Right: black: φi = π and red: φi = 0.

4.2 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₁ = 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 inFig. 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 m, T_h ∕T_c , and 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 ϕ∕ϕ₀, which is always found in a range of values close to the inverse of the cold plasma permittivity ε_c = 1 − Ω⁻² 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 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₀, 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₁ that is lower than f₀ 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 .

Fig. 5 Sensitivity of the simulated response of the probe RPC-MIP, operated in SDL, in a double-Maxwellian evdf plasma. The probe operates with a full emission level and is isolated from the plasma by the vacuum sheath presented in Fig. 3. Top left: influence of R. λDh = 1 m, nh ∕ntot = 0.3 and Th ∕Tc = 20. Top right: influence of λDh. R = 1 m, nh ∕ntot = 0.2 and Th ∕Tc = 10. Bottom left: influence of Th∕Tc. λDh = 1 m, nh ∕ntot = 0.2 and R = 1 m. Bottom right: influence of nh∕ntot. λDh = 1 m, Th ∕Tc = 20, and R = 1 m.

Fig. 6 Evolution of the level of the response of the probe RPC-MIP, operated in SDL vs. the sheath radius R for Ω = 2 in the so-called phase mode (). A fluctuation of ± 0.4 dB (lower than the noise level of RPC-MIP) of the results is observed when the Debye length of the hot electron population, the temperature, and density ratio of the electron populations are varied in the ranges m, , and . The dashed line is the level in vacuum.

Fig. 7 Modulus of the function ϕ∕ϕ0 normalized by the inverse of the cold plasma permittivity εc = 1 − Ω−2 for Ω = 0.6 (top panel) and Ω = 2 (bottom panel) in a double-Maxwellian evdf plasma. r is the distance from the considered pulsating point charge. It corresponds to rij in Eq. (11) and rkm in Eq. (13).

Fig. 8 Response of the probe RPC-MIP, operated in SDL in the so-called phase mode () when the sheath is neglected in the model. λDh = 1m, nh ∕ntot = 0.2, and Th ∕Tc = 20. The influence of the sheath can be appreciated by comparing this figure with Fig. 5 (bottom left).

5 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 concludethat 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 (17)

where the energy E and the temperatures T_c,h are in eV.

Fig. 9 Top panel: experimental spectral power acquired by the probe RPC-MIP operating in SDL mode on 20 March 2016 between 4:10 and 5:10 (UTC). Bottom panel: simulated spectral power with our model considering a vacuum sheath around the probe and a double-Maxwellian plasma. The input parameters for the simulation are nh = 200 cm−3, nc = 800 cm−3, Th = 3.61 eV, Tc = 0.36 eV, and R = 120 cm. The dashed line curves correspond to the model when the sheath is neglected (see Sect. 3.2).

Fig. 10 Fitting our simulations on an experimental spectral power of the signal provided by the probe operating in phase mode () on 6 May 2016. The x-axis has been normalized by 399 kHz. The distances from the Sun of comet 67P and the altitudes of the spacecraft were 2.95 AU and 9.95 km, respectively. The input parameters for the simulation are nc = 1572 cm-3, nh = 393 cm-3, Tc = 0.25 eV, Th = 10.22 eV, and R = 120 cm.

Fig. 11 Fitting our simulations on an experimental spectral power of the signal provided by the probe operating in phase mode () on 24 August 2016. The x-axis has been normalized by 196 kHz. The distances from the Sun of comet 67P and the altitudes of the spacecraft were 3.61 AU and 15.03 km, respectively. The input parameters for the simulation are nc = 560 cm-3, nh = 140 cm-3, Tc = 0.1 eV, Th = 4.95 eV, and R = 130 cm.

Fig. 12 Fitting our simulations on an experimental spectral power of the signal provided by the probe operating in phase mode () on 20 March 2016. The x-axis has been normalized by 343 kHz. The distances from the Sun of comet 67P and the altitudes of the spacecraft were 2.62 AU and 12.07 km, respectively. The input parameters for the simulation are nc = 1162 cm-3, nh = 290 cm-3, Tc = 0.29 eV, Th = 8.85 eV, and R = 110 cm.

Fig. 13 Fitting our simulations on an experimental spectral power of the signal provided by the probe operating in phase mode () on 20 March 2016. The x-axis has been normalized by 427 kHz. The distances from the Sun of comet 67P and the altitudes of the spacecraft were 2.62 AU and 12.07 km, respectively. The input parameters for the simulation are nc = 1800 cm-3, nh = 450 cm-3, Tc = 0.16 eV, Th = 11.70 eV, and R = 110 cm.

Fig. 14 Electron energy distribution functions corresponding to the simulations displayed in Figs. 10–13. The integral of these functions provides the overall electron density ntot.

6 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 displays the estimated electron densities and electron temperatures of the estimated evdf on 20 March 2016 according to the assumptions we considered. Similarly, Fig. 16 displays the same parameters corresponding to 24 August 2016.

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 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.

Fig. 15 Estimate of the evdf on 20 March 2016. Top panel: overall electron density ntot (uncertainty 22%), hot electron density nh (uncertainty 40%), and cold electron density nc (uncertainty 45%). Bottom panel: temperature of the cold (uncertainty 87%) and the hot electron population temperature (uncertainty57%).

Fig. 16 Estimate of the evdf on 24 August 2016. Top panel: overall electron density ntot (uncertainty 22%), hot electron density nh (uncertainty 40%), and cold electron density nc (uncertainty 45%). Bottom panel: temperature of the cold (uncertainty 87%) and the hot electron population temperature (uncertainty57%).

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 electrostaticmodeling. As observed during the Rosetta mission, the plasma signature on the simulated response of the probeis 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 ϕ∕ϕ₀ 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 ϕ∕ϕ₀ 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.

8 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₂O and CO₂).

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.

Appendix A Computation of α_i

A.1 For a cylinder in vacuum

The center of the cylinder facets is located in cylindrical coordinates R, θ, and z, R being the radius of the cylinder. Every facet carries a single charge q_i and is divided in 16 sub-surfaces, each of which carries a sub-charge equal to q_i ∕16 (see Fig. A.1). The potential V_ii induced by the sub-charges at the center of the facet is then

Fig. A.1 Cylinder facet divided into sub-surfaces with associated sub-charges considered to compute αi.

Thus (A.3)

with and

A.2 For a sphere in vacuum

The same approach can be adopted for a spherical object in vacuum. A facet of this object centered on the spherical coordinates R, θ, and φ and carrying a single charge q_i has an area S_i given by

where Δθ and Δφ correspond to the angular steps in spherical coordinates used to build the object. We note that θ is the azimuthal angle, φ is the polar angle, and the origin of the coordinate system is located in the center of the sphere. After dividing the facet into 16 sub-surfaces, the coefficient α_i reads

with d_a,b,c,d,e,f the distances between the sub-charges and the center of the facet, (A.10)

with x = a, b, c, d (A.11)

and S_a,b,c,d,e,f the area of the sub-surfaces of the facets, (A.12) (A.13) (A.14) (A.15)

where the coordinates of the sub-charges reads (A.16) (A.17) (A.18)

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.

Fig. A.2 Point charge in front of an insulated conducting sphere.

Fig. A.3 Electrostatic potential computed along a line parallel to the x-axis according to the DSCD method and by the method of electrostatic images for the configuration of Fig. A.2 with d = 2a, q = 1 C, y = 0, and z = 2a. The potential at the surface of the sphere reaches 25 V according to the computations.

Appendix B Computation of β_k

β_k is analogous to α_i, but for a cylinder or a sphere in a homogeneous infinite plasma. We recall that the function ϕ∕ϕ₀ corresponds to the inverse of the plasma relative permittivity function, which depends on the distance from the corresponding point charge, the electron velocity distribution function (evdf), the electron-to-neutral collisional rate, the magnetic field, and the operating frequency. 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, the functionhas also been tabulated in the same conditions, but considering a double-Maxwellian evdf (Gilet et al. 2017).

B.1 For a cylinder in homogeneous plasma

The cylinder described in Appendix A.1 is now immersed in homogeneous infinite plasma and Eq. (A.2) reads (B.1)

where (B.2)

with and (B.3) (B.4)

B.2 For a sphere in homogeneous plasma

For a sphere immersed in homogeneous plasma the expression of β_i in Eq. (B.1) reads (B.7)

with d_a,b,c,d,e,f the distances between the sub-charges and the center of the facet, (B.8)

with x = a, b, c, d (B.9)

and S_a,b,c,d,e,f the area of the sub-surfaces of the facets, (B.10) (B.11) (B.12) (B.13)

where the coordinates of the sub-charges reads (B.14) (B.15) (B.16)

Appendix C Continuity of the potential and the Maxwell-Gauss equation at the sheath interface

The expression of the potential on both sides of the interface reads (C.1) (C.2)

with V_1k and V_2k the potentials on the internal side (V_1k) and on the external side (V_2k) of an elementary surface S_k located at the plasma-sheath interface. 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 other charges q_1m of the internal side of the plasma-sheath interface. Then the continuity of the potential (V_1k = V_2k) provides Eq. (13).

The sheath edge representing an interface between two dielectric media, the free charge density on the interface is null, leading to a conservation of the perpendicular component of the displacement vector (Maxwell-Gauss equation), (C.3)

with the electric field at the internal side of the interface and the electric field at the external side of this last. is the normal vector at the interface. The subscript v stands for “vacuum”, meaning that q_2k are considered in vacuum (see Béghin & Kolesnikova 1998 for more details). The perpendicular components of the displacement vector are determined according to the difference between the potentials at the internal side and at the external side of the interface divided by 2 × d r. The electric field induced on S_k by its own charge (q_1k for the internal and q_2k for the external side) is also taken into account through − q_1k∕2S_k and q_2k∕2S_k, (C.4) (C.5)

with r_ki1 and r_ki2 the distances between the kth point chargeq_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. Then D_1k⊥ = D_2k⊥ provides Eq. (14).

References

All Tables

All Figures

	Fig. 1 Schematic of the mutual impedance probe mounted on its boom.
In the text

	Fig. 2 Equivalent electrical arrangement of the emitters and receivers. This schematic does not correspond to the geometry considered in the electrostatic simulation (see Fig. 3).
In the text

	Fig. 3 Spatial configuration considered in the electrostatic simulations. 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 presentat the external side of the sheath, which is presented in sectional view in order to show the probe.
In the text

	Fig. 4 Vacuum charge distribution on the probe operating in antiphase mode with a half-level of emission according to the DSCD model. Left: . Right: black: φi = π and red: φi = 0.
In the text

Fig. 5 Sensitivity of the simulated response of the probe RPC-MIP, operated in SDL, in a double-Maxwellian evdf plasma. The probe operates with a full emission level and is isolated from the plasma by the vacuum sheath presented in Fig. 3. Top left: influence of R. λDh = 1 m, nh ∕ntot = 0.3 and Th ∕Tc = 20. Top right: influence of λDh. R = 1 m, nh ∕ntot = 0.2 and Th ∕Tc = 10. Bottom left: influence of Th∕Tc. λDh = 1 m, nh ∕ntot = 0.2 and R = 1 m. Bottom right: influence of nh∕ntot. λDh = 1 m, Th ∕Tc = 20, and R = 1 m.

In the text

Fig. 6 Evolution of the level of the response of the probe RPC-MIP, operated in SDL vs. the sheath radius R for Ω = 2 in the so-called phase mode (). A fluctuation of ± 0.4 dB (lower than the noise level of RPC-MIP) of the results is observed when the Debye length of the hot electron population, the temperature, and density ratio of the electron populations are varied in the ranges m, , and . The dashed line is the level in vacuum.

In the text

	Fig. 7 Modulus of the function ϕ∕ϕ0 normalized by the inverse of the cold plasma permittivity εc = 1 − Ω−2 for Ω = 0.6 (top panel) and Ω = 2 (bottom panel) in a double-Maxwellian evdf plasma. r is the distance from the considered pulsating point charge. It corresponds to rij in Eq. (11) and rkm in Eq. (13).
In the text

	Fig. 8 Response of the probe RPC-MIP, operated in SDL in the so-called phase mode () when the sheath is neglected in the model. λDh = 1m, nh ∕ntot = 0.2, and Th ∕Tc = 20. The influence of the sheath can be appreciated by comparing this figure with Fig. 5 (bottom left).
In the text

Fig. 9 Top panel: experimental spectral power acquired by the probe RPC-MIP operating in SDL mode on 20 March 2016 between 4:10 and 5:10 (UTC). Bottom panel: simulated spectral power with our model considering a vacuum sheath around the probe and a double-Maxwellian plasma. The input parameters for the simulation are nh = 200 cm−3, nc = 800 cm−3, Th = 3.61 eV, Tc = 0.36 eV, and R = 120 cm. The dashed line curves correspond to the model when the sheath is neglected (see Sect. 3.2).

In the text

Fig. 10 Fitting our simulations on an experimental spectral power of the signal provided by the probe operating in phase mode () on 6 May 2016. The x-axis has been normalized by 399 kHz. The distances from the Sun of comet 67P and the altitudes of the spacecraft were 2.95 AU and 9.95 km, respectively. The input parameters for the simulation are nc = 1572 cm-3, nh = 393 cm-3, Tc = 0.25 eV, Th = 10.22 eV, and R = 120 cm.

In the text

Fig. 11 Fitting our simulations on an experimental spectral power of the signal provided by the probe operating in phase mode () on 24 August 2016. The x-axis has been normalized by 196 kHz. The distances from the Sun of comet 67P and the altitudes of the spacecraft were 3.61 AU and 15.03 km, respectively. The input parameters for the simulation are nc = 560 cm-3, nh = 140 cm-3, Tc = 0.1 eV, Th = 4.95 eV, and R = 130 cm.

In the text

Fig. 12 Fitting our simulations on an experimental spectral power of the signal provided by the probe operating in phase mode () on 20 March 2016. The x-axis has been normalized by 343 kHz. The distances from the Sun of comet 67P and the altitudes of the spacecraft were 2.62 AU and 12.07 km, respectively. The input parameters for the simulation are nc = 1162 cm-3, nh = 290 cm-3, Tc = 0.29 eV, Th = 8.85 eV, and R = 110 cm.

In the text

Fig. 13 Fitting our simulations on an experimental spectral power of the signal provided by the probe operating in phase mode () on 20 March 2016. The x-axis has been normalized by 427 kHz. The distances from the Sun of comet 67P and the altitudes of the spacecraft were 2.62 AU and 12.07 km, respectively. The input parameters for the simulation are nc = 1800 cm-3, nh = 450 cm-3, Tc = 0.16 eV, Th = 11.70 eV, and R = 110 cm.

In the text

	Fig. 14 Electron energy distribution functions corresponding to the simulations displayed in Figs. 10–13. The integral of these functions provides the overall electron density ntot.
In the text

	Fig. 15 Estimate of the evdf on 20 March 2016. Top panel: overall electron density ntot (uncertainty 22%), hot electron density nh (uncertainty 40%), and cold electron density nc (uncertainty 45%). Bottom panel: temperature of the cold (uncertainty 87%) and the hot electron population temperature (uncertainty57%).
In the text

	Fig. 16 Estimate of the evdf on 24 August 2016. Top panel: overall electron density ntot (uncertainty 22%), hot electron density nh (uncertainty 40%), and cold electron density nc (uncertainty 45%). Bottom panel: temperature of the cold (uncertainty 87%) and the hot electron population temperature (uncertainty57%).
In the text

	Fig. A.1 Cylinder facet divided into sub-surfaces with associated sub-charges considered to compute αi.
In the text

	Fig. A.2 Point charge in front of an insulated conducting sphere.
In the text

	Fig. A.3 Electrostatic potential computed along a line parallel to the x-axis according to the DSCD method and by the method of electrostatic images for the configuration of Fig. A.2 with d = 2a, q = 1 C, y = 0, and z = 2a. The potential at the surface of the sphere reaches 25 V according to the computations.
In the text