Thrust calculation of electric solar wind sail by particle-in-cell simulation

. In this study, thrust characteristics of an electric solar wind sail were numerically evaluated using full three-dimensional particle-in-cell (PIC) simulation. The thrust obtained from the PIC simulation was lower than the thrust estimations obtained in previous studies. The PIC simulation indicated that ambient electrons strongly shield the electrostatic potential of the tether of the sail, and the strong shield effect causes a greater thrust reduction than has been obtained in previous studies. Additionally, previous expressions of the thrust estimation were modiﬁed by using the shielded potential structure derived from the present simulation re-sults. The modiﬁed thrust estimation agreed very well with the thrust obtained from the PIC simulation.


Introduction
An electric solar wind sail, called the "E-sail," is a recently proposed propulsion device that consists of 50-100 conductive tethers with lengths of 10-20 km and thicknesses of 0.1-1 µm.The E-sail was first proposed by Janhunen (2004).The main body of the spacecraft expands the tethers to form a sail-like structure.The E-sail has electron guns to maintain a positive surface potential on the order of several kilovolts, in order to deflect solar wind protons.The tethers obtain the momentum of these deflected protons via Coulomb scattering and use it as their propulsive force.The system requires electron sources and electrical power for the electron guns to produce thrust.The E-sail is expected to be used as a new propellantless space propulsion device.
The thrust characteristics of the E-sail were first investigated by Janhunen and Sandroos (2007).They performed a one-dimensional (1-D) particle-in-cell (PIC) simulation of a conductive tether with a radius of 1.0 m.They found that an ansatz of the electrostatic potential structure around the tether, which is expressed as agreed very well with the result of their PIC simulation, where r is the distance from the tether, V 0 is the surface potential of the tether, λ D = ε 0 k B T e /e 2 n e is the Debye length of the electron, and r w is the radius of the tether.For evaluating the performance of an E-sail, the thrust per unit length is often used.The total thrust can be calculated from the thrust per unit length by multiplying the number of tethers by the length of one tether.Janhunen and Sandroos (2007) also conducted a two-dimensional (2-D) PIC simulation and suggested that the thrust per unit length acting on the tether is where L is the length of the tether, K is the coefficient of proportionality (K ≈ 3.09 was obtained from their Monte Carlo simulation), m p is the mass of a proton, n 0 is the ambient plasma density, and v d is the drift velocity of the solar wind.Moreover, Janhunen (2009) suggested the thrust of the E-sail may increase because of a lack of electrons around the tether.He proposed the use of multiple tethers to collect ambient electrons so that the electron density around each tether decreases and the ambient electrons cannot completely shield the potential of the tether.According to Janhunen (2009), the thrust with ambient electron removal is 5 times larger than that obtained by Janhunen and Sandroos (2007) without ambient electron removal.Sanchez-Torres (2014) also investigated the thrust of the E-sail considering Coulomb scattering, assuming the absence of trapped electrons around the tether.He expressed the potential structure as where the approximated parameter b is 0.65 for a potential bias of V 0 =10-40 kV and r sh is the sheath radius for a highly positive bias tether and can be calculated from the ambient plasma parameters and V 0 (Sanmartín et al., 2008).Using the solar wind parameters at 1 AU, the thrust per unit length is 407 nN m −1 for a 20 kV charged tether of r w = 20 µm.This is lower than the thrust per unit length of 500 nN m −1 estimated by Janhunen (2009) but higher than the value of 100 nN m −1 estimated by Janhunen and Sandroos (2007) because Eq. ( 3) yields greater values of the potential than Eq.( 1) under almost all conditions.The study by Sanchez-Torres (2014) was purely analytical; no plasma simulations were performed.However, Hoshi et al. (2016) showed that the actual electrostatic potential structure around the tether was lower than those given by Eqs.(1) and (3) based on the results of a full three-dimensional (3-D) PIC simulation with V 0 = 240 V.The thrust of the E-sail is generated from the deflection of solar wind protons by the electrostatic potential.If the electrostatic potential derived from the tethers is greatly shielded by ambient electrons, the actual thrust is lower than that estimated in previous studies.Hoshi et al. (2016) did not consider the thrust because a potential of 240 V was not sufficient to deflect solar wind protons with a drift velocity of approximately 400 km s −1 (≈ 0.8 keV).
In the present paper, we performed 3-D full PIC simulations to simulate a transient of the thrust of the E-sail.The thrust found in this paper is lower than that obtained in previous studies, with a sufficiently high potential to deflect ambient protons (V 0 > 1 kV).Section 2 discusses the 3-D PIC simulation with V 0 ≤ 4.0 kV that was performed to confirm that the potential was lower than those obtained from Eqs. (1) and (3).The propulsive force that acts on the tether is also calculated in the PIC simulation.In Sect.3, the thrust is numerically estimated and compared with the PIC results.Two estimation procedures employed in previous studies (Janhunen and Sandroos, 2007;Sanchez-Torres, 2014) are modified to contain a shielded potential structure derived from the  present simulation results.The estimated thrust and the PIC results are found to be in good agreement with each other and lower than those estimated in previous studies.
2 Full PIC simulation of E-sail

