Numerical investigation of internal plasma currents
in a magnetic nozzle

Electrodeless plasma thrusters have been developed worldwide and are expected to have a long lifetime.[1] They mainly consist of a radiofrequency (rf) antenna, a dielectric tube, and a solenoid. The solenoid generates both an axial magnetic field and a magnetic nozzle downstream of the dielectric tube, which transports the plasma axially along the field lines and accelerates it.[2, 3, 4, 5, 6, 7] The physics underlying the acceleration of the magnetic nozzle has also been investigated in astrophysics,[8] and the performance of electrodeless plasma thrusters has been improved as an engineering application. In a recent experiment, a thrust of 70 mN was delivered and a thruster efficiency of 20% for an rf power of 6 kW was achieved.[9]

In the magnetic nozzle, a Lorentz force, attributed to the azimuthal current and the magnetic field, is exerted on the plasma, which produces thrust in the axial direction. Previous studies implied that the azimuthal current in electrodeless plasma thrusters was mainly the result of a diamagnetic effect, both theoretically and experimentally, when varying the magnetic field strength.[10, 11, 12] In addition, the azimuthal current in the electrodeless plasma thruster was measured experimentally, with the results indicating that the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ drift current decreased with increasing magnetic field strength and the diamagnetic effect became dominant.[13] The magnetic field gradient has also been demonstrated to significantly affect the electron transport.[14, 15] However, the experiment did not accurately evaluate other drifts such as grad-$B$ and curvature drifts.

The diamagnetic effect in the magnetic nozzle was analyzed using a two-fluid model.[16] Moreover, particle-in-cell simulations were also conducted to investigate potential structures in the magnetic nozzle.[17, 18] These numerical simulations focused on the downstream region of the thruster and did not simulate the plasma source; however, the connection between the two different regions of the plasma source and the magnetic nozzle is the key element that would need to be analyzed to understand the plasma dynamics in an electrodeless plasma thruster that uses a magnetic nozzle. In our previous studies, we conducted particle-in-cell simulations of an electrodeless plasma thruster with Monte Carlo collisions (PIC-MCCs) to investigate the axial momentum loss to the lateral dielectric wall and the effects of neutral distributions. These simulations focused on the source region by employing a closed boundary for the source tube exit.[19, 20] However, PIC-MCC simulations that take into account both the plasma source and the magnetic nozzle acceleration have not been performed yet. These simulations would obviate the need for assumptions about the density and velocity distributions for the plasma injection at the magnetic nozzle.

In this study, we conducted PIC-MCC simulations that include both the source and magnetic nozzle regions of the electrodeless plasma thruster and investigated the internal plasma current. Here, PIC-MCC simulations can accurately and self-consistently derive the azimuthal current induced by charged particles in the magnetic nozzle. The $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ and diamagnetic drift currents can also be computed from the potential and pressure profiles of the electrons and are compared with the net current profile; the results demonstrate that the internal azimuthal current is diamagnetic in the case of a strong magnetic field. The numerical simulation additionally demonstrated that the Lorentz force induced by the net current and the magnetic nozzle can be enhanced by increasing the magnetic field strength, which is consistent with the previous observation.[12, 13]

In Sec. II, the numerical model employed in this study is briefly described. Sec. III presents the calculation of the net electron current density using the electron motions for three solenoid currents. These results are compared with the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ drift current resulting from the electrostatic field and the diamagnetic drift current due to the electron pressure gradient. As the solenoid current increases, the net electron current density is increasingly dominated by the diamagnetic effect. In addition, the Lorentz force induced by the net electron current density is discussed.

