Particle Simulation of Plasma Drag Force Generation in the Magnetic Plasma Deorbit
Abstract
The magnetic plasma deorbit method has been proposed for micro and nanosatellites in the low earth orbit (LEO). The magnetic plasma deorbit utilizes the plasma drag force generated by the interaction between space plasma and the magnetic field surrounding magnetic torquers (MTQs). In this study, the plasma drag force generated by the magnetic plasma deorbit method is estimated by using a plasma flow simulation based on the fully kinetic model. The dependences of the plasma drag force on the MTQ angle, atmospheric plasma density, and MTQ magnetic moment are investigated. The simulation results show that the structures of the bow shock and magnetosphere are formed in the vicinity of the MTQs. The predicted plasma drag force is 105.2 nN for an MTQ of 10 A m when the MTQ magnetic moment is parallel to the satellite velocity direction. The plasma drag force shows a strong dependence on the magnetic moment angle. In addition, the plasma drag force is proportional to the atmospheric plasma density and the MTQ magnetic moment. Steady drag force generation in the magnetic plasma deorbit is validated, and the sensitivities of the drag force to the relevant parameters are elucidated.
1
Particle Simulation of Plasma Drag Force Generation in the Magnetic Plasma Deorbit
Rei Kawashima and Junhwi Bak 
The University of Tokyo, Tokyo 1138656, Japan 
and 
Shinji Matsuzawa and Takaya Inamori 
Nagoya University, Aichi 4648603, Japan 
JSR˙kawashima.nls
I Introduction
Because the development of nano and microsatellites is fast and lowcost, a number of them were successfully launched for various types of missions. Most of these satellites are used in the low earth orbit (LEO). In general, these satellites complete their missions within ten years; furthermore it must be noted that many of them do not stay in orbit for very long periods of time. The lifetime of CubeSatclass satellites is less than 200 days on average [1]. There also exist longlifetime nanosatellites. An example of a longlifetime nanosatellite is CubeSat XIV [2]. This CubeSat was launched in October 2005, and it is still operational in December 2017. One of the reasons for the reduced lifetime of nanosatellites is inadequate thermal or electronic designs. If the technologies of nanosatellite development advance further, the average lifetime of nanosatellites is expected to increase in the future. However, the increasing development and use of nanosatellites will inevitably cause the issue of space debris. Propulsion systems are generally not installed on nanosatellites because of strict constraints on mass and volume. Hence, deorbit methods for nano and microsatellites have been of concern in recent years to prevent the accumulation of space debris.
As a deorbit method for nano and microsatellites, a method referred to as the magnetic plasma deorbit has been proposed by the authors [3]. The concept of the method is illustrated in Fig. 1. This method uses the plasma drag force generated by the interaction between space plasma and the magnetic field surrounding a satellite. The magnetic field is generated by a magnetic torquer (MTQ), which is generally installed on nanosatellites for attitude control. Thus, the advantage of this method is that it requires no additional equipment such as propulsion systems or large structures for deorbit.
The fundamental physics of the drag force generation in the magnetic plasma deorbit is supposed to be similar to that of magnetic sails [4]. In the research of magnetic sails, Funaki et al. proposed an analytical model for the prediction of the plasma drag force as follows [5]:
(1) 
where is the characteristic length of the magnetosphere, and it is estimated by considering the balance of magnetic and plasma pressures as follows:
(2) 
In Ref. [3], to evaluate the plasma drag force generated in the magnetic plasma deorbit method, Eqs. (1) and (2) are used. However, this analytic model was used for magnetic sails where solar wind plasma was involved, and the model was not validated for the magnetic plasma deorbit in LEO plasma. It is noted that the concept of the magnetic plasma deorbit in Fig. 1 is made based on the plasma flow characteristics of magnetic sails. Thus, it is also unknown whether various structures, such as the bow shock and magnetosphere, are formed in the case of the magnetic plasma deorbit.
Similarly, the sensitivities of the plasma drag force in the magnetic plasma deorbit to the relevant parameters were not elucidated. First, the effect of the magnetic moment angle on the plasma drag force is concerned. In the case of magnetic sails, it is known that the has an influence on the plasma drag force [6], and hence is supposed to affect the performance of the magnetic plasma deorbit. Moreover, the effect of the atmospheric plasma density and magnetic moment of an MTQ is concerned, since these parameters vary in wide ranges in the cases of nanosatellites in the LEO. The aforementioned model yields the relations of and . The validity of these relations must be checked in the case of the magnetic plasma deorbit.
In this study, a plasma flow simulation based on the fully kinetic model is conducted to analyze the plasma drag force generated in the magnetic plasma deorbit. The primary objective is to confirm steady plasma drag force generation in the magnetic plasma deorbit method. Further, the simulated plasma drag force is compared with that of the analytic model in Eqs. (1) and (2), and the validity of the model is examined. Additionally, a parametric study is performed to investigate the effect of , , and on the plasma drag force. Concerning the relationships between the plasma drag force and these parameters, the simulation results are compared with those of the analytic model.
Ii Basics of the Electrostatic Fully Kinetic Model
ii.a Plasma Flow
The fully kinetic particle model is used for the description of the magnetized plasma flow. The motions of ions and electrons are computed by using the twodimensional (2D) threevelocity particleincell (PIC) method. Several millions of particles are bundled into one macroparticle, and each macroparticle is accelerated by outer forces. The equation of motion for each macroparticle is formulated as follows:
(3) 
(4) 
is the electromagnetic (EM) force and subscript s denotes each species. It is noted that is the relative speed to the satellite as depicted in Fig. 1. The EM force acting on each macroparticle can be expressed as follows:
(5) 
where is equal to for ions and for electrons. The magnetic field is calculated by the superposition of the MTQ magnetic field and induced magnetic field: .
An appropriate plasma flow model was discussed in previous researches on magnetic sails [7, 8]. The criterion for determining the plasma flow model is the ratio of the Larmor radius to the magnetosphere size . In the case of LEO plasma and an MTQ strength of 1–10 A m, the ratio of the ion Larmor radius to is supposed to be on the order of 10, and the ratio of the electron Larmor radius to is on the order of 1. In this case, the kinetic model should be used for both ions and electrons. Thus, the fully kinetic model is used in the modeling of the LEO plasma flow surrounding the nanosatellite MTQs. Interparticle collisions are neglected in the present model. All of the meanfree paths of ion–ion Coulomb, electron–ion Coulomb, ion–neutral momentumexchange, and electron–neutral momentumexchange collisions are estimated to be greater than 10 m under atmospheric conditions at high altitudes where the magnetic plasma deorbit method is used. Very few particles experience collisions at the spatial scale of the magnetic plasma deorbit since it is 1 m as shown later.
ii.b EM Field
In this study, the electric field is assumed electrostatic, and the space potential is obtained through Gauss’s law. In the 2D simulation in the plane, the uniformity of space potential is assumed in the direction. The 2D Poisson equation of Gauss’s law is formulated as follows:
(6) 
After deriving the space potential, the electric field is calculated by the relation: .
The magnetic field is also assumed static, and it is calculated through the vector potentials of the MTQ and induced magnetic field. Concerning the MTQ, a line dipole is assumed for the 2D simulation in the plane. The vector potential of a line dipole is expressed as follows:
(7) 
where the definition of the coordinate and are as illustrated in Fig. 2. The vector potential for the induced magnetic field is calculated by using Poisson’s equation for a vector potential as follows:
(8) 
Note that the current in the direction is induced by the effects of magnetization. As discussed in Sec. IV.B., ions are not magnetized with the magnetic field strength concerned in the magnetic plasma deorbit method, and hence the induced magnetic field is mainly caused by the electron current. After deriving the magnetic vector potentials, the magnetic field in the plane is calculated as follows:
(9)  
(10) 
In the present model, the assumed electrostatic and magnetostatic fields that do not include the effects of displacement current. Here, we call this model the electrostatic–magnetostatic (ESMS) model. Information about the displacement current propagates in space with the speed of light, which means almost instantaneously at the time scale of the ion flow. Note that the ratio of the speed of light to the speed of the ion flow surrounding a satellite in the LEO (i.e., condition number) is . The time derivative of the electric field would be small at the ion time scale, and the ESMS model is supposed to be a good approximation for the simulation of the magnetic plasma deorbit method. Another approach to obtaining the electric and magnetic fields is to directly solve Faraday’s law and Ampere’s law in Maxwell’s equations [9]. The model is called the EM model, and it was used in the full particle and particlefluid hybrid modeling of magnetic sails [8, 10]. It was reported that the ESMS model yields almost the same thrust as the EM model in the full particle simulations of a 2D magnetic sail, for which the condition number is 600 [11]. This fact supports the validity of the ESMS model for steadystate simulations of the magnetic plasma deorbit method.
Further, the Earth geomagnetic field (GMF) is neglected in this study. In the research of magnetic sails, Nishida and Funaki conducted a numerical simulation of a magnetic sail with the interplanetary magnetic field (IMF) and observed that the IMF affects the thrust characteristics [6]. It was reported that the IMF enhances the thrust and induces a pitching moment on the magnetic sail. The detailed effect of the background magnetic field is complicated even in the sense of qualitative trends since two angles of the satellite magnetic moment and IMF are involved. The GMF is also supposed to influence the performances of the magnetic plasma deorbit method. In the case of the magnetic plasma deorbit method used in the LEO, the direction of the GMF changes depending on the position of satellites, i.e., longitude and latitude. Hence, a detailed analysis with various combinations of the MTQ and GMF fields is required to understand the GMF effect comprehensively. In this paper, we focus on the case of no GMF for simplicity to benchmark the basic characteristics of the magnetic plasma deorbit method.
ii.c Artificial Electron Mass Model
In particle simulations of plasma, there are mainly three restrictions for the time step: 1) resolution of the plasma frequency, 2) Courant–Friedrichs–Lewy (CFL) condition for electron velocity, and 3) CFL condition for ion velocity. Here, the CFL conditions are related to the speed of information propagation in electrostatic electron or ion waves [12]. Among the restrictions above, the CFL condition for electron velocity typically requires the smallest time step in the case of the LEO plasma flow. The use of a small time step results in a large number of time steps to achieve the steady state at the time scale of ions. This issue of numerical stiffness arising from a large mass ratio of appears in many full PIC simulations, and concessions must be made to speed up the convergence [13, 14]. In the present modeling, to accelerate the simulation, the artificial electron mass (AEM) model is used. In the AEM model, the electron mass is artificially increased as follows:
(11) 
Increasing the electron mass by a factor reduces the electron velocity by and accelerates the convergence by .
The magnetic plasma deorbit uses charge separation caused by the difference in the Larmor radii between ions and electrons. Thus, the electron gyroradius must be maintained while is introduced. The electron gyroradius is expressed as follows:
(12) 
In this simulation, the magnetic flux density for electrons is increased as follows:
(13) 
The AEM factor is selected as for all the simulation cases in this study. The validity of the selection of this value is discussed in the appendix. After the application of the AEM model, the time step is typically set to s.
Iii Calculation Condition and Numerical Method
iii.a Calculation Condition
In the present study, satellites at an altitude of 300–1000 km are concerned. The plasma environment in this altitude range is referred to as the exosphere [15]. This region is inside the magnetopause generated by the geomagnetic field, and the effect of highvelocity solar wind is small. Table 1 lists the assumed parameters for the simulation of the magnetic plasma deorbit in the LEO. The properties of plasma are calculated through the International Reference Ionosphere (IRI) 2012 model as functions of the altitude and time [16]. In this simulation, the altitude of km and daytime are assumed as the calculation conditions. All ions are assumed to be O. This is because the fraction of O is greater than 0.9 for both daytime and nighttime, and the contributions of other ion species to the drag force are supposed to be small. The ion and electron temperatures are also obtained from the IRI model. The satellite velocity is 8,000 m s. In the simulation, the ion flow velocity is assumed to be equal to the satellite velocity in the control volume surrounding a satellite. In the exosphere region, the plasma wind tends to blow horizontally from the subsolar heated side to the night cold side. The typical wind speed ranges from 100 m/s to 300 m/s for quiet geomagnetic conditions [15]. Since the wind speed is smaller than the satellite velocity by one order of magnitude, the effect of plasma wind is neglected in the present study. The calculation domain is illustrated in Fig. 2. The control volume surrounding the satellite is used, and the plasma flow is introduced to the domain from the lefthand side boundary. It is known that the typical satellite velocity is larger than the ion thermal velocity in the LEO environments [17]. Hence, the number flux density of ions from the lefthand side boundary is given by using the satellite velocity as . However, for the top and bottom boundaries, the number flux density of ions is given by a random flux. The random flux can be written as follows:
(14) 
where denotes each species. Concerning electrons, the thermal velocity is much greater than the satellite velocity. Thus, the number flux density of electrons based on the random flux is introduced to the calculation domain from the left, right, top, and bottom boundaries. Because the satellite velocity is much smaller than the speed of light, the Lorentz transformation is not considered in this calculation. A satellite is located at the center of the calculation field, and the dipole magnetic field is generated around the satellite.