Simulation settings
This section describes the PIC simulation configurations of a positively charged tether in the solar wind environment.The simulation code HiPIC, which was developed by the Japan Aerospace Exploration Agency's Engineering Digital Innovation Center (JEDI) (Muranaka et al., 2011), was used to perform this simulation.HiPIC is an electrostatic code that models 3-D rectangular cells in space and uses the full PIC method to calculate collisionless kinetic plasma.HiPIC solves Newton's equations of motion for each particle using the Buneman-Boris method and solves Poisson's equations to obtain the electric potential structure in the computational domain using a discrete sine transformation.HiPIC can be used to calculate the interaction between plasmas and the spacecraft, which is modeled with rectangular internal boundaries.A detailed description of the performance of the code is given in Muranaka et al. (2011).
The simulation and physical parameters are given in Table 1.The electron density n e and the proton density n p are n e = n p = n 0 = 1.0×10 7 m −3 , and their temperatures are k B T e = 100 eV and k B T p = 12 eV, respectively.The background plasmas have a solar wind drift velocity of v d = 400 km s −1 along the x axis.All of the edges of the simulation domain were fixed to V = 0V (Dirichlet boundary condition).The background plasma particles were injected from all the domain boundaries in each time step as many times as the number of outgoing particles in previous time step.
A tether-like rectangular model was set in the center of the computational domain, as shown in Fig. 1.The size of the model is 1×200×1 cells, which is equivalent to a tether with dimensions of r w = 15 cm and L = 60 m.Tethers with L = 30 and 75 m were also simulated to confirm the influence of the length of the tether on the thrust per unit length.
Because of the limitation of the calculation resources, we cannot include the emission of the electron beam from the tether's edge and simulate the self-charging of the tether.Instead of emitting the electron, the surface potential V 0 of the tether was fixed to an inputted value.Hoshi et al. (2016) showed that effects of emitted electrons on the potential structure were small, so we consider that the absence of emitted electrons do not cause significant differences in the force acting on the tether.
V 0 was varied from 0 to 4.0 kV, and the thrust acting on the tether F x was calculated.F x is the x component of the thrust and was calculated as the sum of the total momentum of the particles impinging on the tether during each time step divided by t and the Coulomb force calculated from the Maxwell stress tensor.The simulation progressed with time steps of t = 20 ns until the time variation of the external force became zero.In almost all the cases, the total iteration was 10 000 steps (= 0.2 ms).An additional 2000 steps were calculated for the V 0 = 4.0 kV case.
To perform the computation, we applied an MPI parallelization and an OpenMP parallelization.Each case used 2048 cores on Cray XE4 for calculation (1024 processes for MPI, two threads for OpenMP), requiring approximately 12 h to run.