We decided to employ a two-dimensional and symmetric calculation model to avoid singularities on the central axis and to reduce the calculation cost compared with the axisymmetric model. Figure 1 shows a schematic of the area of the electrodeless plasma thruster that was the subject of this study. The calculation model consisted of an rf antenna, a dielectric plate, and a solenoid. The rf antenna supplies power to the plasma, and the dielectric plate confines the plasma to the $y$-direction. The solenoid produces a magnetostatic field that forms a magnetic nozzle that accelerates the plasma. The size of the area included in the calculation is 2.5 cm $\times$ 0.56 cm; this area contains a dielectric with a length of 1.5 cm and a thickness of 0.05 cm. Only the solenoid magnetic field is solved using a calculation area that is ten times larger. It should be noted that the height of the discharge chamber employed in our simulation is 1 cm, which is approximately one-sixth of the diameter of the thruster (6.4 cm in this experiment).[13] The calculation model is symmetric about the $x$- and $y$-axes such that the plasma is extracted bidirectionally in this simulation. Although the bidirectional thruster differs from the normal unidirectional thruster, the plasma physics in the magnetic nozzle within $x$ = 1.5–2.5 cm would remain essentially unchanged. In addition, bidirectional thrusters have also been proposed and have been demonstrated for space debris removal.[21] The calculation area is divided into 50 $\mu$m $\times$ 50 $\mu$m cells. In Sec. III, we focus on the downstream region of the thruster to investigate the internal plasma currents in the magnetic nozzle and we present only the results that pertain to the area demarcated by the dashed red rectangle in Fig. 1. The simulation results presented in Sec. III are averaged over 30 $\mu$s after reaching the steady state.

Figure 2 shows the calculation flow for the PIC-MCC simulation. Our PIC-MCC simulation mainly considers the kinetics of charged particles, electron-neutral collisions, weighting for the charge, and an electrostatic field. Details of the PIC-MCC simulation appear in our previous papers[19, 20]; hence, we only include a brief description of our PIC-MCC simulation in this paper.

The equations of motion of the charged particles are written as

The potential generated by true charges (charged particles and surface charge) $\phi_{t}$ satisfies Poisson’s equation as follows:

The complex electric field induced by the rf antenna $\tilde{E}_{em}$ is written as

The solenoid magnetic field $B$ is written as

We treat elastic, excitation, and ionization collisions as electron-neutral collisions in the simulation using the null-collision method.[26] Here, ion-neutral collisions are neglected for simplicity because the mean free path for ion-neutral collisions is estimated to be 5 cm using a cross-section of 1 $\times$ 10^-18 m^-3 and a neutral density of 2 $\times$ 10¹⁹ m^-3. Many ions do not collide with neutral particles in the calculation area of 2.5 cm $\times$ 0.56 cm. The density of these neutral particles was set to $2\times 10^{19}$ m^-3, and the temperature of these particles was set to 300 K.

The PIC simulation can directly provide the internal plasma current density by considering the velocity components of the charged particles. The plasma current determined by the fluid model was discussed in a previous report.[13] To compare the simulation results with the results of the experiment, the electron pressure $p_{e}$ is calculated from

The net electron current density in the $z$-direction $j_{e,z}$ is calculated from the electron motions in the PIC-MCC simulation as

We obtain the Lorentz force density in the $x$-direction $f_{x}$, which is exerted on the electrons and calculated as

Figures 3 and 4 show the $x$-$y$ profiles of the electron number density $n_{e}$ and the electron temperature $T_{e}$, respectively, for the three solenoid currents of 0.1, 0.4, and 2.0 kA. Then, the electron pressure $p_{e}$ is calculated from the electron number density $n_{e}$ in Fig. 3 and the electron temperature $T_{e}$ in Fig. 4. Figure 5 shows the $x$-$y$ profiles of the electron pressure $p_{e}$ for the three solenoid currents of 0.1, 0.4, and 2.0 kA. The $y$ profiles of the electron pressure $p_{e}$ for the solenoid currents of 0.1 and 0.4 kA indicate center-peaked plasmas as shown in Fig. 5(a) and 5(b), but that for the solenoid current of 2.0 kA indicates a bimodal shape as shown in Fig. 5(c). The reason for the bimodal shape is that energetic electrons are generated in a peripheral region near the rf antenna and the strong magnetic field prevents the plasma from being transported on the central $x$-axis, well reproducing the previous experimental results.[27, 28, 29, 30, 31] The maximum electron pressure increases as the solenoid current $I_{B}$ increases because the strong magnetic field confines the plasma more efficiently. A marked diamagnetic effect would be observed in regions where the electron pressure gradient $\nabla p_{e}$ is large. Here, because both the cylindrical thruster in experiments and the Cartesian one in simulations give rise to plasma with a bimodal shape, the cylindrical and Cartesian coordinates are considered not to have any essential differences.