Parameter  Symbol  Value 
Ion number density, m  1.0 10  
Ion species  –  O 
Ion temperature, K  1230  
Electron temperature, K  2000  
Satellite velocity, m s  8000  
Magnetic moment, A m  10  
Angle of MTQ  0–  

iii.b Numerical Method
The MTQ magnetic field is input at the beginning of the simulation. In addition, one million macroparticles of ions and electrons are uniformly distributed in the calculation region. The basic computational flow in the time marching consists of the electrostatic and magnetic field calculations and particle motion calculations. Possion’s equation in Eq. (6) is solved for the scalar potential by applying the secondorder central differencing. An iterative method of successive line overrelaxation method is used with a truncation error of 10. After deriving the scalar potential, the electric field is calculated by applying the secondorder central differencing. The vector potential and the magnetic field are also calculated by applying the secondorder central differencing to Eqs. (8) and (10).
The overall time advancement is implemented in an explicit manner. The time advancement of particle motions described by Eqs. (3) and (4) is computed by using the Buneman–Boris method [18]. After the motion calculations of ions and electrons, the charge density and the current density are calculated by weighting the particle mass and momentum to each cell. The weighting function in the PIC method is the secondorder spline function. The number of macroparticles for each species in the calculation field is normally one million, which corresponds to 25 particles per cell. The inflow of particles from the boundaries is implemented by assuming the Maxwellian particle energy distribution with the ion and electron temperatures listed in Table 1. The plasma drag force is computed by integrating the momentum changes of particles in the xdirection, while they traverse the calculation field.
Iv Results and Discussion
iv.a Confirmation of Steady Drag Force Generation
After several millions of time steps, the steady state is confirmed in the plasma flow. A typical iteration history of the calculated plasma drag force is shown in Fig. 3. Here, the case of = 10 m, = 10 A m, and = 0 is simulated. In this case, the quasisteady state is reached at ms with the tolerable statistical noise stemming from the particle model. The steady plasma drag force is estimated by averaging the history for the last 0.5 ms. In what follows, this timeaveraged steady drag force is discussed in each case.
The 2D distributions of the ion number density and space potential are shown in Fig. 4. The bow shock and magnetosphere are observed as shown in Fig. 4(a). Here, the boundaries of the bow shock and magnetosphere are defined by the contour lines of and , respectively. The coefficient of 0.5 is chosen for clear visualization of the boundary. As shown in Sec. IV.B., the ion number density steeply changes at the boundary of the magnetosphere. Thus, even if this coefficient is changed between 0.2 and 0.8, the location of the boundary is almost unchanged. As shown in Fig. 4(b), the highpotential region is formed at the upstream of the MTQ center. The highpotential region is formed by the charge separation between the ion and electron flows, and this charge separation is generated owing to the magnetization of electrons by the MTQ magnetic field. The highpotential region brakes and deflects ions approaching the MTQ center. The structures of the bow shock and magnetosphere are formed as consequences of the ion braking. As a conclusion, the drag force generation mechanism as shown in Fig. 1 is confirmed by the simulation.
In the case of the result in Fig. 3, the simulated steady plasma drag force is 105.2 nN. If one substitutes this plasma drag force into the analytical model in Eq. (1), the corresponding drag coefficient is calculated as . In the experiments and simulations of a magnetic sail, the drag coefficient is estimated as 0.5–2.0 [5, 20]. The very small in the magnetic plasma deorbit is attributed to the small magnetosphere size . Ashida et al. showed that is no longer constant against in the cases of small magnetospheres, in which the ratio of the ion Larmor radius to the magnetosphere size, , is greater than unity [8]. Under the calculation condition of the magnetic plasma deorbit, the ratio of is approximately 10. Hence, the analytic model of magnetic sails is not valid for the prediction of the plasma drag force in the magnetic plasma deorbit method. The validity of the simulated plasma drag forces is discussed in Appendix B via a numerical simulation of a smallscale magnetic sail.
iv.b Dependence of the Plasma Drag Force on the Magnetic Moment Angle of an MTQ
Figure 5 shows the variation of the plasma drag force when the MTQ magnetic moment angle is changed from 0 to with an increment of . The maximum is observed when , i.e., the MTQ magnetic moment is parallel to the direction of the satellite velocity. decreases as increases, and the minimum is generated when , i.e., the MTQ magnetic moment is perpendicular to the direction of the satellite velocity. In what follows, the plasma drag forces in the cases of and are referred to as and , respectively. The similar trend of was also confirmed in the magnetohydrodynamic analysis of a magnetic sail [6]. Therefore, the dependence of on is qualitatively consistent with that of the magnetic sail.
However, in the quantitative sense, the dependence of on does not agree with that of the magnetic sail. In the case of the magnetic plasma deorbit, and are 105.2 nN and 24.2 nN, respectively, and the ratio of these drag forces is about 3.5. However, this ratio is approximately 1.5 in the case of a magnetic sail [20]. The plasma drag force generated in the magnetic plasma deorbit shows a strong dependence on the magnetic moment angle, compared with the magnetic sails. One of the possible reasons is the assumption of the dipole magnetic field in this simulation. In the magnetic sail simulation in Ref. [20], the initial magnetic field was generated by a ring coil. If the magnetic field is generated by ring coils, finitecoil radius effects appear and may change the dependence of the drag force on the magnetic moment angle.
Another possible cause of the strong dependence on the magnetic moment angle is considered using the 2D distributions. The steadystate distributions of the ion number density and space potential in the case of are shown in Fig. 6. In the analytical model in Eq. (1), the plasma drag force is associated with the magnetosphere size . However, is obtained by Eq. (1) as m in the case of A m, which is much larger than the magnetosphere sizes observed in this simulation. Hence, another definition is used for the magnetosphere size to reflect the simulation results. In Figs. 4 and 6, the boundary of the magnetosphere is defined by the contour lines of m. Based on the simulation results, is defined as the magnetosphere size in the direction at m. is estimated as 0.44 m and 0.40 m in the cases of of 0 and , respectively. The difference in is not so large as the difference between and . Thus, the difference in the magnetosphere size cannot be the reason for the strong dependence on the magnetic moment angle of the magnetic plasma deorbit.
In the case of , the ion reflection from the center of the MTQ to the direction is observed, as shown in Fig. 4(a). The reflected ions are not observed in the case of . The reflected ions may contribute to the strong plasma drag force of . The drag force generated by the reflected ions can be roughly estimated by using the distribution. The critical space potential , which is required for the reflection of incident ions, can be calculated as follows:
(15) 
In the case of = 8.0 km s, is calculated as 5.3 V.
The region where is visualized in Fig. 4(b). Here, this region is called the highpotential region. Incident ions cannot penetrate the highpotential region, and they are reflected to the direction. Note that the highdensity region does not exist in the result of the case. In the case of , the size of the highpotential region in the direction, , is estimated as m. The ion number flux flowing in the direction and approaching the highpotential region is given as
(16) 
where is the ion number flux entering the highpotential region. Note that the unit length is considered in the direction. If one assumes that ions approaching the highpotential region are reflected to the direction horizontally, the velocity of these ions changes from to . Thus, the drag force generated by the reflected ions can be calculated from the momentum changes of these ions as follows:
(17) 
where is the estimated drag force generated by the reflected ions. is calculated as 100.3 nN when m. This value agrees well with the plasma drag force in the case of . Therefore, the cause of the strong dependence of the plasma drag force on the magnetic moment angle is attributed to the strong drag force by the reflected ions as the case of .
The electric field also plays an important role in the drag force generation in the case of , even if there is no ion reflection. This fact can be derived by comparing the bow shock thickness with the Larmor radii of ions and electrons. Figure 7(a) shows the onedimensional distribution of the ion number density along the line of m in Fig. 6. of the incoming flow is given by , and it is discontinuously increased in the bow shock region. The plasma cavity exists behind the shock wave, and almost all the ions are evacuated around the MTQ dipole center. The peak of is observed at m. If the boundary of the bow shock is defined by the contour line of , the bow shock thickness is estimated as m. The distributions of the ion Larmor radius and electron Larmor radius are plotted in Fig. 7(b). Note that the magnetic field generated by the MTQ dipole is considered in this analysis. and at the location of the peak are 43.3 m and 0.05 m, respectively. is close to the bow shock thickness , and hence electrons are magnetized at the location of the shock wave. On the other hand, is greater than by two orders of magnitude. This means that ions are unmagnetized and mainly accelerated by electrostatic forces.
iv.c Parametric Study on the Atmospheric Plasma Density and MTQ Magnetic Moment
A parametric study was conducted to investigate the dependence of the plasma drag force on the atmospheric plasma density and magnetic moment of the MTQ. The dependence of the simulated plasma drag force on the atmospheric plasma density is shown in Fig. 8. In this analysis, the atmospheric plasma density in the orbit is varied between 0.3–3.0 m. This plasma density range corresponds to the altitude ranges of 550–1200 km during daytime and 250–650 km during nighttime. Additionally, the dependence of the plasma drag force on the MTQ magnetic moment is shown in Fig. 9. The magnetic moment is varied in the range of 3–30 A m. First, the relation of is always confirmed in this parametric study. In the researches on magnetic sails, the relation of is also observed in the cases of large , for which the ratios of the Larmor radius to the magnetosphere size for ions and electrons become and , respectively [7]. However, of 30 A m is already large for the MTQs in nanosatellites, and it would be difficult to use a sufficiently large to make the conditions of and . Therefore, the relation of is supposed to be valid within the range of used in the magnetic plasma deorbit method. Concerning the , the relation of is confirmed, and is also observed when A m. According to the analytic model in Eqs. (1) and (2), the plasma drag force is supposed to change according to the relations of and . Thus, the analytic model is not valid in predicting the effect of variations in and on the plasma drag force in the magnetic plasma deorbit. Note that Ashida et al. also showed that is proportional to in the cases of smallscale magnetic sails where finiteLarmor radius effects become significant [21].
iv.d Deorbit Analysis
The deorbit duration achieved by the magnetic plasma deorbit method is estimated by using the plasma drag forces obtained by the present simulation. First, the effect of the aerodynamic force is considered. The aerodynamic drag force is generally formulated as follows:
(18) 
The atmospheric gas density at each altitude can be calculated by using the NRLMSISE00 atmosphere model [22]. Numerous models have been proposed to model the drag coefficient [23, 24, 25]. It is known that a satellite of a relatively simple shape like a sphere has a drag coefficient ranging 2.0–2.4 in the LEO [26]. Here, a constant drag coefficient of is assumed for simplicity. Concerning the plasma drag force, the case of the magnetic moment parallel to the satellite velocity is considered, i.e., is assumed as the plasma drag force. Figure 10 shows the dependences of the plasma and aerodynamic drag forces on the altitude. with face areas of m and m are plotted, whereas with MTQ magnetic moments of A m and A m are shown. Both of and are weakened as the altitude becomes higher, because both the atmospheric gas and plasma densities decrease exponentially with the altitude. At a relatively low altitude of 500 km, the aerodynamic force is predominant. However, decreases more rapidly with the altitude than , and becomes comparable with within 500 km 600 km. At an altitude exceeding 600 km, the plasma drag force is more effective than the aerodynamic force. Therefore, the magnetic plasma deorbit method is expected to be beneficial for the satellite deorbit from high altitudes of 600 km.
An analytic model to predict the operation period of the magnetic plasma deorbit method is proposed in Ref. [3]. In this model, the plasma density dependency on the altitude is approximated as follows:
(19) 
where is a constant fitting parameter, which is assumed to be m [3]. The satellite semimajor axis is expressed by a perturbation equation for a circular orbit. The deorbit duration is obtained by analytically integrating the perturbation equation as follows:
(20) 
where
(21) 
and are the standard gravitational parameter and Earth’s radius, respectively. The detailed derivation process of the deorbit duration is explained in Ref. [3]. Only the plasma drag force is taken into consideration in this model, and Eq. (20) gives an instantaneous estimation of the deorbit duration. By using Eq. (1), Eq. (20) can be rewritten with the plasma drag force . Here, is assumed as the plasma drag force, and is further assumed according to the discussion in Sec. IV.C. Owing to these assumptions, the plasma drag force is expressed as , and Eq. (20) is rewritten as follows:
(22) 
Note that in Eq. (1) is replaced by in the case of the magnetic plasma deorbit. The reference values are determined based on the condition of km, as follows: km, A m, and nN. A constant satellite velocity is assumed as m s. The deorbit duration is estimated for 1 kg and 3 kg satellites with the initial altitude of 800 km. The results are plotted with various MTQ magnetic moments in Fig. 11. For instance, of 30 A m is required for the deorbit of a 3 kg satellite within 25 y. The effectiveness of the magnetic plasma deorbit method is quantitatively evaluated by this analysis. In order to predict the deorbit duration more accurately, one needs to consider several effects, such as the satellite attitude, geomagnetic field, finite satellite body and coil, hour of day, and so on. In addition, the optimization of the deorbit scenario considering these effects would also be one of the themes of subsequent papers.
V Conclusion
A plasma flow simulation was conducted to predict plasma drag forces generated by the magnetic plasma deorbit method. Static electric and magnetic fields are assumed, and the plasma flow is described by using the fully kinetic model. The simulation is continued until a steady state is reached, and the steady plasma drag force is calculated. The simulated plasma drag forces are compared with those predicted by the analytic model of magnetic sails, and the validity of the model is examined. The findings are summarized as follows:

