Boltzmann electron PIC simulation of the E-sail effect

The solar wind electric sail (E-sail) is a planned in-space propulsion device that uses the natural solar wind momentum flux for spacecraft propulsion with the help of long, charged, centrifugally stretched tethers. The problem of accurately predicting the E-sail thrust is still somewhat open, however, due to a possible electron population trapped by the tether. Here we develop a new type of particle-in-cell (PIC) simulation for predicting E-sail thrust. In the new simulation, electrons are modelled as a fluid, hence resembling hydrid simulation, but in contrast to normal hybrid simulation, the Poisson equation is used as in normal PIC to calculate the self-consistent electrostatic field. For electron-repulsive parts of the potential, the Boltzmann relation is used. For electron-attractive parts of the potential we employ a power law which contains a parameter that can be used to control the number of trapped electrons. We perform a set of runs varying the parameter and select the one with the smallest number of trapped electrons which still behaves in a physically meaningful way in the sense of producing not more than one solar wind ion deflection shock upstream of the tether. By this prescription we obtain thrust per tether length values that are in line with earlier estimates, although somewhat smaller. We conclude that the Boltzmann PIC simulation is a new tool for simulating the E-sail thrust. This tool enables us to calculate solutions rapidly and allows to easily study different scenarios for trapped electrons.