We compare the electron pressure $p_{e}$ computed with our PIC-MCC simulation with the experimental results from a previous study, which were measured using a Langmuir probe and is shown in Fig. 6(d) of Ref. [13]. Figure 6 shows the $\hat{y}$ profiles of the normalized electron pressure $\hat{p}_{e}$ for the two solenoid currents of 0.4 and 2.0 kA, where the coordinate $y$ is normalized by the maximum position $y_{max}$ as $\hat{y}=y/y_{max}$ and the electron pressure $p_{e}$ is normalized by the maximum pressure $p_{e,max}$ as $\hat{p}_{e}=p_{e}/p_{e,max}$. The electron pressure, $p_{e}$, in Fig. 6 was calculated using Eq. (14), where $n_{e}$ and $T_{e}$ in the experiments are obtained from the current-voltage characteristics assuming a Maxwellian electron energy distribution function. Here, the simulation result is the cross-section of $x$ = 2.0 cm at $I_{B}$ = 0.4 and 2.0 kA, and the experimental result is that of $z$ = 5 cm at $I_{B}$ = 2 and 12 A, respectively. The position of $x$ = 2.0 cm in the simulation roughly corresponds to that of $z$ = 5 cm in the experiment because the size of the calculation model is approximately one-sixth as described in Sec. II. The normalized electron pressure $\hat{p}_{e}$ obtained from our simulation is in good agreement with the experimental results. The simulation results shown in Fig. 6(b) indicates the low pressure at the center of $\hat{y}$ = 0 compared with the experimental result because the calculation model is two-dimensional ($x$-$y$), whereas actual thrusters have a cylindrical shape. The cylindrical discharge chamber confines the plasma to the radial direction. However, the two-dimensional model with its rectangular coordinates is less capable of structurally confining the plasma to the center because free vertical motions exist in the $z$-direction. The electron pressure gradient $\nabla p_{e}$ becomes large compared with the experimental result; thus, our two-dimensional model slightly overestimates the diamagnetic effects in the central region. However, the normalized electron pressure $\hat{p}_{e}$ obtained from our simulation is qualitatively consistent with the experimental result, indicating a confirmation of the validity of the simulation results, although it may not quantitatively capture the physics of cylindrical thrusters.

In these simulations, the ratio $r_{L}/L$ is roughly adjusted between the simulations and the experiments, where $r_{L}$ is the Larmor radius of ions, $L$ is the characteristic length, and $L$ is set to the thruster height (diameter in cylindrical coordinates). However, the parameters including the ratio $r_{L}/L$ are not scaled up or down accurately in the simulations. Instead of scaling, we consider the distribution of the electron pressure to be a key parameter for validating the experimental results. If the plasma pressure distributions are the same, the simulations are expected to reproduce the internal plasma current and force quantitatively. In addition to the plasma profiles, the Larmor radii were roughly the same between the simulations and the experiments, and the Hall parameters were in excess of a few hundred in both cases, indicating that there were no essential differences between the simulations and the experiments.

Figure 7 shows the $x$-$y$ profiles of the net electron current density $j_{e,z}$ for the three solenoid currents of 0.1, 0.4, and 2.0 kA, where the net electron current density $j_{e,z}$ is calculated using Eq. (16). The net electron current densities in the upper and lower half regions are opposite to each other in Figs. 7(a) and 7(b). However, in Fig. 7(c), the opposite net electron current density $j_{e,z}$ only exists in the upper half region. The reason for this unusual phenomenon is mentioned later in this paper. Although the net electron current density does not follow the magnetic field lines in the case of the weak magnetic field strength as shown in Fig. 7(a), they are fairly aligned with the magnetic field lines for $I_{B}$ = 0.4 kA as shown in Fig. 7(b). Figure 7(c) shows that the net electron current density is clearly distributed along the magnetic field lines because of the strong magnetic field.