The structures of the bow shock and magnetosphere are observed in the vicinity of the MTQ, and a steady plasma drag force is generated. The basic concept of the drag force generation is confirmed.

The predicted is 105.2 nN in the case of m, km s, A m, and . is estimated as 0.2. This small indicates that the analytic model of magnetic sails is not applicable to the prediction of plasma drag forces in the magnetic plasma deorbit.

The plasma drag force shows a strong dependence on the MTQ magnetic moment angle. The drag force in the case of the magnetic moment parallel to the satellite velocity (i.e., ) is 3.5 times larger than that in the case of the magnetic moment perpendicular to the satellite velocity (i.e., ). This strong dependence is attributed to a strong drag force in the case of , which is mainly generated by reflected ions. The relation of is valid within the conditions of the magnetic plasma deorbit method.

The dependence of the plasma drag force in the case of the magnetic moment parallel to the satellite velocity () on the atmospheric plasma density and the MTQ magnetic moment was linear.
Appendix A: Validity of the AEM Model
As stated in Sec. II.C, the AEM model is used in the simulation to mitigate the issue of numerical stiffness. The validity of the AEM model is examined by confirming the convergence of the plasma drag force when the AEM coefficient is changed from 150 to 500. ( case) and ( case) calculated with various are shown in Fig. 12. In this figure, and are normalized by 100 nN and 25 nN, respectively. Note that means that the real electron mass is used and the electron motion approaches the real situation. In the case of , and are 105.1 nN and 24.6 nN, respectively. As is increased, the calculation time is reduced, but the difference between the AEM and the real electron mass becomes large. In the case of , and are 102.4 nN and 24.1 nN, respectively. The differences in the plasma drag forces between the cases of and are and for and , respectively. Therefore, the variations in and are considered to be small when . It is important to find a point of compromise that satisfies both reasonable calculation time and accuracy. In the simulation of the magnetic plasma deorbit, is selected based on several test calculations.
Appendix B: Model Validation by a Magnetic Sail Simulation
To check the validity of the present model, the thrust of a magnetic sail is simulated, and the result is compared with previous numerical simulation results. A magnetic sail is one of the propulsion systems for interplanetary transportation originally proposed by Zubrin and Andrews [4]. This system obtains its thrust by the interaction between solar wind plasma and the magnetic field generated by an onboard magnetic coil. Although the plasma properties of the solar wind are different from those of the LEO plasma, a numerical simulation of a magnetic sail would yield a quantitative evaluation regarding the adequacy of the present model.
As written in Sec. II.A., the plasma flow surrounding a dipole magnetic field can be characterized by the ratio of the Larmor radius to the magnetosphere size . In the case of the magnetic plasma deorbit, the criteria for ions and electrons are and , respectively. The thrust of magnetic sails under the conditions of and was evaluated by threedimensional full particle simulations in Ref. [21]. The properties of the solar wind assumed in this simulation are as follows: plasma density is 510 m, solar wind velocity is 510 m s, plasma temperature is 10 eV, and all ions are assumed to be H. The interplanetary magnetic field is ignored, and the dipole magnetic moment is assumed.
A magnetic sail with a dipole of A m is simulated by using the present model in the calculation domain of 800 m800 m with a grid of 300300. The dipole is located at the center of the domain and with a magnetic moment angle of . The 2D distribution of the ion number density is shown in Fig. 13. The solar wind is deflected by the magnetic field, and a plasma cavity is formed behind the dipole. The simulated thrust at the steady state is given as nN. An empirical equation to predict the thrust of a smallscale magnetic sail is presented based on the results of the threedimensional simulations as follows [21]:
(23) 
In the case of A m, the thrust predicted by this empirical equation is nN. The difference between and is 14.2 nN, which is 24% of . Thus, the simulation using the present model may include this level of uncertainty. Note that this comparison is made between the present 2D simulation and the previous threedimensional simulations. A more suitable validation method is to compare simulation results with experimental data. This is one of the themes of future works.
References
 [1] Swartwout, M., “The First One Hundred CubeSats: A Statistical Look,” Journal of Small Satellites, Vol. 2, No. 2, 2013, pp. 213–233.
 [2] Funase, R., Takei, E., Nakamura, Y., Nagai, M., Enokuchi, A., Yuliang, C., Nakada, K., Nojiri, Y., Sasaki, F., Funane, T., Eishima, T., and Nakasuka, S., “Technology demonstration on University of Tokyo’s picosatellite “XIV”and its effective operation result using ground station network,” Acta Astronautica, Vol. 61, No. 7, 2007, pp. 707–711.
 [3] Inamori, T., Kawashima, R., Saisutjarit, P., Sako, N., and Ohsaki, H., “Magnetic plasma deorbit system for nano and microsatellites using magnetic torquer interference with space plasma in low Earth orbit,” Acta Astronautica, Vol. 112, 2015, pp. 192–199.
 [4] Zubrin, R. M. and Andrews, D. G., “Magnetic sails and interplanetary travel,” Journal of Spacecraft and Rockets, Vol. 28, No. 2, 1991, pp. 197–203.
 [5] Funaki, I., Kojima, H., Yamakawa, H., Nakayama, Y., and Shimizu, Y., “Laboratory Experiment of Plasma Flow Around Magnetic Sail,” Astrophysics and Space Science, Vol. 307, No. 1, 2007, pp. 63–68.
 [6] Nishida, H. and Funaki, I., “Analysis of Thrust Characteristics of a Magnetic Sail in a Magnetized Solar Wind,” Journal of Propulsion and Power, Vol. 28, No. 3, 2012, pp. 636–641.
 [7] Kajimura, Y., Funaki, I., Matsumoto, M., Shinohara, I., Usui, H., and Yamakawa, H., “Thrust and Attitude Evaluation of Magnetic Sail by ThreeDimensional Hybrid ParticleinCell Code,” Journal of Propulsion and Power, Vol. 28, No. 3, 2012, pp. 652–663.
 [8] Ashida, Y., Funaki, I., Yamakawa, H., Usui, H., Kajimura, Y., and Kojima, H., “TwoDimensional ParticleInCell Simulation of Magnetic Sails,” Journal of Propulsion and Power, Vol. 30, No. 1, 2013, pp. 233–245.
 [9] Harned, D. S., “Quasineutral hybrid simulation of macroscopic plasma phenomena,” Journal of Computational Physics, Vol. 47, No. 3, 1982, pp. 452–462.
 [10] Kajimura, Y., Usui, H., Funaki, I., Ueno, K., Nunami, M., Shinohara, I., Nakamura, M., and Yamakawa, H., “Hybrid ParticleinCell Simulations of Magnetic Sail in Laboratory Experiment,” Journal of Propulsion and Power, Vol. 26, No. 1, 2017/12/20 2010, pp. 159–166.
 [11] Ashida, Y., “Study on Propulsive Characteristics of Magnetic Sail and Magneto Plasma Sail by Plasma Particle Simulations,” Ph.D. Thesis, Kyoto University, 2014.
 [12] Chen, F. F., Introduction to plasma physics and controlled fusion, Plenum Press, New York, 2nd ed., 1984.
 [13] Cho, S., Komurasaki, K., and Arakawa, Y., “Kinetic particle simulation of discharge and wall erosion of a Hall thruster,” Physics of Plasmas, Vol. 20, No. 6, 2013, pp. 063501.
 [14] Szabo, J., Warner, N., MartinezSanchez, M., and Batishchev, O., ‘‘Full ParticleInCell Simulation Methodology for Axisymmetric Hall Effect Thrusters,” Journal of Propulsion and Power, Vol. 30, No. 1, 2013, pp. 197–208.
 [15] Schunk, R. and Nagy, A., Ionospheres: physics, plasma physics, and chemistry, Cambridge University Press, 2009.
 [16] Bilitza, D., “International Reference Ionosphere 1990,” NASA Technical Report, 1990, pp. 0–84.
 [17] Hastings, D. E., “A review of plasma interactions with spacecraft in low Earth orbit,” Journal of Geophysical Research: Space Physics, Vol. 100, No. A8, 1995, pp. 14457–14483.
 [18] Birdsall, C. K. and Langdon, A. B., Plasma physics via computer simulation, CRC Press, 2004.
 [19] Mason, R. J., “Implicit moment particle simulation of plasmas,” Journal of Computational Physics, Vol. 41, No. 2, 1981, pp. 233–244.
 [20] Nishida, H., Ogawa, H., Funaki, I., Fujita, K., Yamakawa, H., and Nakayama, Y., “TwoDimensional Magnetohydrodynamic Simulation of a Magnetic Sail,” Journal of Spacecraft and Rockets, Vol. 43, No. 3, 2006, pp. 667–672.
 [21] Ashida, Y., Yamakawa, H., Funaki, I., Usui, H., Kajimura, Y., and Kojima, H., “Thrust Evaluation of SmallScale Magnetic Sail Spacecraft by ThreeDimensional ParticleinCell Simulation,” Journal of Propulsion and Power, Vol. 30, No. 1, 2014, pp. 186–196.
 [22] Picone, J. M., Hedin, A. E., Drob, D. P., and Aikin, A. C., “NRLMSISE00 empirical model of the atmosphere: Statistical comparisons and scientific issues,” Journal of Geophysical Research: Space Physics, Vol. 107, No. A12, 2002, pp. SIA 15–1–SIA 15–16, 1468.
 [23] Moe, M. M., Wallace, S. D., and Moe, K., “Refinements in determining satellite drag coefficients  Method for resolving density discrepancies,” Journal of Guidance, Control, and Dynamics, Vol. 16, No. 3, 1993, pp. 441–445.
 [24] Mostaza Prieto, D., Graziano, B. P., and Roberts, P. C. E., ‘‘Spacecraft drag modelling,” Progress in Aerospace Sciences, Vol. 64, No. Supplement C, 2014, pp. 56–65.
 [25] Vallado, D. A. and Finkleman, D., “A critical assessment of satellite drag and atmospheric density modeling,” Acta Astronautica, Vol. 95, No. Supplement C, 2014, pp. 141–165.
 [26] Moe, K., Moe, M. M., and Wallace, S. D., “Improved Satellite Drag Coefficient Calculations from Orbital Measurements of Energy Accommodation,” Journal of Spacecraft and Rockets, Vol. 35, No. 3, 1998, pp. 266–272.