Simulation results
Figure 2a shows the time variation of F x at V 0 = 1.3 kV.At the beginning of the simulation, F x was almost zero because the ambient particles had not yet begun to respond to the electrostatic potential of the tether.F x then increased with time and converged to a specific value.At the end of the simulation, F x was 0.35, 0.76, and 0.96 µN for L = 30, 60, and 75 m, respectively.These values represent the total force acting on the tether, including the sum of the particles hitting the whole tether and the force calculated from the Maxwell stress tensor.However, the thrust per unit length indicates the performance of the E-sail; thus, F x /L was calculated, as shown in Fig. 2b.F x /L was 11.7, 12.7, and 12.8 nN m −1 for L = 30, 60, and 75 m, respectively.The thrust per unit length for L = 30 m was slightly smaller than that for 60 and 75 m.This lower value of F x /L may be due to the end effect or the effect which arises when L is small in comparison with λ D .At L = 60 and 75 m, the values of F x /L were almost equal; thus, L = 60 m is considered to be sufficiently long to simulate an infinite tether.
The kinks in the thrust between t = 10 and 20 µs shown in Fig. 2 correspond to collections of ambient electrons.Figure 3 shows the time history of the current on the surface of the tether with L = 60 m.The ambient electron current (purple line) varied dramatically between t = 10 and 20 µs.This is an initial response of the ambient electrons to the potential of the tether.The electron plasma frequency ω p was approximately 178 kHz so the response time of ambient electrons were approximately 5.6 µs.Collected electrons do not directly contribute to the thrust, because their momentum is small.Instead, they temporarily shield the potential structure more strongly than in a steady state, causing the thrust to stop increasing with time, as shown in Fig. 2.
Figure 3 also reveals that the dominant source of the force acting on the tether was the force from the Maxwell stress tensor, not from the protons hitting the tether because the proton current was almost zero.This fact shows that the computation successfully simulates the thrust generation by proton deflection.
To compare the present results with previous thrust estimations, F x /L was calculated at various V 0 values, as shown in Fig. 4. The blue line in Fig. 4 is the result of our PIC simulation with L = 60 m and r w = 15 cm.The black and red dotted lines show the thrust estimated by Janhunen and Sandroos (2007) and Sanchez-Torres (2014), respectively, with r w = 15 cm.In the present simulation, F x /L was 67.1 nN m −1 at V 0 = 4.0 kV; in contrast, Janhunen and Sandroos (2007) and Sanchez-Torres (2014) estimated F x /L to be 123 and 246 nN m −1 , respectively.These results indicate that the thrust characteristics of the E-sail are different from those of conventional estimations.This difference is likely due to the difference in the potential structure around the tether.
Figure 5  tained from the present PIC simulations was lower than the potentials obtained using Eqs.( 1) and (3). Figure 5 indicates that potential shielding by ambient electrons is not appropriately included in Eqs. ( 1) and (3). Figure 5 also shows that the sheath length assumed by Eq. ( 3) is consistent with that obtained by PIC simulation.
Figure 6 shows the proton density structure at t = 0.2 ms.At V 0 = 2.0, 3.0, and 4.0 kV, a zero-density (n p = 0 m −3 ) region was present in front of the tether.No protons impinged on the surface of the tether.The momentum of the protons was transferred to the tether through the Coulomb force acting on the tether.In contrast, there was no zero-density region in front of the tether at V 0 = 1.0 kV.This is because the radius r w of the tether is relatively large.The drift velocity of the proton (400 km s −1 ) is equivalent to 0.83keV, but the figure indicates that V 0 = 1.0 kV is not sufficient to deflect all of the protons for r w = 15 cm.Although a wake region is present at V 0 = 1.0 kV, it was formed by protons impinging on the surface of the tether.At a very small r w , a zero-density region may appear for V 0 = 1.0 kV.
This study then considered the high-density proton region in front of the tether.Figure 7 shows the electron and proton density structure along the x axis at V 0 = 2.0 kV.As with Fig. 3a in Janhunen and Sandroos (2007), Fig. 7b reveals a high-density region in front of the tether.The maximum proton density was 2.0 × 10 7 m −3 at V 0 = 2.0 kV.The maximum electron density was approximately 1.0×10 8 m −3 at V 0 = 2.0 kV.The maximum positive charge density of the high-density region is one fifth of the negative charge density of the electrons.The high-density region does not compensate for the potential; that is, it does not reduce the shielding effect of the electrons.Trapped electron removal, which was discussed by Janhunen (2009), was not observed in the present simulation.Thus, the thrust models given by Eqs. ( 1) and (3) are inappropriate for estimating the thrust of the E-sail, and a new model considering appropriate potential shielding must be developed.

Numerical estimation of the thrust of the E-sail
This section describes the proposed method of estimating the thrust of the E-sail and presents the estimation results.First, a semi-analytical solver of the electrostatic potential structure around an infinite cylinder in plasma is introduced.This solver was used to obtain a realistic estimation of the electrostatic potential using a 2-D inverse fast Fourier transform (FFT).The two thrust estimation procedures used to evaluate the thrust characteristics of the E-sail are then described.The first is the effective radius method by Janhunen and Sandroos (2007), and the second is the Coulomb scattering method by Sanchez-Torres (2014).These two methods are hereafter called Methods 1 and 2, respectively.These two evaluation procedures were then used in combination with   the realistic electrostatic potential given by the 2-D FFT.The thrust estimation script written by Python 2.7 is found in the Supplement (available online).