Here, we define the characteristic magnetic field strength as the magnetic field strength under the solenoid at $x$ = 1.1 cm and on the $x$-axis. For solenoid currents of 0.1, 0.4, and 2.0 kA, the characteristic magnetic field strengths are approximately 5, 20, and 100 mT, respectively. Then, the respective Larmor radii become 1.70, 0.43, and 0.08 mm at the electron temperature of 5 eV. In addition, the Hall parameters were estimated to be more than one hundred for any solenoid current. Therefore, electrons in the magnetic nozzle are not likely to be transported across the magnetic field lines by elastic collisions. Moreover, because the Larmor radius for the solenoid current of 2.0 kA is sufficiently small compared with the thruster size, the plasma profile is split into the two regions.

We derived the potential distributions to investigate the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ effect. Figure 8 shows the $x$-$y$ profiles of the potential $\phi$ for the three solenoid currents of 0.1, 0.4, and 2.0 kA. The potential decreases from the plasma source on the left side to the downstream on the right side at all solenoid currents. The potential gradient $\nabla\phi$ and the solenoid magnetic field $B$ would be expected to induce an $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ drift current. It should be noted that there are sheaths near the right boundary at $x$ = 2.3–2.5 cm and the top and bottom boundaries at $y$ = $\pm$(0.3–0.56) cm because the potential $\phi$ is solved with the Dirichlet boundary condition as $\phi$ = 0.

We calculate the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ drift current density $j_{E\times Be}$ and the diamagnetic drift current density $j_{De}$ for comparison with the net electron current density $j_{e,z}$, which are expressed in terms of the two-dimensional fluid model as

Figure 11 shows the $y$ profiles of the net electron current density $j_{e,z}$, the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ drift current density $j_{E\times Be}$, and the diamagnetic drift current density $j_{De}$ at $x$ = 1.8 cm for the three solenoid currents of 0.1, 0.4, and 2.0 kA. Although both the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ and the diamagnetic effects exist in all the cases, they are canceled out because they exist in opposite directions to each other, as is obvious from the net electron current density $j_{e,z}$. In the case of the solenoid current of 0.1 kA, the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ and diamagnetic current densities are between $\pm$0.25 kA/m², which is small compared with those of 0.4 and 2.0 kA. Although the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ drift current density $j_{E\times Be}$ is larger than the diamagnetic one $j_{De}$ at $I_{B}=0.1$ kA, the diamagnetic drift current density $j_{De}$ is larger than the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ one $j_{E\times Be}$ for both the solenoid currents of 0.4 and 2.0 kA. The effect of the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ drift current density $j_{E\times Be}$ decreases as the solenoid current increases, and the diamagnetic drift current density $j_{De}$ becomes dominant. In addition, other drifts such as grad-$B$ and curvature drifts are negligibly small in the presence of the strong magnetic field because the net electron current density $j_{e,z}$ is almost equal to the diamagnetic current density $j_{De}$ as shown in Fig. 11(c).

The Dirichlet boundary condition at $x$ = 2.5 cm and $y$ = $\pm$0.56 cm, which imposes a zero potential, affects the simulation results. Here, a marked sheath effect is expected to be seen in the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ current in Fig. 9, because it is directly determined by the electrostatic field $\mbox{\boldmath$E$}_{es}$. In Fig. 9, anomalous current densities exist in the regions $x$ = 2.3–2.5 cm and $|y|$ ¿ 0.3 cm, and these are not suitable for validating the experimental results. However, anomalous current densities do not exist in $x$ = 1.5–2.3 cm and $|y|$ ¡ 0.3 cm, indicating that the sheath effect is sufficiently small in these regions. Because the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ current density is larger when the solenoid current is 0.1 kA as shown in Fig. 11(a), the sheath effect may be relatively large. However, when the solenoid current is 2.0 kA, the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ current density is negligibly small, as shown in Fig. 11(c), and the sheath does not affect the result, i.e., the diamagnetic effect dominates the net current density.

The $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ drift current density $j_{E\times Be}$ is proportional to the electrostatic field $E_{es}$ and inversely proportional to the magnetic field strength $B$ as indicated in Eq. (19): The electrostatic field $\mbox{\boldmath$E$}_{es}$ tends to remain unchanged with increasing solenoid current $I_{B}$, as shown in Fig. 8, and thus the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ drift current density $j_{E\times Be}$ decreases as the solenoid current $I_{B}$ increases.