Introduction
The electric solar wind sail (electric sail, E-sail) is a planned device for producing interplanetary spacecraft propulsion from the natural solar wind by a set of centrifugally stretched thin metallic tethers that are kept artificially at a high positive potential . In order to engineer E-sail devices in detail, one should predict the magnitude of the Coulomb drag that the flowing solar wind exerts on the charged tether. Thus far this problem has been studied using particle-in-cell (PIC) simulations (Janhunen and Sandroos, 2007;Janhunen, , 2012, Vlasov simulations (Sánchez-Arriaga and Pastor-Moreno, 2014) and other methods (Sanmartín et al., 2008;Sanchez-Torres, 2014). However, a fully satisfactory way of estimating E-sail thrust has not yet emerged, except for negative polarity tethers (Janhunen, 2014a), which are however more suitable to use in low Earth orbit (LEO) as a deorbiting plasma brake device  than in the solar wind . Independent laboratory measurements of the width of the electron sheath around a positively biased tether placed in streaming plasma were made by Siguier et al. (2013). The laboratory results were produced in a plasma mimicking conditions in LEO and they are in good agreement with PIC thrust predictions (Janhunen, 2014b).
In this paper we develop a new variant of the PIC simulation of Janhunen (2012) which treats electrons as a fluid obeying the Boltzmann relation, with two additional tricks to be detailed below. The motivation is to have a code which is relatively fast and easy to run and which makes it possible for the user to control the assumed amount of trapped electrons in a simple way. We present the simulation code and use it to derive a thrust estimate in a case which has relevance to the E-sail. Making more comprehensive set of runs with different voltages and in different solar wind conditions is outside of the scope of the paper.

Physics of tether Coulomb drag
The negative polarity Coulomb drag effect can be successfully simulated by a PIC code (Janhunen, 2014a). The negative polarity effect is less challenging to simulate than the positive polarity one, because in the negative polarity case, trapped particles do not form. Electrons are not trapped because they are repelled by the tether's negative potential structure. Ions (protons) are also not trapped, because in rel-2 P. Janhunen: E-sail Boltzmann PIC simulation evant cases the surrounding plasma flow is supersonic so that in a coordinate frame where the tether is stationary, ions enter the potential structure with significant kinetic energy and hence cannot be easily trapped by the potential.
In the positive polarity case which is the subject of the present paper, electrons may become trapped because the tether attracts them and because the electron flow is subsonic so that there are some electrons which enter the potential structure with nearly zero initial energy. The issue of trapped electrons was identified by Janhunen and Sandroos (2007) and simulations were later made where chaotisation of trapped electron orbits when scattering at the end of the tether and at the spacecraft was explicitly included in the PIC code, as well as occasional removal of trapped electrons due to collisions with the tether (Janhunen, 2012). However, in that calculation the timescale where trapped electrons were removed was made much faster than in reality, in order to complete the simulation run in a reasonable time. Basically, processes controlling the size of the trapped electron population are poorly known, because we think that they happen at longer timescales than what is accessible by PIC simulation.
Consider a 2-D positive potential structure which thus forms a negative potential energy well (attractive structure) for electrons. If the potential is stationary and hosts no trapped electrons and if the ambient electron distribution is a non-drifting Maxwellian and the plasma is tenuous enough to be collisionless, then it can be shown using Liouville's theorem that the local electron density can nowhere exceed the background plasma density (Laframboise and Parker, 1973;). In fact, under the stated conditions the electron density is typically significantly less than the background density in most regions covered by the attractive potential structure. Physically the result is due to the fact that although the tether attracts electrons and therefore focusses them to move through its vicinity, the electron density remains low because accelerating the electron reduces the time that the electron spends inside the attractive potential structure. Conservation of the electron's originally nonzero angular momentum also limits the minimum distance that the electron can have from the tether.
That the electron density nowhere exceeds the background density leads to an interesting dilemma (Éric Choinière, private communication) that is associated with trapped electrons because the ion density, on the other hand, must be higher than the background density inside at least in some regions of the upstream side ion deflection shock. Then, a significant positive charge density can build up at the upstream ion deflection region, which spreads out the electronattractive potential structure and thus decreases the local electron density further which, unless some process intervenes, makes the charge imbalance worse. What exactly happens remains not well known, but conceivably, instabilities which are able to scatter electrons might develop and take place until enough trapped electrons have been formed so that proper charge neutralisation in the ion deflection shock region has been restored. A simple approach to model this kind of general behaviour may be to reduce the number of trapped electrons in the simulation until the result becomes nonphysical in some way or another and postulate that the most realistic run is the physical run with least amount of trapped electrons. This is what the simulation presented next attempts to accomplish.

Simulation code
We use a special version of electrostatic PIC where electrons are not modelled as particles, but as a fluid whose local density n e (x, y) is an assumed function of the local potential V (x, y). When the potential repels electrons (i.e. when V (x, y) < 0), the Boltzmann relation is assumed, where n o is the background plasma density and T e = 10 eV is the background electron temperature in energy units. The Boltzmann relation is supposedly a good approximation for the density of particles that try to penetrate thermally into a repulsive potential structure. For attractive parts of the potential, the electron density depends on the number and distribution function of trapped electrons and thus extra assumptions are needed. That one needs such assumptions is natural because we are not modelling electron dynamics using first principles. We want to have a functional relationship which contains a parameter that can be used to tune the number of trapped electrons in a simple way.
In this paper we use the following functional form for electron density in attractive parts of the potential: where ν > 0 is the trapped electron control parameter.
In Figure 1 we show how the local electron density depends on the local potential according to Eqs. (1) and (2).
The code that we use in this paper is a descendant of a PIC code that we used earlier (Janhunen and Sandroos, 2007;Janhunen, 2012Janhunen, , 2014a. It is an electrostatic code with linear particle weighting and additional grid-level smoothing which is implemented in the Fourier domain. The code supports 2-D and 3-D simulations and optionally includes a constant external magnetic field. The Poisson solver is implemented by fast Fourier transform and the code is parallelised with OpenMP threads and the Message Passing Interface (MPI) library. The code is also vectorised, so that for a processor with 256-bit wide AVX2 instruction set, each processor core calculates eight particle updates simultaneously. For particle coordinates and velocities, 32-bit floating point accuracy is used. To avoid loss of accuracy, particle coordinates are expressed relative to the centroid of the grid cell where the particle resides. Particles are stored in lists belonging to grid cells which enables fast accumulation of the charge density (no irregular memory addressing patterns needed) and enables good data locality and cache friendliness. When a particle crosses a grid cell boundary, it is removed from the old list and entered into a new one.
An additional trick is needed to avoid solutions that look unphysical and oscillate drastically in time. Namely, the electron density must not react instantaneously to changed electric potential, but with some delay. If n * e (x, y) is the preliminary instantaneous electron density computed by Eqs. (1) and (2), then we propagate the real electron density n e (x, y) by the differential equation where τ is a time constant. The solution n e (t) of (3) is effectively a low-pass filtered version of n * e (t) where frequencies above ∼ 1/τ are filtered out. We have tested different values of τ and found that in order to produce good results (in the sense of avoiding unphysical oscillations), τ should be proportional to the ion plasma oscillation period, because shorter timescale phenomena are not included in a model with electron fluid. Specifically, we use where ω pi is the ion plasma frequency By numerical experimentation, C ≤ 1 produces unphysical oscillations which disappear if C > 1. We use C = 1.5 to have some extra margin of safety. On the other hand there is no reason to increase C further. In summary, the code pushes ions forward by their equations of motion (in this paper with B=0), accumulates ion density n i from the ion positions where S is the linear area weighing function (Birdsall and Langdon, 1991). Then the code applies Fourier domain filtering to smooth the ion density spatially, calculates the provisional electron density n * e by Eqs. (1) and (2), updates the electron density n e by Eq. (3), computes the charge density ρ e = e(n i − n e ), solves the Poisson equation computes the gridded electric field E(x) from φ(x) by finite differencing and repeats the process for the next timestep. The particle equations of motion (6) are discretised in time using the second order accurate leapfrog method as is common in momentum conservative explicit PIC simulations (Birdsall and Langdon, 1991). Equation (3) is discretised by forward Euler method. This is equivalent to linearly mixing the portion ∆t/τ of n * e into n e at each timestep, where ∆t is the length of the timestep. Table 1 summarises the simulation parameters.
The timestep ∆t in our code is ∼ 50 times longer than what it would be in a full electron-ion PIC code using the same grid resolution. In addition to that, it seems that the code can produce useful results with a smaller number of ions per grid cell than a full PIC code so that the improvement factor in practical computational efficiency is probably several hundred. The latter property is probably due to averaging over timescale τ ≈ 210 ∆t which effectively reduces the numerical noise.

Results
First we made a family of runs at low resolution, varying the trapped electron parameter ν between 0.16 and 0.08. The grid spacing was 10 m and the box side length 1.92 km so that the grid was 192 × 192. The timestep was 4 µs, corresponding to 2500 km/s maximum signal speed, resulting in Courant-Friedrichs-Lewy (CFL) parameter of 0.16 with respect to the used solar wind speed of 400 km/s ( Table 1). The solar wind plasma density is 7.3 cm −3 and for simplicity, no interplanetary magnetic field is used. The ion and electron temperature is 10 eV. Pure proton plasma is used as the solar wind. The plasma parameters correspond to average solar wind properties at 1 au distance from the sun.
The real solar wind also has typically 10 nT magnetic field. However, the resulting background electron Larmor radius of ∼1 km is much larger than the diameter of the electron sheath ∼100-200 m. Thus the magnetic field is expected to play only a minor role as far as the magnitude of the E-sail thrust is concerned, so that assuming zero magnetic field is not a bad first approximation.
For each run we determined the thrust by two complementary methods, "Mom" and "Coul". In method Mom we keep track of the momentum x component of all particles entering and leaving the simulation box and deduce that the momentum which is left in the box is the one which gets transferred to the tether and the spacecraft. In method Coul we evaluate the negative gradient of the plasma potential at the tether's location and use Coulomb's law, dF/dz = λE where dF/dz is the thrust per tether length, λ is the tether's linear charge density and E is the electric field caused by the plasma and evaluated at the tether. The line charge λ is given as an input parameter (in the runs performed, λ = 8.64·10 −8 As/m), and the corresponding tether voltage V 0 is calculated from the run results after the end of the run by summing the tether's vacuum potential and the potential created by the plasma. Because the plasma potential varies somewhat from run to run, the resulting tether voltage V 0 is also slightly different in each run. For each run, we made a qualitative visual check of the result to determine if the solution was physically reasonable in the sense that it contained not more than one ion deflection shock.
We then made an identical family of runs with halved grid spacing of 5 m. The computing time is about 8 times longer than in the low resolution set. Results of both families of runs are summarised in Table 2.
Run ν = 0.12 is selected as the baseline because our visual inspection indicates that ν = 0.12 is the smallest value of ν for which we get not more than one ion deflection shock. Figure 2 shows the electron and ion density and the total potential in the baseline ν = 0.12 run as well as in another run where ν = 0.08 and which exhibits two ion shocks. In addition, in cases where two shocks emerge in the solution, the shocks also tend to slowly propagate upstream so that full convergence is not reached. The high resolution run data are used, but for viewing speed convenience the spatial grid resolution has been halved by neighbour averaging before plotting in Fig. 2. The solar wind arrives from the right in Fig. 2.
In Fig. 2 one clearly sees that the run ν = 0.12 looks physically meaningful and well-behaved, while run ν = 0.08 exhibits two shocks and is thus physically not acceptable. Similar visual judgement can be made also for the other performed runs, and information presented for each run in Table  2 was compiled in this way. Figure 3 shows the thrust history of the ν = 0.12 run (high resolution version). After natural initial transients have died out, the Coul and Mom methods are in good agreement with each other regarding the thrust per length produced by the E-sail. By construction, the Mom method reacts to changes slower than Coul, and during the first ∼ 10 ms both methods give unreliable results.

Discussion and conclusions
Using Boltzmann electrons in PIC simulation with two additional tricks enables us to rapidly compute reasonable looking solutions for the positive polarity Coulomb drag problem, applicable to the E-sail in the solar wind. The first trick is that while the Boltzmann relation (1) is a natural choice for the electron density in repulsive parts of the potential, in the attractive parts we use an a power law relationship (2) containing one free parameter that can be used to control the amount of trapped electrons. The second trick is that to avoid spurious temporal oscillations, the electron density is not updated instantaneously, but frequencies higher than ∼ ω pi are filtered out by Eq. (3).
The calculated solutions depend on parameter ν that must be specified by the user and that parametrises the number of trapped electrons. Once ν is fixed, the simulation produces a prediction for the thrust per tether length. We calculated the thrust per length by two complementary methods and found good mutual agreement. For too small values of ν the result looks unphysical, however, because it contains two ion deflection shocks. Our strategy is that by finding the smallest value of ν which gives a physically reasonable single shock solution, we obtain a solution which is expectedly the best approximation to reality within this framework.
For voltage V 0 of 17.7 kV and nominal solar wind at 1 au (density 7.3 cm −3 , speed 400 km/s), the thrust predicted by the framework is ∼ 320 nN/m. In our earlier paper (Janhunen, 2009, Figure 12) we arrived at estimate 500 nN/m at V 0 = 20 kV if there are no trapped electrons, 400 nN/m if trapped density equals background density, and 320 nN/m if the average trapped density in the sheath is four times higher than background. If one assumes that the thrust is approximately linearly proportional to V 0 , after scaling to 20 kV the new estimate becomes ∼ 360 nN/m. The result is thus in quantitative agreement with Janhunen (2009) when one takes into account that the average trapped density is ∼ 1.5 in the new simulation. In summary, if one assumes that the number of trapped electrons increases until the solution has only one ion deflection shock, then one ends up with ∼ 360 nN/m thrust estimate for 20 kV and average 1 au solar wind.
One must be somewhat cautious when interpreting the numerical result, because the electron Boltzmann simulation in any case contains less physics than the earlier full PIC simulations, and, for example, the particular functional form (2) used to describe the electron density in attractive parts of the potential was chosen mostly by mathematical convenience. In any case, the Boltzmann simulation is a new addition in the toolbox of the E-sail thrust analyst. Its benefits are fast computation compared to traditional PIC and easy possibility to test different assumptions concerning trapped electrons.