Semi-analytical solver of electrostatic potential structure
The Poisson equation in plasma is expressed as where ρ is the space charge density.From the velocity distribution function of electrons, we assume that the electron's density distribution becomes the Boltzmann distribution: Assuming that the normalized potential eV (r)/k B T e becomes small in the distance and also assuming that n p = n 0 , Eq. ( 4) becomes with the first-order approximation.Defining k D = e 2 n 0 /ε 0 k B T e , we obtain The 1-D solution of Eq. ( 7) is well known and is called the Yukawa potential.However, the 2-D analytical solution of Eq. ( 7), which would represent the shielded potential structure around an infinite tether, remains unknown.This is why previous studies had to assume an artificial potential structure, such as those given by Eqs. ( 1) and (3).Hoshi et al. (2016) developed a numerical method of calculating the potential structure around a tether in plasma, adding the term h 2 k 2 D to the solution of the difference equation corresponding to the differential equation given by Eq. ( 7) as follows space, and h is the cell width used in the Fourier transformation.The 2-D inverse FFT of Eq. ( 8) was taken to obtain the shielded electrostatic potential in plasma, and the solution was found to be consistent with the potential given by the full PIC simulation.
The proposed estimation method was used to obtain a realistic potential without performing the full PIC simulation.
To realize an equipotential within the radius of the tether, the capacity matrix method was also used (Hockney and Eastwood, 1981).Figure 8 compares the potential structures obtained using the 2-D FFT method with the PIC results.A cell width of h = 0.03 and N = 8192 cells in 2-D space were used for the Fourier transformation.The potential structures estimated using the proposed method were consistent with those obtained from the PIC simulation results (Fig. 5).At V 0 = 1 kV, there was a small difference (approximately 20 V) between the potential structures estimated using the proposed method and the PIC simulation results behind the tether, but this difference did not cause a difference in the thrust, because only the potential in front of the tether contributes to the thrust.

Effective radius method (Method 1)
Method 1 considers the effective radius r s that satisfies the following equation: where r s is the distance from the tether at which the electrostatic potential energy is equal to the kinetic energy of the drifting proton.Janhunen and Sandroos (2007) assumed that the scattering cross section is proportional to r s with a coefficient of proportionality of K, meaning the thrust per unit  length dF /dL can be expressed as where P dyn = m p n 0 v 2 d is the dynamic pressure of a solar wind proton.
In Janhunen and Sandroos (2007), V (r) is given by Eq. ( 1).In the present study, V (r) was replaced with the numerical solution obtained using the 2-D FFT method, which considers potential shielding by ambient electrons.For the numerical calculation, the value of r that minimizes the difference between the potential energy and the kinetic energy f 1 , which is expressed as is used as r s instead the value of r.The coefficient of proportionality was set to K = 3.09, as obtained by Janhunen and Sandroos (2007).

Coulomb scattering method (Method 2)
The thrust modeling method by Sanchez-Torres ( 2014) is based on Coulomb scattering.The following equation is formed by adding an angular momentum term to Eq. ( 9): where L m = m p v d ρ is the angular momentum and ρ is the impact parameter.Defining Additionally, l and U eff were defined as U eff (l) is equal to the left-hand side of Eq. ( 13), meaning the distance r min that minimizes f 2 ( ρ r ) = 1 − U eff ( ρ r ) is equivalent to the distance r s obtained using the effective radius method.When the scattering angle is defined as χ (ρ) = π − δ(ρ), as shown in Fig. 4 of Sanchez-Torres (2014), the following equations give χ (ρ): where l max = ρ/r min and l min = ρ/r max .r max denotes the maximum length that scattered particles can be affected by an electrostatic potential.r max was defined as r sh e −b by Sanchez-Torres (2014) and has been calculated as follows by Sanmartín et al. (2008): Then, the thrust per unit length can be expressed as In the present estimation, V (r) in Eq. ( 15) was replaced with the numerical solution obtained using the 2-D FFT method, and the maximum affection length r max = 2.0λ D was used instead of the value of r max = r sh e −b used by Sanchez-Torres (2014).

Estimation results
The thrust per unit length was obtained for V 0 = 0 to 4.0 kV with r w = 15 cm. Figure 9 shows the thrust estimated using the effective radius method (purple line) and the Coulomb scattering method (green line).These modified thrust estimations are similar to the thrust obtained from the PIC simulation, meaning they are also significantly lower than estimations from previous studies.The modified estimated thrusts obtained using Methods 1 and 2 are 79.1 and 62.0 nN m −1 , respectively.The thrust estimated using Method 2 at 4.0 kV is in very good agreement with the PIC result of 67.1 nN m −1 .It should be noted that the difference between the modified and previous estimation methods is simply the electrostatic potential structure.The effective proton deflection area with shielding is smaller than that without the shielding.Thus, the inclusion of potential shielding results in reduced thrust.
These two modified estimation methods yielded similar estimated values of the thrust per unit length.Thus, the two estimation procedures were shown to be essentially similar, and the difference between the original estimated thrust values was revealed to have been caused by differences in the potential structures.The effective radius method contains an approximation coefficient K, whereas the Coulomb scattering method does not have contain an approximation coeffi-   cient; thus, the Coulomb scattering method was considered to be more consistent with the PIC results.