The diamagnetic drift current density $j_{De}$ is proportional to the electron pressure gradient $\nabla p_{e}$ and inversely proportional to the solenoid magnetic field $B$, as indicated in Eq. (20). The electron pressure gradient $\nabla p_{e}$ increases significantly as the solenoid current increases from 0.1 to 0.4 kA because the maximum electron pressure increases tenfold as shown in Figs. 5(a) and 5(b). This higher rate of increase in $\nabla p_{e}$ than $B$ results in a more than two-fold increase in $j_{De}$. However, the increase in the maximum electron pressure from 0.4 to 2.0 kA is only twofold, and the plasma width becomes roughly half in the $y$-direction because of the bimodal shape, approximately quadrupling the electron pressure gradient. Thus, the increasing rate of $\nabla p_{e}$ is approximately equal to that of $B$, leading to a slight increase in $j_{De}$.

The net electron current density $j_{e,z}$ is dominated by the diamagnetic drift current density $j_{De}$ at $I_{B}$ = 2.0 kA as shown in Fig. 11(c). A similar trend was observed in a previous experiment (see Figs. 6(e) and 6(f) in Ref. [13]). Our simulation results were in good agreement with the experimental results.

Figures 12–14 show the $y$-profiles of the sum of the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ and diamagnetic drift current density $j_{E\times Be}+j_{De}$ and the net electron current density $j_{e,z}$ at $x$ = 1.8, 2.0, and 2.2 cm for the three solenoid currents of 0.1, 0.4, and 2.0 kA. Note that $j_{E\times Be}$ + $j_{De}$ for the solenoid current of 0.1 kA at $x$ = 2.0 and 2.2 cm include anomalous results near $y$ = 0 cm. For the solenoid current of 0.1 kA, $j_{E\times Be}+j_{De}$ is equal to $j_{e,z}$ at $x$ = 1.8 cm for $|y|$ ¡ 0.3 cm, but they do not correspond at $x$ = 2.0 and 2.2 cm. Therefore, $j_{e,z}$ is composed of not only $j_{E\times Be}$ and $j_{De}$ but also other drift currents, e.g., grad-$B$ and curvature drifts. For the solenoid current of 0.4 kA, $j_{E\times Be}+j_{De}$ is roughly equal to $j_{e,z}$ at all locations. In addition, $j_{E\times Be}+j_{De}$ is completely equal to $j_{e,z}$ for the solenoid current of 2.0 kA except in the peripheral region within —$y$— ¿ 0.3 cm, where the sheath is expected to affect the plasma. Because the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ effect decreases as the solenoid current increases as shown in Fig. 11, increasing the solenoid current suppresses the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ effect and other effects such as grad-$B$ and curvature drifts, and the diamagnetic effect becomes dominant instead.

We obtain the Lorentz force density exerted on the electrons in the $x$-direction ($f_{x}$), which is calculated from the net electron current density $j_{e,z}$ and the solenoid magnetic field in the $y$-direction $B_{y}$, as shown in Eq. (17). Figure 15 shows the $x$-$y$ profiles of the Lorentz force density in the $x$-direction ($f_{x}$) exerted on the electrons for the three solenoid currents of 0.1, 0.4, and 2.0 kA. The Lorentz force density in the $x$-direction ($f_{x}$) increases dramatically with an increase in the solenoid current $I_{B}$. However, the magnitude of the net electron current density $j_{e,z}$ does not change significantly when the solenoid current $I_{B}$ increases, especially in the case of the solenoid currents of 0.4 and 2.0 kA as shown in Figs. 11(b) and 11(c). Here, the Lorentz force density in the $x$-direction ($f_{x}$) is the product of the net electron current density $j_{e,z}$ and the solenoid magnetic field in the $y$-direction $B_{y}$, as shown in Eq. (17); additionally, the magnetic field strength $B$ of the solenoid increases with an increase in the solenoid current $I_{B}$. Therefore, the strong magnetic field pushes electrons downstream electromagnetically, although the electron current density does not change significantly with the solenoid current $I_{B}$.