Discussion
In this study, a full PIC simulation of the E-sail was conducted.The most significant difference between the results obtained in this paper and those obtained in previous studies is that the electrostatic potential structure around the tether was considered, which yielded different values of the thrust.
In the PIC simulation by Janhunen and Sandroos (2007), the potential shielding effect of the ambient electrons was not significant.In the analytical estimation by Sanchez-Torres (2014), they assumed the absence of ambient electrons so their estimations of the potential (Eq. 3) were higher than those obtained using Eq.(1) at high positive potentials.In the present PIC simulation, the ambient electrons shield the potential of the tether more effectively than those in Janhunen and Sandroos (2007) and Sanchez-Torres (2014).This study focused on the source of the differences between the present PIC simulation and the PIC simulations performed by Janhunen and Sandroos (2007), which are shown in Table 2.
There are several differences among the simulations, particularly regarding the tether dimensions, plasma conditions, simulation duration, and cell width.We consider that the differences caused the disagreement of the simulation result.Note that the electron temperature k B T e = 100 eV adopted in this paper is the typical value in solar wind at 0.5 AU and is about 1 order of magnitude higher than the typical value at 1 AU (applied in previous studies).
To the best of the authors' knowledge, this study is the first to perform a full 3-D PIC simulation of the E-sail without any approximations.The radius of the tether (r w = 15 cm) was large in comparison with several tens of micrometers; thus, the E-sail was not completely simulated, but the difference is negligible because 15 cm is much smaller than λ D (approximately several meters to a few tens of meters in in-terplanetary space), so the ambient electron collection is not significantly different.
The increase in thrust caused by the removal of electrons, which was discussed by Janhunen (2009), was not investigated in this study; the present simulation did not consider multiple tethers, and the simulation duration (0.2 ms) was not long enough to describe such a effect.If any efficient trapped electron removal mechanisms exist, the thrust of Esail may increase asymptotically, so we must remark that the thrust characteristics obtained in this paper do not consider the effect of the trapped electron removal.However, the presented simulation successfully described the response of ambient electrons as shown in Fig. 3. Hence, our results reveal at least the minimum thrust characteristics of the E-sail.
A modified thrust estimation proposed in this paper, which is obtained by replacing the electrostatic potential structures used in the estimations in the previous studies, is a better reference model of the minimum thrust of the E-sail.The proposed estimation method can be easily used to calculate the minimum thrust of E-sail.

Conclusions
In this study, the first full 3-D PIC simulation of the tether of the E-sail was performed, and the transient of its thrust was numerically calculated.At V 0 ≤ 4.0 kV, the thrust obtained by PIC simulation was almost half of the thrust estimated in previous studies.This difference is caused by the electrostatic potential structure around the tether.The potential structure in the present simulation differed greatly from the structure used in previous estimations due to the strong potential shielding by ambient electrons.
Additionally, a modified thrust estimation method with a shielded potential structure was proposed.In this new method, the potential structure employed in previous estimations was replaced with the potential structure derived from our simulation result.The estimated thrust obtained using the modified method agreed very well with the PIC simulation

Figure 1 .
Figure 1.Definition of the tether model.A tether-like rectangle is located in the center of the computational domain.The solar wind originates from x = 0.

Figure 2 .Figure 3 .
Figure 2. Time history of the thrust with V 0 = 1.3 kV: (a) x component of the total external force; (b) x component of the thrust per unit length.

Figure 4 .Figure 5 .
Figure 4. Comparison of present thrust simulation results with previous estimations.The blue line shows the thrust obtained from the present PIC simulation with L = 60 m and r w = 15 cm.The black and red dashed lines show the thrust estimated byJanhunen and Sandroos (2007) and Sanchez-Torres (2014), respectively, with r w = 15 cm.The estimation procedures used in these two previous studies are described in Sect.3.

Figure 8 .
Figure 8.Comparison of potential structures at surface potentials of (a) V 0 = 1 kV, (b) V 0 = 2 kV, (c) V 0 = 3 kV, and (d) V 0 = 4 kV.Blue lines show the potential structure obtained using the proposed method.Purple lines show the results of the PIC simulations.

Figure 9 .
Figure 9.Comparison of thrust estimations.The purple line shows the modified estimation obtained using the effective radius method.The green line shows the modified estimation based on the Coulomb scattering method.The blue and black lines are the same as those in Fig. 4.