Figure 16 shows the $y$ profiles of the Lorentz force density in the $x$-direction ($f_{x}$) exerted on electrons at $x$ = 1.8 cm, 2.0 cm, and 2.2 cm for the three solenoid currents of 0.1, 0.4, and 2.0 kA. The Lorentz force density in the $x$-direction ($f_{x}$) is almost zero for the solenoid current $I_{B}$ = 0.1 kA at all locations, whereas a positive Lorentz force density is obtained for the solenoid current of $I_{B}$ = 0.4 kA. In the case of $x$ = 1.8 cm and the solenoid current of 2.0 kA, the positive Lorentz force density in the $x$-direction is 6.43 N/m³ at maximum, whereas the negative Lorentz force density in the $x$-direction is exerted on the electrons between $y$ = $\pm$ 0.1 cm. Therefore, the positive Lorentz force density in the $x$-direction increases as the solenoid current $I_{B}$ increases, whereas the negative density emerges in the central region.

Table 1 shows the magnetic nozzle thrust $F_{x}$ calculated by Eq. (18), where the Lorentz force density in the $x$-direction $f_{x}$ is integrated in 1.5 cm ¡ $x$ ¡ 2.3 cm and $|y|$ ¡ 0.3 cm to reduce the sheath effects. The magnetic nozzle thrust $F_{x}$ increases as the solenoid current $I_{B}$ increases. At the solenoid current $I_{B}$ = 0.1 kA, the magnetic nozzle thrust $F_{x}$ is $-$1.22 $\mu$N/m, indicating that the magnetic nozzle does not accelerate the plasma in the downstream direction. The negative magnetic nozzle thrust is due to the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ drift current because it occurs in the direction opposite to the diamagnetic current and $-j_{E\times Be}B_{y}$ is in the negative $x$-direction. Note that the total thrust is not negative because the ions accelerated by the electrostatic force impart a positive thrust and the magnetic nozzle thrust is composed of only electron momentum. At the solenoid currents of 0.4 and 2.0 kA, positive magnetic nozzle thrusts are obtained. In addition, the magnetic nozzle thrust at $I_{B}$ = 2.0 kA is about 1.41 times larger than that at $I_{B}$ = 0.4 kA. Therefore, the strong magnetic field generates larger thrust, even though the thrust is also produced in the negative thrust region. The net electron current density $j_{e,z}$ remains almost unchanged as the solenoid current increases, suggesting that the strong magnetic field would increase the magnetic nozzle thrust $F_{x}$.

The magnetic nozzle thrust $F_{x}$ does not depend directly on the mass ejection, as indicated in the Tsiolkovsky equation. In the magnetic nozzle, the electron pressure is converted to ion momentum by electrostatic fields such as a current-free double layer or an ambipolar electric field, and the thrust is finally obtained by the ion ejection [32]. However, the conversion of momentum between ions and electrons has not been clarified in detail and would have to be investigated in the future.

We conducted two-dimensional PIC-MCC simulations of an electrodeless plasma thruster using a magnetic nozzle for three solenoid currents and investigated the dependence of the internal plasma currents on the magnetic field strength. The diamagnetic drift current density becomes dominant with increasing solenoid current (i.e., magnetic field strength), whereas the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ drift current density decreases, indicating that the net electron current is caused due to the diamagnetic effect. Other drifts such as grad-$B$ and curvature drifts are observed for a solenoid current of 0.1 kA in the simulation, but their effects diminish as the solenoid current increases. These results are consistent with the previous experimental results, which implied that the dominant drift changed from the $\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}$ effect to the diamagnetic one as the magnetic field became stronger. The Lorentz force density has positive and negative regions at the highest solenoid current that was examined ($I_{B}$ = 2.0 kA), but the net magnetic nozzle thrust is positive and increases with increasing solenoid current. Therefore, our simulations demonstrated that the diamagnetic effect significantly contributes to thrust generation in the magnetic nozzle. Because the simulations with the sheath under the Dirichlet boundary condition reproduced the internal plasma currents measured in previous experiments, the sheath effect is not expected to be essential for the diamagnetic effect and thrust generation in the magnetic nozzle. However, the Dirichlet boundary condition is not suitable for validating the plasma in space operations. In the future, it would also be necessary to validate the internal plasma currents for operations in space with open boundary conditions, for example, the Neumann boundary condition.