Contents
CHAPTER 1: INTRODUCTION ························ 1
1-1 History of space explorations and future target ································································· 1
1-2 Application of magnetic nozzle for space propulsion system ··················································· 2
1-3 Principle of magnetic nozzle and plasma detachment··························································· 2
1-4 Plasma detachment scenario ························································································ 4
1-5 Problem in previous research on plasma detachment ··························································· 5
1-6 Research object ········································································································ 6

CHAPTER 2: NUMERICAL CALCULATION METHOD ············· 7
2-1 The basic formulas of numerical analysis ········································································· 7
2-2 Calculation flow ······································································································· 7
2-3 Calculation of the fields ······························································································ 8
2-3-1 Maxwell’s equations ······························ 8
2-3-2 Finite-Difference-Time-Domain (FDTD) scheme ··················· 9
2-3-3 Boundary condition ······························· 11
2-3-4 Initial calculation for preparing a background magnetic field ············· 11
2-4 Calculation of particles····························································································· 12
2-4-1 Particle shape function ····························· 12
2-4-2 Equation of motion······························ 13

2-4-3 Boundary condition ······························· 15
2-5 Calculation of the current ························································································· 15
2-6 Super particle ········································································································ 15
2-7 Time step and grid length ·························································································· 15
2-8 Parallel calculation method ······················································································· 16
2-9 Parallelization in EPOCH code ··················································································· 16

CHAPTER 3: TEST CALCULATION RESULTS ················ 17
3-1 Calculation model ··································································································· 17
3-2 Verification of EPOCH code ······················································································ 18
3-3 Validation of a magnetic nozzle working on EPOCH code ··················································· 21
3-4 Investigation of calculation time and use of parallelization ·················································· 23

CHAPTER 4: SIMULATION RESULTS AND DISCUSSION ··········· 27
4-1 Calculation condition ······························································································· 27
4-2 Time evolution of plasma number density ······································································ 28
4-4 Analysis about the cyclotron radius ·············································································· 35
4-4-1 Estimation of thermal velocity ·························· 35
4-4-2 Analysis results ································ 36
4-5 Discussion ············································································································ 43

CHAPTER 5: CONCLUSION ························ 47
REFERENCES ······························ 48
ACKNOWLEDGEMENTS ························· 50

Chapter 1

Chapter 1: Introduction
1-1 History of space explorations and future target
The space exploration has brought about great clues to understand the origins of life and thecreation of the solar system. The first lunar landing was achieved in the Apollo program of the
National Aeronautics and Space Administration (NASA) in 1969. The research group found that
there is lunar regolith on the surface. Lunar regolith consists of oxygen, silicon, iron, calcium,aluminum, magnesium, and others [1]. The existence of these elements explains that oxygen andhydrogen oxide may be produced from cosmic materials. Hayabusa [2], a spacecraft of the Japan
Aerospace Exploration Agency (JAXA), successfully brought back a sample of materials from small
near-earth asteroid named Itokawa in 2010. The Hayabusa used ion engines to gain a propulsion
force (< 10 mN), and the sample return mission took 7 years. The research group defined that
Itokawa consists from some components same as meteorites [3]. NASA has planned to explore Mars
as a manned space mission [4]. Mars has carbon dioxide and iced water, meaning that human could
live on it in future. Therefore, Mars manned mission is worth exploring to expand human’s habitats.
In addition, European Space Agency (ESA) have planed the Jupiter Icy Moons Explorer (JUICE)[5]. The satellites of the survey object have iced water, so that an organism may exist.
The thrust and fuel consumption rate is a trade-off relation in the existing thrusters as shownin Fig. 1.1. A chemical propulsion has the highest thrust in existing thrusters, so it has been used tolaunch from a surface of the earth to space. However, the chemical propulsion shows low specificimpulse that is defined as ejection velocity over gravitational acceleration (Isp = Vejection/g sec), so itneeds a lot of fuels to travel a long roundtrip such as Mars and Jupiter missions. In contrast, anelectric population such as ion thruster has high specific impulse, so it has been used for unmanned
space exploration mission. However, the mission span takes many years, astronauts will takesomatic damages from cosmic rays and mental damages from being restrained in a small living
space of spaceship when the existing electric propulsion systems are used for manned missions.
Hence, an advanced propulsion system that has large thrust flux density (> 10 5 Nm-3, it iscomparable to the chemical propulsion) and low fuel consumption (Isp > 104 sec, it is comparable tothe electric propulsion) is indispensable to achieve manned deep space missions such as Mars andJupiter.

Figure 1.1: The relationship between thrust flux density and specific impulse [6]. The trust flux density is defined as
thrust per unit a volume of injection outlet, Nm-3. A chemical propulsion has the highest range on this value, so it
has been used to launch from a surface of the earth to space. The specific impulse is defined as ejection velocity
over gravitational acceleration (Isp = Vejection/ g sec). A high value of this means that a fuel consumption efficiencyof the thruster is low, it can travel a long distance with a few propellants such as an ion thruster. The thrust andspecific impulse is trade-off relationship in the current propulsion systems.

1Chapter 1

1-2 Application of magnetic nozzle for space propulsion system
Chemical rocket engines obtain large thrust from a chemical reaction simply. The thermalenergy is converted to a kinetic energy on a rear part of the engines by using a Laval nozzle as shownin Fig. 1.2. The nozzle converts a subsonic flow to a supersonic flow to extract fuel from the nozzle.
A concept of magnetic nozzle is similar to the Laval nozzle as a magnetic field forms the shape ofnozzle instead of a solid wall. This magnetic wall captures charged particles and pushes them back.
This basic effect is explained by a Lorentz force induced from the magnetic field. For example, amagnetic mirror that is one of magnetic confinement methods to generate a fusion power in fusion
research to trap a high temperature plasma using magnetic fields. The principle of magnetic nozzle
is explained in the next section. The magnetic nozzle has been planned to apply to some plasma
rockets: Magneto-Plasma-Dynamic (MPD) thruster [7-8], Helicon plasma thruster [9], and Laser
Fusion Rocket (LFR) [10-14]. A LFR has high ejection velocity and high ejection mass rate because it
uses laser-fusion that generate a large energy and high density plasma. Therefore, the LFR has apotential to realize manned explanatory explorations such as Mars mission.

Figure 1.2: The Laval nozzle and magnetic nozzle [15]. The Laval nozzle transforms a subsonic ﬂow into asupersonic ﬂow using a throat in an interior ﬂow. The magnetic nozzle is applied to a plasma thruster. The plasma
flow direction is guided to the magnetic field, and a thermal energy is converted to a kinetic energy.

1-3 Principle of magnetic nozzle and plasma detachment
The magnetic nozzle controls a plasma velocity using a dipole magnetic field. A charged
particle receives the Lorentz force: qVB where an external magnetic field is applied. The force
causes a particle rotation, makes a swirl in a plasma. This behavior is called as the cyclotron motion,
following Eq. (1-1) shows the cyclotron radius of charged particles as follows:𝑚𝑉

𝑟c = |𝑞|𝐵⊥ ,(1-1)

where m is the particle mass, V⊥ is the perpendicular particle velocity to the magnetic field line, q
is the particle charge, and B is the magnetic flux density. The cyclotron radius is inversely
proportional to B, therefore, rc becomes smaller where a magnetic field becomes stronger since V⊥
increases to a rotational angle by the law of conservation of angular momentum. These radius andvelocity change on the upstream of magnetic nozzle, especially, near the throat. In addition, a plasma
receives a JB force induced from a reaction of the swirl current and a magnetic field, and it is2

Chapter 1
accelerated to axial direction: ejection direction. The rotational and axial kinetic energy of a plasma
are conserved in the magnetic nozzle, so V⊥ is converted to V‖ in the downstream. In other word,
the thermal energy of plasma is converted to the kinetic energy through the magnetic nozzle.
The LFR uses this mechanism that converts the plasma velocity to gain thrust. As shown in
Fig.1.3, this system uses inertial confinement fusion (ICF) plasma generated by a laser-fusion andan external dipole magnetic field generated by an electromagnetic coil. First, the plasma is generated
in the magnetic field and expands isotopically. Second, the plasma forms a magnetic cavity insideand compresses the magnetic field outside. After a balance of energy between the plasma and themagnetic field, the compressed magnetic field pushes back the plasma in the right direction in the
figure. As a result, the coil is accelerated to the left direction in the figure by a reaction force.
However, there is an unclear point in the magnetic nozzle when it is applied to a space propulsion.
Each particle such as electron and ion moves along the magnetic field keeping the cyclotron motion.
This means that the ejected plasma may come back to the spaceship along with a looped magnetic
field line. As shown in Eq. (1-1), rc is directly proportional to the particle mass. It means that theion cyclotron radius rci is 1800 times larger than the electron cyclotron radius rce. The difference ofrc implies that a magnetic field is more effective to electrons than ions, and it may cause a different
charge distribution in space. If only ions detach from a magnetic field and electrons move along amagnetic field line as shown in Fig. 1.4(a), an electric field is induced between the particles asshown in Fig. 1.4(b). Furthermore, if the ion is pulled back by the Coulomb force as shown in Fig.
1.4(c), the ejected plasma returns to the magnetic nozzle. It means the magnetic nozzle isinappropriate for a space propulsion system. Although we need to find a condition that both ion andelectron detach from a magnetic field: Plasma Detachment condition, as shown in Fig.1.4(d), it has
not been elucidated.

Figure 1.3: Plasma ejection mechanism used magnetic nozzle in LFR. (a) Plasma is generated in the magnetic
field and expands isotopically. (b) Plasma forms a magnetic cavity inside and compresses the magnetic fieldoutside. (c) Induced current in plasma and magnetic field generate the Lorentz force JB in the right direction andplasma is ejected by the force. As a result, the coil is accelerated in the left direction by a reaction force.

3

Chapter 1

Adapted Kawabuchi et al., 2008[16]Figure 1.4: A problem caused by the deference of mass between ion and electron. (a) Only ion moves separatingfrom a magnetic flux line and an electron moves along a magnetic flux line. (b) An electric field is induced between
the particles. (c) The ion is pulled back by the Coulomb force induced the electric field, the ejected plasma returns
to the magnetic nozzle along magnetic field line. (d)Both ion and electron should detach from a magnetic field line;
Plasma Detachment.

1-4 Plasma detachment scenario
One of the candidate that causes a plasma detachment is demagnetization or loss of adiabatic
invariant. As shown in Fig. 1.5, charged particles show the cyclotron motion in a magnetizes plasma
as explained in Eq. (1-1). If the cyclotron orbit changes drastically over a single cyclotron period,
the cyclotron orbit becomes demagnetized when ∆rc  rc. Also, the magnetic moment is lost when
the local gradient of B becomes large on the length that a plasma walks over a single cyclotron
period. This can be expressed as follows:
,

(1-2)

where ωc is the cyclotron frequency and Δ refers to a time change in a single cyclotron period.
Since the particles spiral along the magnetic flux line until  is no longer conserved, Δωc/ωc canbe estimated using the path along the magnetic field line[17] as follows:
,

(1-3)

where Δs is the distance that a particle travels along the field line during the single cyclotron period,and b is the unit vector along the magnetic flux line. If most of the energy is converted fromperpendicular to parallel while the magnetic moment is still conserved, the distance can beapproximated as follows:.
(1-4)
Using this result together with Eq. (1-2), the approximate range where adiabatic invariance is valid
is estimated as follows [18]:4Chapter 1
.

(1-5)

Since the magnetic moment is conserved when the second term in Eq. (1-5) is less than unity, theadiabatic invariant decreases when the product of the cyclotron radius and the magnetic scale length:
exceeds unity.

Figure 1.5: Detachment by breakdown of the adiabatic invariant [15]. The particle cyclotron radius increases in an
expanding dipole magnetic field and the perpendicular velocity to the magnetic field is converted into parallel
velocity. When the cyclotron radius changes from an orbit to the next in a single cyclotron period, the particles
become demagnetized and continue on the same trajectory at the breakdown of moment.

1-5 Problem in previous research on plasma detachment
Terasaka et al. have measured the flow velocity of ions in an electron cyclotron resonance(ECR) plasma in a diverging magnetic field, where the electron is magnetized in the whole plasma
region while the ions are slightly magnetized in the downstream region [19, 20]. They obtained theflow velocity and also measured the detachment of ion from the magnetic field. They concluded thatthe ion detachment takes place when the parameter |fciLB/Vi| which equals to (rci|B|/B)-1 in Eq. (15), becomes order of unity. Olsen et al. have applied an experiment using the VASIMR VX-200
device operating at a power level of 100 kW to map the behavior of plasma ﬂowing through amagnetic nozzle to investigate the detachment process. [18] They concluded that plasma detachment
is well described by two processes. The first part involves the separation of the ion from the magnetic
ﬁeld through the breakdown of the ion magnetic moment. The second part involves turbulence,created by electrostatic two-stream instabilities, where a ﬂuctuating electric ﬁeld facilitatescompeting interactions between detached ion and magnetized electrons. As the magnetic ﬁeld
weakens, the detached ions begin to dominate and the electrons respond to the turbulent electricﬁeld. Therefore, the electrons cross the magnetic ﬁeld lines pursuing the ion until the turbulent
electric ﬁeld dissipates, this transport may cause an electron detachment.
These two research groups have explained the ion detachment as the breakdown of the ion
magnetic moment in the experiments. In contrast, the mechanism of electron detachment is notclarified. A numerical simulation as well as experiment has conducted to investigate a condition of
the electron detachment. Kawabuchi et al. have conducted a numerical simulation using a twodimensional electromagnetic fully Particle-In-Cell (PIC) code TRISTAN in order to clarify theplasma behavior and the possibility of the electron detachment from a magnetic nozzle. [14] They
have shown the particle trajectories that ions move with the breakdown of the magnetic moment,
electrons also move to follow the ions. However, they assumed that the mass ratio of ion and electron
mi/me was 10. The total calculation region size was too small since the calculation costs of PIC
simulation is high and they could not handle the high cost calculation. Therefore, they could not find5

Chapter 1
the difference of ions and electrons in the calculation and have not clarified the plasma detachment
physics. Hence, a numerical simulation model is needed that solve with the both of a high calculation
cost and an accuracy of physics to understand the plasma detachment mechanism, especially, about
electron detachment.

1-6 Research object
This research object is to develop a calculation model in order to understand the physics ofplasma detachment, especially, the loss of adiabatic invariant of electron. Chapter 2 explains about
the numerical calculation method, problems in a current calculation, and solutions of these problems.
Chapter 3 shows the calculation results for a checking calculation code and an evaluation of parallel
calculation method. Chapter 4 shows simulation results, the analysis, and the discussion about the
behavior of the plasma in a magnetic nozzle. Finally, summary and conclusion are shown in chapter
5.

6Chapter 2

Chapter 2: Numerical Calculation Method
Here, an investigation of individual particle motions is needed to analysis the behavior of ions and electrons
in a magnetic field in this research. Therefore, a PIC simulation is used to calculate on first principles due to hightemperature and low density plasma. The mean collision period is long compared to the calculation period, so we
ignored the effect of collision. PIC simulation is a self-consistent calculation method of particle movements in thecalculation space, and resultant current, as well as magnetic field and the electric field. In this research, EPOCH
code [24] was used for PIC simulation code. However, a calculation cost of PIC simulation increases by the total
number of particles in the calculation, it is difference to take long calculation time. Therefore, a parallel calculation
using many CPU core was needed to solve this problem. In this chapter, the fundamental equations for plasma are
explained. In addition, the equations used in the PIC code and the finite-difference methods are explained. Finally,
the method of calculation cost reduction is explained in the following sections.

2-1 The basic formulas of numerical analysis
The equation of motion is used as a basic formula for analyzing the behavior of plasma. The
formula is shown as follows:
,
𝑑𝒙
𝑑𝑡

=𝑽 ,

(2-1)
(2-2)

where m is the particle mass, V is the velocity of the particle, x is the position of the particle, t is
the time, q is the charge of the particle, E is the electric field, and B is the magnetic flux density.

2-2 Calculation flow
EPOCH code mainly solves five variables; E, B, x, particle momenta P, and current J. Figure
2.1 shows the calculation ﬂow diagram on a single time step, n to n+1. First, the electric ﬁeld E n
and the magnetic field B n are advanced with the help of the current J n by half a time step to E
n+1/2
and B n+1/2. Second, the particle position x n is advanced with the help of the particle momenta
P n by half a time step to x. Third, 𝑷𝑛 is advanced with the help of the ﬁelds 𝑬𝑛+1/2 and𝑩𝑛+1/2 and the information 𝒙𝑛+1/2 to 𝑷𝑛+1 . The current 𝑱𝑛+1 is calculated. Finally, 𝑬𝑛+1/2
and 𝑩𝑛+1/2 are advanced to 𝑬𝑛+1 and 𝑩𝑛+1 .

7Chapter 2

Figure 2.1: Calculation flow in a single time step: the time step changes n step to n+1 step. First, the electric ﬁeld
E n and the magnetic field B n are advanced with the help of the current J n by half a time step to E n+1/2 and B n+1/2.
Second, the particle position x n is advanced with the help of the particle momenta P n by half a time step to x.
Third, 𝑷𝑛 is advanced with the help of the ﬁelds 𝑬𝑛+1/2 and 𝑩𝑛+1/2 and the information 𝒙𝑛+1/2 to 𝑷𝑛+1 .
The current 𝑱𝑛+1 is calculated. Finally, 𝑬𝑛+1/2 and 𝑩𝑛+1/2 are advanced to 𝑬𝑛+1 and 𝑩𝑛+1 .

2-3 Calculation of the fields
2-3-1 Maxwell’s equations
The change of the magnetic and electric field confirms to the Maxwell’s differential equations
as follows;
,

(2-3)

,

(2-4)

,

(2-5)

where 0 is the magnetic permeability of vacuum, 0 is the dielectric constant of vacuum, and  isthe charge density. Eq. (2-3) is rewritten to Eq. (2-6) using the relation c = (00)-1/2, here c is the
speed of light, as follows:.
In a numerical analysis, the magnetic field and the electric field are advanced by solving thedifference equations of Eqs. (2-4), (2-5), and (2-6).

8

(2-6)

Chapter 2
2-3-2 Finite-Difference-Time-Domain (FDTD) scheme
As shown in Eqs. (2-4) and (2-6), the change of electric field in time is dependent on therotation of magnetic field and also the change of magnetic field in time is dependent on therotation of electric field. Therefore, the field values cannot be solved in the same time. The FiniteDifference-Time-Domain (FDTD) scheme [25] is a method for solving the Maxwell’s differential
equations to stagger the position of field values and the time step for each values. The distributed
grid is called a Yee grid [25]. For an electric field component, it moves a half grid point in the field
points. For a magnetic field component, it moves a grid point in the field points. The twodimensional Yee grid is shown in Fig. 2.2.

Figure 2.2: Two dimensional Yee grid. For an electric field component, it moves the field half agrid point in the rotational direction. For a magnetic field component, it moves the field half a grid
point in all directions except the one it points.

The Maxwell’s differential equations of Eqs. (2-4), (2-5), and (2-6) are indicated as follows:
𝑑𝐸𝑥
𝑑𝑡
𝑑𝐸𝑦
𝑑𝑡
𝑑𝐸𝑧
𝑑𝑡
𝑑𝐵𝑥
𝑑𝑡
𝑑𝐵𝑦
𝑑𝑡
𝑑𝐵𝑧
𝑑𝑡

= 𝑐2 (

𝑑𝐵𝑧
𝑑𝑦
𝑑𝐵𝑥

𝑑𝐵𝑦

−

𝑑𝑧
𝑑𝐵𝑧

𝐽

)− 𝑥

𝜀0
𝐽𝑦

= 𝑐 2 ( 𝑑𝑧 − 𝑑𝑥 ) − 𝜀 ,

(2-7)

0

𝑑𝐵

𝑑𝐵

𝐽

= 𝑐 2 ( 𝑑𝑥𝑦 − 𝑑𝑦𝑥 ) − 𝜀𝑧

0

𝑑𝐸

𝑑𝐸

𝑑𝐸

𝑑𝐸

𝑑𝐸

𝑑𝐸

= − ( 𝑑𝑦𝑧 − 𝑑𝑧𝑦 )= − ( 𝑑𝑧𝑥 − 𝑑𝑥𝑧 ) ,(2-8)= − ( 𝑑𝑥𝑦 − 𝑑𝑦𝑥)
.

(2-9)

Next, the discretization of spatial derivatives is introduced as follows:,
9(2-10)

Chapter 2,
where(2-11). The step sizes along the coordinate directions and time aredenoted by ∆x, ∆y, ∆z, and ∆t. The ﬁeld F represents each of the ﬁeld components of E and B
with the properties,.

(2-12)The discrete Maxwell equations are written as follows:,(2-13),(2-14),(2-15)) ,(2-16)

where
𝑛+

1

1
2
1
𝑗+ 𝑘𝑙
2

𝑛+

𝑬𝑗𝑘𝑙 2 = ((𝐸𝑥 )

1
2
1
𝑗𝑘+ 𝑙
2

, (𝐸𝑦 )𝑛+𝑛+

, (𝐸𝑧 )

1
2

𝑗𝑘𝑙+

1
2

𝑛

𝑛
𝑛
𝑛
𝑩𝑗𝑘𝑙
= ((𝐵𝑥 )𝑗𝑘+
1 1 , (𝐵𝑦 )1 ) ,1
1 , (𝐵𝑧 ) 1𝑙+𝑗+ 𝑘+ 𝑙𝑗+ 𝑘𝑙+
2

2

2

2

2

(2-17)

2

𝑛+1

𝑛+1
𝑛+1
𝑱𝑛+1
𝑗𝑘𝑙 = ((𝐽𝑥 )𝑗+1𝑘𝑙 , (𝐽𝑦 )𝑗𝑘+1𝑙 , (𝐽𝑧 )𝑗𝑘𝑙+1 ) .
2(2-18)

2

2

In the equation of motion, the ﬁelds at half time steps are needed. Hence, Maxwell’s equations
are split as follows:,(2-19),(2-20),(2-21),(2-22).

(2-23)

When we use a two dimensional model, the z derivative was omitted as follows:
𝑑𝐸𝑥
𝑑𝑡
𝑑𝐸𝑦

𝑑𝐵

𝐽

= 𝑐 2 𝑑𝑦𝑧 − 𝜀𝑥0𝐽𝑑𝐵
= −𝑐 2 𝑑𝑥𝑧 − 𝜀𝑦
𝑑𝑡
0
𝑑𝐸𝑧
𝑑𝐵𝑥
2 𝑑𝐵𝑦

𝐽

𝑑𝑡

0

,

= 𝑐 ( 𝑑𝑥 − 𝑑𝑦 ) − 𝜀𝑧10(2-24)

Chapter 2
𝑑𝐵𝑥

𝑑𝑡
𝑑𝐵𝑧
𝑑𝑡

𝑑𝐸

= − 𝑑𝑦𝑧

𝑑𝑡
𝑑𝐵𝑦

=

𝑑𝐸𝑧

𝑑𝑥
𝑑𝐸𝑦
𝑑𝐸𝑥

,

(2-25)= − ( 𝑑𝑥 − 𝑑𝑦 ).

(2-26)

2-3-3 Boundary condition
The boundary conditions for the magnetic and electric field are to clamp the field given by thefield variable to zero on the boundary. For example, when the value of Ex is clamped to be zero on
the boundary, then just set Ex(0) = 0 since Ex(0) lies on the boundary as shown in Fig. 2.2. Bx(0) =- Bx(1) has to be set to do the same for Bx on the boundary. This is because the point between
Bx(0) and Bx(1) has the value Bx(1/2) = (Bx(1) + Bx(0))/2 = 0 using a linear reconstruction of Bx.
The electro-magnetic waves outﬂowing characteristics propagate through the boundary on thiscondition.
2-3-4 Initial calculation for preparing a background magnetic field
Setting of the background magnetic field was needed to simulate a magnetic nozzle. In thisstudy, two currents Jz in opposite directions were used for setting the dipole magnetic field asshown in Fig. 2.3. The direction of the magnetic field produced by the current is determined by theright-hand screw rule. In addition, due to the character of the EPOCH code, an electric field was
generated according to the change of the magnetic field, so each current was set at first, then thesolving particles motion was started in a state where the magnetic field becomes stable.

Figure 2.3: Initial magnetic field topology. Two currents in opposite directions (red) is used forsetting the dipole magnetic field. The black arrows show the magnetic flux line. The calculation of
initial magnetic field is executed before the main calculation starts. This magnetic field is obtained
the calculation of Eqs. (2-19) - (2-23), 40000 cycles.

11

Chapter 2

2-4 Calculation of particles
2-4-1 Particle shape function [26]The particle shape function is important for controlling the smoothness of the solution of a
PIC simulation since the number of particles in an actual plasma is 1024 m-3 in the experimental
plasma case, CPU cannot handle all of the particles. The function describes the assumed
distribution of the real particles to replace a super particle. The simplest method is to assume thatthe super particles uniformly spread in the cell where the super particle is located. This has theadvantages of calculation speed and simplicity, however also has the disadvantages of inducing
very large numerical error. The next simplest method is to assume a triangular shape function withthe peak of the triangle located at the position of the super particle and a width of 2∆x. Figure 2.4
shows the overlap between the distributed force and the form factor of a super particle with atriangular shape. The force distribution is indicated by the rectangle and the super particle by thetriangle. The x j denotes the grid at which the electric and magnetic fields are defined. The positionof the quasi-particle is xi. The electric and magnetic fields reside on staggered grids in case the
FDTD scheme is used. The bold line indicates the cell on the field grid that the particle is in. The
thick dots define the grid that is used to localize the particles. This method is a good trade-oﬀ
between cleanness of solution and calculation speed, so that the EPOCH code use a triangularshape function.

xi

xj - 2

xj - 1

xj

xj+1

xj+2

xj+3

Adapted from Hartmut Ruhl [26]Figure 2.4: The overlap between the distributed force and the form factor of a super particle with atriangular shape. The force distribution is indicated by the rectangle and the super particle by thetriangle.

The plasma behavior is obtained from solving the Vlasov equation. In the PIC simulation, thesolution is obtained by propagating super particles according to the equation of motion Eqs. (2-1)
and (2-2). For this, it is necessary to set the force distribution around the calculation grids. The
force distribution around a calculation grid x j is assumed as follows:(2-27)In the following, xj -∆x/2  x  xj +∆x/2 is assumed. Under this provision, the shape function of thesuper particle at xj overlaps with the forces distributed around the grids x j-1, x j, and xj+1 as shown in
Fig.2.4. The integral is split into four parts: a part of the shape function which overlaps with thecell x j-1, a part of the shape function from the left boundary of xi to the point of the triangle, a partof the shape function from the point of the triangle to the right hand edge of xi, and ﬁnally a part ofthe shape function which overlaps cell xj+1. Hence, the acceleration of the super particle is shownas follows:
12

Chapter 2.

(2-28)

Performing the integrations,+.

(2-29)

Eq. (2-29) shows that a particular interpolation scheme for the force equation has to be applied forthe shape function of super particles.
2-4-2 Equation of motion [26]The equation of motion Eq. (2-1) is rewritten as follows:
𝑑𝑃𝑥
𝑑𝑡
𝑑𝑃𝑦
𝑑𝑡
𝑑𝑃𝑧
𝑑𝑡

= 𝑞(𝐸𝑥 + 𝑉𝑦 𝐵𝑧 − 𝑉𝑧 𝐵𝑦 ),= 𝑞(𝐸𝑦 + 𝑉𝑧 𝐵𝑥 − 𝑉𝑥 𝐵𝑧 ),= 𝑞(𝐸𝑧 + 𝑉𝑥 𝐵𝑦 − 𝑉𝑦 𝐵𝑥 ).

(2-30)The integration of Eq. (2-2) up to second order accuracy is obtained by first integrating for half atime step up to first order accuracy as follows:,,.

(2-31)In addition, the momentum integration is performed up to second order accuracy as follows:(2-32)
Finally, Eq. (2-2) is integrated again over half a time step with first order accuracy. However, thevelocities V(t +∆t/2) are known yielding overall integration up to second order accuracy as
follows:,,
13

Chapter 2.

(2-33)

It is convenient to introduce the following variables [27] as follows:,,(2-34)

where Px/y/z denotes the momenta Px, Py, and Pz. A similar convention applies to the electric field.
Inserting Eq. (2-34) into Eq. (2-32) yields,,
𝑃𝑦− − 𝑃𝑦+ = −𝑞∁𝐸𝑉𝑧 𝐵𝑥 + 𝑞∁𝐸𝑉𝑥 𝐵𝑧,,(2-35)

where it is found up first order accuracy in ∆t,
𝑉𝑥⁄𝑦 ⁄𝑧 =𝑃𝑥+⁄𝑦⁄𝑧 +𝑃𝑥−⁄𝑦⁄𝑧
2
2
− 2
𝑃+
𝑃+ +𝑃−
𝑃+ +𝑃−𝑦 +𝑃𝑦𝑥) +() +( 𝑧 𝑧 )
2
2
2

.

(2-36)

2√1+( 𝑥Thus it follows that,,,

(2-35).

(2-36)

where τx, τy, and τz are defined as

Solving for 𝑷+
𝑥 yields
𝑃𝑥+
𝑃𝑥−+(𝑃𝑦 ) = 𝑨 (𝑃𝑦− ),𝑃𝑧−𝑃𝑧+(2-37)

where.

Eq. (2-38) represents a rotation. Therefore, it is found up to first order accuracy as follows:

14(2-38)

Chapter 2
𝑃 + +𝑃𝑥−(𝑥22

𝑃 + +𝑃𝑦−) +( 𝑦22

𝑃 + +𝑃𝑧−) +( 𝑧22) = 𝑃𝑥− 2 + 𝑃𝑦− 2 + 𝑃𝑧− 2 .

(2-39)

2-4-3 Boundary condition
The boundary condition for super particles is an open boundary. Super particles are transmitted
through the boundary and removed from the calculation without a reflection. Therefore, the energy
conservation of particles in the total calculation region is not assured after any particles are
removed from the calculation region.

2-5 Calculation of the current
EPOCH uses the Villasenor and Buneman current calculating scheme [28] which solves Eq.
(2-5) to calculate the current at each time step. The main advantage of this scheme is that itconserves a charge on the grid rather than just globally conserving charge on the particles. Thismeans that the error in the solution of Poisson’s equation is conserved, so if Poisson’s equation issatisﬁed for t = 0, it remains satisﬁed for all time. The Villasenor and Buneman scheme works
because exactly the same charge added to one cell is subtracted from another cell, which in turnmeans that exactly the same current added to one cell is subtracted from another cell. This iscorrect since a particle crossing a cell boundary would represent the loss of that particle’s
contribution to the current from the source cell and the gain of that particle’s contribution to thecurrent by the destination cell. In fact, this simple type of cell boundary crossing current
calculation was used in classical Buneman type PIC codes. The scheme is simple. After the main
particle push, the particle is advanced a further half time step into the future to ﬁrst order using thevelocities calculated at the end of the particle push. The particle positions at t +∆t/2 are stored
earlier, and combined with the newly calculated particle position at t +3∆t/2. This allows a time
centered evaluation of d/dt, meaning that the current update is second order accurate in time. The
spatial order of the scheme matches the spatial order of the particle shape function. The shapefunctions for transferring particle properties onto the grid at the two time steps are calculated
including a shift when necessary to allow for the particle having crossed a cell boundary. Since thecharge associated with the particle is spatially distributed using the shape function, all that isnecessary to calculate d/dt is to subtract the two functions, multiply by the charge on the superparticle and the super particle weight and ﬁnally divide by dt. The spatial derivative of J is thenconverted to a one sided ﬁnite diﬀerence form and solved directly. In multiple dimensions, this isslightly complicated by the eﬀects of oﬀsets in directions other than the direction that a given
current component is pointing in, with this adding additional weight factors based on the overlapof the shape functions in other directions.

2-6 Super particle
A laser produced plasma contains over 1024 particles per m3, and handling all of this particle
is impossible even with the latest supercomputer, which is not practical. Therefore, super particles
are used as particles to be handled in the mass, while keeping the charge mass ratio of the actual
particles constant while summarizing the charge and mass of many particles into one in particle
simulation.

2-7 Time step and grid length
The time step corresponds to the plasma frequency. If it is less than 0.1

ωpe, a numerical instability does not occur. Also, the Courant-Friedrichs-Lewy (CFL) condition is
15

Chapter 2
needed to satisfy. In the CFL condition, the time step is decided to satisfy dt = dx/c. In thissimulation, the time step is decided to the minimum value to satisfy the both condition. In order tosimulate the phenomenon associated with the cyclotron motion of charged particles in a magnetic
field, the time step is needed to be under electron cyclotron frequency. This condition is alreadysatisfied because the electron cyclotron frequency is smaller than the plasma frequency when ne =1024 m-3. The grids length for calculation is decided to be the value under the initial rc.

2-8 Parallel calculation method
A difficult problem of PIC simulation is a large calculation cost. The largest part of thecalculation cost is the particle calculation due to use over 1 million particles. There are generally
two methods to reduce the calculation cost. The first method is a particle division parallelization
method: it allocates each particles calculation to many CPU core depending on the total number ofparticles. The second method is a spatial division parallelization method: it allocates the divisional
calculation regions to many CPU core depending on the total parallelization processes. In our
simulation model, the initial super particles are concentrated to the near origin of calculation
region because we assume an isotropic expansion plasma initially. Hence, the spatial divisionparallelization method is suitable to reduce the calculation cost using Message Passage
Interface(MPI).

2-9 Parallelization in EPOCH code
EPOCH code uses a load balancing [29] for the limiting factors in the spatial divisionparallelization method of a PIC simulation. Due to the synchronization of the currents required forthe update of the electromagnetic fields, the code runs at the speed of the slowest process. The
load balancer is used to move the inter-processor boundaries, so that the number of particles is as
close to the same on each processor as possible. The load balancer is invoked at the start of thecode and when the ratio of the least loaded processor to the most loaded processor falls below a
specified critical point. The load balancer in EPOCH code works by rearranging the processor
boundaries in one dimensional sweeps in each direction, rather than attempting to perform
multidimensional optimization. Also, the MPI in EPOCH code requires each processor to be
simply connected at every point, so it must have one processor to the left, front, and right,
respectively. That introduces a further restriction on the load balancer. Otherwise, the loadbalancer is fully general. The load balancing sweep is shown in Fig. 2.5. The specified parameteris used to force the code to perform a load balancing sweep even when it would normallydetermine that the imbalance is not large enough to force a load balancing sweep.

Figure 2.5: The load balancing sweep. The specified parameter is used to force the code toperform a load balancing sweep even when it would normally determine that the imbalance is notlarge enough to force a load balancing sweep.
16Chapter 3

Chapter 3: Test Calculation Results
3-1 Calculation model
The two dimensional calculation model in this study is shown in Fig. 3.1. The origin of
calculation region is located on the left-side of the region to focus on a downstream of a magnetic
nozzle (right-side). The initial ions and electrons are distributed in a plasma radius Rp (the red circle
as shown in Fig. 3.1). The initial plasma velocity is expressed as V = |rV0/Rp|r, where r is √𝑥 2 + 𝑦 2 ,V0 is the maximum velocity of plasma. Here, we assume that the initial plasma expands isotopically.
The absolute value of the initial plasma velocity |V| is set to correspond to the initial particle position
since the edge particles have higher velocity than other particles that located near points ofcalculation origin, when the initial distribution of plasma has expanded isotopically. The coil current
is set on two grids to produce the dipole magnetic field. Particle collision with the coil is notconsidered in this simulation.

Figure 3.1: Two dimensional calculation model. The initial plasma distribution is decided to fill thecircle(red). The opposing coil currents are located on the two grids to make a dipole magnetic field.

17

Chapter 3

3-2 Verification of EPOCH code
First, a simple calculation was executed for a verification of EPOCH code. The calculation
condition was set as follows. The number of grids was set to 500 in x and y directions, respectively(totally 500 500). The grid widths ∆x and ∆y were set to 5 m. The initial plasma expansion
radius, Rp were set to 10∆x. The initial maximum velocity of particles, V0 was set to 105 m/s toradial direction isotopically. The total number of super particles was set to 10 6. The number
density of the initial plasma, np was 1024 m-3. The mass ratio of ion and electron, mi/me was set to
100 to reduce the difference of moving time scale that change the cyclotron frequency of ion andelectron since ion almost stay in initial positions when mi/me is real ratio in a numericalcalculation. The time step ∆t was set to 10-14 s. The magnetic field was set to Bx = 0 or 1 T as aparallel magnetic field to verify the effect of magnetic field.
Figures 3.2 and 3.3 show the calculation results with B = 0 and 1 T, respectively. Fig. 3.2
shows the results of Bx =0 T. The spatial distribution of number density in Fig. 3.2 shows that bothions and electrons expand isotopically keeping initial velocities. The expansion radius changes
from 50 m to 60 m at 0.1 ns which is comparable to the free-streaming velocity V0t =105 m/s 
0.1 ns = 10 m (∆Rp). The spatial distribution of x component of electric field shows that local
electric field is created at the edge of plasma because some electron escape from plasma cloud andthose make the potential difference at the edge. This means that the difference of mass between
ion and electron works correctly. Figure 3.3 shows the results of Bx = 1 T. The spatial distributionof number density in Fig. 3.3 shows that ions and electrons keep expanding isotopically until 0.1
ns. Only electrons drift to x direction at 0.2 ns because of cyclotron motion. Ex distribution showsthat electric field is created at the right edge of the plasma (x ~50 - 100 m and y ~-50 – 50 m).
The Ex generated by the plasma is bipolar seeing from macro-viewpoint.
Two cases Figs. 3.2 and 3.3 show different field structures. This difference is from themagnetic field setting: no external magnetic field and parallel magnetic field along x direction.
The EPOCH code indicated that it can calculate the electron cyclotron motion in a magnetic fieldcorrectly.

18

Chapter 3

Figure 3.2: The spatial distribution of number density of ions and electrons, and x component ofelectric field with B = 0 T, respectively. Ions and electrons keep expansion isotopically. Ex
distribution shows that the electric field is created at the edge of the plasma because of thepotential difference between ions and electrons.

19Chapter 3

Figure 3.3: The spatial distribution of number density of ions and electrons, and x component ofelectric field with B = 1 T, respectively. Ions and electrons keep expansion isotopically until 0.1
ns. Electrons drift to x direction at 0.2 ns because of the cyclotron motion. Ex distribution showsthat electric field is created on the x edge of plasma because of the potential difference between
ions and electrons.

20Chapter 3

3-3 Validation of a magnetic nozzle working on EPOCH code
Next, a test calculation was executed to validate a calculation of a magnetic nozzle using
EPOCH code. Figures 3.4–6 show plasma number density, magnetic flux density, and JB force
obtained from EPOCH calculation, respectively. The conditions of simulation were decided fromprevious research [14]. The number of grids was set to 500 on each x and y direction (totally 500
500). The coil radius (Rc, the distance from a coil current grid to the coil center) was set to 30∆x,where ∆x is grid width = c/ωpe ~ 5 m. The distance between the coil center to the center of initialplasma distribution RL, was set to 3RC that is same condition as previous research [14]. The initialplasma expansion radius Rp, was set to 0.42 RC and the initial particle drift velocity V0, was set to1.41105 m/s in the same way [14] . The total number of super particles was set to 106 and thenumerical weight for each a super particle was set to 1019 since the number density of initial plasma
np, is 1024 m-3. The mass ratio of ion and electron mi/me, was set to 100. The time step ∆t was set to
0.1/fpe ~10-14 s. The magnetic flux density at the center of initial plasma B0, was set to 0.1 T. The
cyclotron radius of electron with B ~0.1 T has same scale as ∆x and the cyclotron radius of ion has
100 times larger scale than electron. Therefore, ions are weakly magnetized and electrons aremagnetized, so that the plasma is expected to receive the Lorentz force induced from the magnetic
field in the magnetic nozzle.
Figures 3.4(a) and 3.5(a) show the spatial distribution of plasma number density and absolutevalue of magnetic flux density at t = 0.1 ns. Rp is 80 m on the time, so that the plasma expands 15
m from the initial position. This means that the plasma expansion is correctly calculated because
V0t ~ 15 m. Also, the magnetic flux density is decreased from 0.1 T to 0.08 T in the plasma region.
This means that the diamagnetism of the plasma is correctly calculated and the plasma countervails
the initial magnetic field. Therefore, a magnetic field cavity is produced where plasma locates asshown in Fig. 3.5(a). Figures 3.4(b) and 3.5(b) show the number density and magnetic flux density,
respectively at t = 3 ns. The plasma distribution is asymmetry in the y axis because the plasma that
initially locates in x > 0 has positive Vx keeping free expansion. On the other hand, the plasma that
initially locates in x < 0 has negative Vx. It obtains Lorentz force to be decelerated. The left edge ofplasma is decelerated by the force, then it is piled up. As the results the number density of left edge
is increased as shown in Fig. 3.4(b). The plasma cavity grows up to form a funnel like structuresince the plasma expand without the deceleration when those only have a parallel velocity to themagnetic field line (x ~ -350– -200 m and y ~ -30–30 m) as shown in Fig. 3.5(b). Figures 3.4(c),3.4(d), 3.5(c), and 3.5(d) show the distributions at t = 6 ns and 9 ns, respectively. The plasma is notcompletely decelerated in -x direction and it escapes from the calculation region being guided by
the magnetic field as shown in Figs. 3.4(c) and 3.5(c). Later in time, high density plasma locatesnear the origin as shown in Fig. 3.4(d). This density plasma consists of slow particles which initiallylocate near origin. The magnetic flux density is increased between x = -400 to -350 m at 9 ns after
the particle escaping because the magnetic cavity becomes shallow by the decrease of plasma
number density as shown in Fig. 3.5(d).
Figure 3.6 shows the spatial distribution of Lorenz force. Figures 3.6(a) and (b) show the xcomponent of Lorentz force: (JB)x = (JyBz - JzBy) [Am-2 T = Nm-3] at t = 3 ns and 6 ns. Figures
3.6(c) and (d) show the y component of Lorentz force; (JB)y = (JzBx - JxBz) Nm-3 at 3 ns and 6 ns.
It is confirmed that a deceleration force in x direction is large on the plasma edge ( x,y = -200, ∓200m) as shown in Fig. 3.6(a), and also the force in y direction is large in x < -300 m as shown in
Fig. 3.6(c).The magnetic field diverges as x increase and this is why the plasma expands in ydirection, different from a parallel magnetic field as shown in Figs. 3.6(b) and (d). These calculation
results showed the Lorentz force is correctly calculated in the EPOCH code.
21

Chapter 3

Figure 3.4: The spatial distribution of plasma number density plasma np m-3. The sequence of (a)–(d) shows the changing over time. The snapshot time are at 0.1 ns(a), 3 ns(b), 6 ns(c), and 9 ns(d).

Figure 3.5: The spatial distribution of magnetic flux density. The contour shows the absolute value|B| T. The sequence of (a)–(d) shows the changing over time. The snapshot time are at 0.1 ns(a), 3
ns(b), 6 ns(c), and 9 ns(d).
22Chapter 3

Figure 3.6: The spatial distribution of JB force. The contours (a) and (b) show the x component,(c) and (d) show the y component of the force. The snapshot time are at 3 ns (a) and(c), at 6 ns (c)and (d).

3-4 Investigation of calculation time and use of parallelization
A huge number of particles are required in particle simulation, in general, because enoughnumber of particles are needed in each cell to calculate electric and magnetic field accurately. As
the number of particles is increased, the total calculation cost is larger. To estimate the calculation
time, we have checked the dump time of calculation that depends on the number of particles. The
calculation model is same as explained in the section 3-2, and the total number of particles waschanged from 106 to 107. As shown in Table 3.1, it takes six hours and two minutes calculating until
1 ns on the 106 case. On the other hand, it takes 57 hours and 25 minutes calculating until 1 ns on
the 107 case. As we expected, the calculation time is directly proportional to the total number ofparticles. This simulation dump is short to investigate the behavior of plasma in the magnetic nozzle.
We need five times longer time to simulate plasmas in downstream region, and the parallelizationis necessary to calculate in realistic time.
Table 3.1: Comparison of calculation time
Total number of particles
Simulation dump
Calculation machine dump
106
1 ns
6 h and 2 min
7
10
1 ns
57 h and 25 min
As explained in the section 2.9 in the chapter 2, the EPOCH code uses the load balancer to
optimize the parallelization. As shown in Fig. 3.1, the initial particles are located near the calculation
origin and most of the regions are handled as the vacuum. Figure 3.7 shows the division patternusing the load balance method. The initial plasma region is subdivided in 4 regions for plasma (4for total), 8 regions for plasma (8 for total), 30 regions for plasma (36 for total), and 49 regions for23

Chapter 3
plasma (for total 144, 288, and 576). It found that the maximum division number for plasma regionis 49. Figure 3.8 shows the comparison of calculation machine time. The calculation model wassame as explained in the section 3-2. The total number of particles was 106 and simulation dump
was set to 5 ns. Here, MPI process number 1 means no parallelization and the calculation machinetime is not directly 5 times of Table 3.1 because some particles escape from the calculation region
during the calculation and it makes the calculation time shorter. The calculation machine time using
144 processers is 43.4 % of that without parallelization. This means that the speedup in latency
Slatency ~ t no parallelization / t parallelization ~ 2.3. However, the expected calculation time on 288 and 576
processers cases are not decreased, because the division balance is determined by the initial particle
distribution and the division number of initial plasma region is same in the cases of 144, 288, 576
processers. This Slatency is also expressed with a law of Amdahl [30]:𝑆letency =1𝑃𝑁
where P is the proportion of parallelization in the calculation code, N is the number of processers
used. Therefore, the proportion P in our case is estimated as 0.567. In each region, the iterations ofparticle calculations are not parallelized in our code. To increase P, the parallelization in practicecalculation would be required. Hence, N = 144 was the optimum number for calculating the plasma
behavior in the magnetic nozzle in the present simulation.
(1 − 𝑃) +24

Chapter 3

4 processers

8 processers

36 processers

144 processers

288 processers
576 processers
Figure 3.7: Allocation of divided regions on CPU cores using spatial division method on EPOCH
code. The division pattern is determined by the load balance method. The initial plasma region is
subdivided in 4 regions for plasma (4 for total), 8 regions for plasma (8 for total), 30 regions forplasma (36 for total), and 49 regions for plasma (for total 144, 288, and 576).
25

Chapter 3
1600
1, 1448

Calculation dump time [min]1400

36, 1284

1200
1000
800
600

288, 648

144, 629

576, 563

400
200

0
0

100

200

300

400

500

600

700

Spatial division number [MPI process]Figure 3.8: Comparison of calculation machine dump time. The total number of particles was 10 6
and simulation dump was 5 ns. The calculation machine time using 144 processers is 629 minutes.
It is the 43.4 % of no parallelization case. The decreasing effect is comparable in 144, 288 and 576
processers cases.

Figure 3.9: The relationship of number of processer and speedup in latency by parallelization [30].
For example, the theoretical limit of speedup in latency becomes 2 when the proportion ofparallelization on code is 0.5; Slatancy = 1/(0.5+(1-0.5)/N)).

26

Chapter 4

Chapter 4: Simulation results and discussion
The simulation results are discussed in this chapter. In the first part, the five cases with different magnetic
field flux density are explained. Then, the simulation results are shown and explained. After that, the behavior of theplasma in the simulation is discussed. Finally, the electron detachment in these conditions are discussed.

4-1 Calculation condition
Calculation model was same as explained in the sections 3-1 and 3-3 in the chapter 3. The
behavior of the plasma in the magnetic nozzle is changed by magnetic field strength, therefore,
five calculation cases with different magnetic flux densities were considered as shown in Table
4.1. These were defined from the criterion value b [31] :
𝜀b =

𝑟ci

,

𝑟𝑏

(4-1)

where rci is cyclotron radius of ion. rb is calculated as follows:1

𝑟b = (

3𝜇0 𝐸0 3),
2𝜋𝐵0 2(4-2)

where E0 is the plasma kinetic energy and B0 is the magnetic field density at plasma location.
When b  1, both of ion and electron in a plasma are magnetized, meaning that both particles
move with the magnetic field. In this study, b were from 1.3 to 6 in the initial conditions.

Table 4.1: Calculation condition
Number of calculation grids (x,y)
4000, 4000
∆x and ∆y510-7 m (< rce)∆t

10-15 s (< 0.01/fpe)

rC

300∆x

rL

3rC

rP

0.42rCnp

1024 m-3

V0

1.41105 m/s

E02.7310-13 J

B00.014, 0.14, 0.35, 0.7, 1.4 T

b6, 2.8, 2, 1.6, 1.3

Number of particles

106

mi/me

100

27

Chapter 4

4-2 Time evolution of plasma number density
Figures 4.14.5 show the time evolution of plasma number density in the case of B0 = 0.014
T, 0.14 T, 0.35 T, 0.7 T, and 1.4 T at t = 0.52 ns. As shown in Fig. 4.1, the plasma expands
isotopically to the direction of initial velocity in the weak magnetic field of 0.014 T. The number
density decreases with passage of time because the weak magnetic field hardly generates anacceleration force to the plasma. A high density plasma locates on the center of calculation region.
This plasma does not have any expansion velocity by the initial condition: V0 = |rV0/Rp|r ~ 0, so it
stays the region. Therefore, the density cavity is created in 0.1 mm < (x2+y2)1/2 < 0.2 mm at t = 2ns.

0.5 ns

1 ns

1.5 ns2 ns

Figure 4.1: The number density of plasma with B0 = 0.014 T and b = 6 at 0.5 ns, 1ns, 1.5 ns, and
2.0 ns.

28

Chapter 4
Figure 4.2 shows ten times stronger magnetic field than 0.014 T as shown in Fig.4.1. The
density increases near the coil currents at t = 1 ns since the plasma receives JB force as explainedin the section 3.3. As the result, the plasma is reflected by the magnetic field near the coil. The
front of the reflected flow expands on x = 0.5 mm at t = 2 ns.

0.5 ns

1 ns

1.5 ns2 ns

Figure 4.2: The number density of plasma with B0 = 0.14 T and b = 2.8 at 0.5 ns, 1ns, 1.5 ns, and
2.0 ns.

29

Chapter 4

Figures 4.3 and 4.4 show 1.5 and 5 times stronger magnetic field case than Fig. 4.2: 0.35 T, and
0.7 T. In these cases, the plasma reflection points become far from the coil currents since themagnetic flux density becomes larger. Therefore, the front of the reflected flow expands on x = 0.8mm at t = 2 ns as shown in Fig. 4.3 and on x = 1.1 mm at t = 2 ns as shown in Fig.4.4.
0.5 ns

1 ns

1.5 ns2 ns

Figure 4.3: The number density of plasma with B0 = 0.35 T and b = 2 at 0.5 ns, 1ns, 1.5 ns, and
2.0 ns.

30

Chapter 4

0.5 ns

1 ns

1.5 ns2 ns

Figure 4.4: The number density of plasma with B0 = 0.7 T and b = 1.6 at 0.5 ns, 1ns, 1.5 ns, and
2.0 ns.

31

Chapter 4

Figure 4.5 shows strongest magnetic field of 1.4 T in this simulation. The plasma is restricted
to expand to y direction at t = 0.5 ns because of strong magnetic field. The plasma can expand to y
direction in the downstream where the magnetic field is weak.

0.5 ns

1 ns

1.5 ns2 ns

Figure 4.5: The number density of plasma with B0 = 1.4 T and b = 1.3 at 0.5 ns, 1ns, 1.5 ns, and
2.0 ns.

32

Chapter 4

4-3 Evaluation of total plasma momentum
Figure 4.7 shows the density of (JB)x force Nm-3 with B0 = 1.4 T. This force accelerates theplasma in x direction. The density becomes high in the plasma edge ( -0.2 mm < y < 0.2 mm) at 0.5
ns because the gradient of plasma density becomes high in the region. (JB)x decreases as time
passes from 1 ns to 2 ns because the plasma moves to downstream, and the gradient of plasma
density becomes low, so the diamagnetic current in the plasma becomes small.

0.5 ns

1 ns

1.5 ns2 ns

Figure 4.6: The density of (JB)x force Nm-3 with B0 = 1.4 T and b = 1.3 at 0.5 ns, 1ns, 1.5 ns, and
2.0 ns.

33

Chapter 4
Figure 4.7 shows the time evolution of the total (JB)x. The peaks are at 0.7 ns 0.9 ns with
0.14 T 1.4 T as shown in Fig.4.7. Therefore, the acceleration by (JB)x continues to the peak inthis simulation of magnetic nozzle.

Figure 4.7: The time evolution of total JB force in x direction.
Figures 4.8 shows the total of the plasma momentum in x direction:  mVx. The y and z
components of plasma momentum are almost zero because the plasma expands symmetricallyalong y axis as shown in Figs 4.14.5.  mVx is almost zero until t = 0.3 ns because the plasma
expands isotopically. After t = 0.3 ns, the plasma which approaches to the coil current receives
the JB force, therefore, mVx increases. After the peak at t ~ 1.3 ns, the  mVx decreases sincesome particles arrive at the boundary of calculation region and we ignored the momentum ofthem. The maximum momentum is obtained in the case of 1.4 T since the magnetic field energy islarge and it produces the large JB force.

Figure 4.8: The time evolution of total plasma momentum in x direction. The maximum
momentum is obtained in the case of the initial magnetic field of 1.4 T.
34

Chapter 4

We can judge that the magnetic nozzle in this simulation can obtain thrust when the plasma
detaches from the magnetic field. As explained in the section 4.1, both of ion and electron in aplasma are magnetized when b  1. On the other hand, ions are weakly magnetized or unmagnetized and electrons are magnetized in b > 1. Therefore, we need to analyze the b toinvestigate the plasma condition.

4-4 Analysis about the cyclotron radius
4-4-1 Estimation of thermal velocity
We analyze the ion and electron cyclotron radiuses statistically in order to investigate aplasma detachment condition. The cyclotron radius is calculated from the temperature of theplasma perpendicular to the magnetic field. The temperature should be calculated from thevelocity distribution in a single cell. However, the temperature was not defined precisely because
the total number of particles in a cell was less than 20. Therefore, we evaluated the velocity
distribution perpendicular to magnetic flux line in 50 50 cells. The velocity distribution waspotted as a histogram of V⊥B particle, where V⊥B particle is calculated as follows:
.

(4-3)

After that, assuming that the velocity distribution is based on Maxwell’s velocity distribution, theGaussian fitting was operated using follow equation,
2

2

𝐹(𝑣) = 𝑎𝑒 (−(𝑉−𝑉flow ) /(𝑉thermal ⁄√2) ) ,(4-4)
where a, Vflow, and Vthermal are determined by the fitting calculation. For example, Fig. 4.9 shows
the distribution of ion velocity perpendicular to the magnetic field where the averaged region at-0.025 mm  x  0 mm and -0.025 mm  y  0 mm (50  50 cells) in B0 = 0.35 T at t = 1.7 ns. The
ion thermal velocity perpendicular to the magnetic field (Vion thermal⊥B) is 3.27105 m/s as shown in
Fig. 4.7. The temperature is calculated from the relation: Vthermal = (2kBTi/m)1/2, and Ti becomes
1115 eV when Vion thermal⊥B = 3.27105 m/s. The ion cyclotron radius is calculated using theequation: rci = Vion thermal /ωpi ~ 1.02×10-4×Ti1/2B−1 m [32] ~ 9.88 mm. The electron cyclotron radius
is calculated in the same way as rce = Velectron thermal /ωpe ~ 2.38×10-6×Te1/2B−1 m [32].

Figure 4.9: The distribution of ion velocity perpendicular to the magnetic field in -0.025 mm  x 
0 mm and -0.025 mm  y  0 mm. The purple points show the number of particles each V⊥. The
22green line shows the fitted result: 𝐹(𝑣) = 𝑎𝑒 (−(𝑉−𝑉flow ) /(𝑉thermal⁄√2) ) determined a, Vflow, and
Vthermal.
35

Chapter 44-4-2 Analysis results
Figures 4.104.14 show the time evolution of ion and electron cyclotron radii estimated on x
axis: -0.025 mm  y  0 mm in the case of B0 = 0.014 T, 0.14 T, 0.35 T, 0.7 T, and 1.4 T at t =0.52 ns. Both of ion and electron cyclotron radii become large as time passes in the downstream(x > 0 mm). The rate of change of ion and electron cyclotron radius shows qualitative agreement.
Ions are unmagnetized in x < -0.4 mm with B0  0.35 T because ri > 1 mm ~ plasma scale length
as shown in Figs. 4.104.12. On the other hand, ions are weakly magnetized in x < -0.4 mm
because ri < 1 mm as shown in Figs. 4.134.14. In the other region, ions are unmagnetized.
Electrons are also categorized as same as the ion, the border is x = 0 mm.
Ion

Electron

Figure 4.10: Ion and electron cyclotron radii estimated on x axis in -0.025 mm  y  0 mm with B0= 0.014 T (b = 6).
Ion
Electron

Figure 4.11: Ion and electron cyclotron radii estimated on x axis in -0.025 mm  y  0 mm with B0= 0.14 T (b = 2.8).
Ion
Electron

Figure 4.12: Ion and electron cyclotron radii estimated on x axis in -0.025 mm  y  0 mm with B0= 0.35 T (b = 2).
36

Chapter 4

Ion

Electron

Figure 4.13: Ion and electron cyclotron radii estimated on x axis in -0.025 mm  y  0 mm with B0= 0.7 T (b = 1.6).
Ion

Electron

Figure 4.14: Ion and electron cyclotron radii estimated on x axis in -0.025 mm  y  0 mm with B0= 1.4 T (b = 1.3).

Figures 4.154.19 show the time evolution of ion and electron cyclotron radii in the case of
B0 = 0.014 T, 0.14 T, 0.35 T, 0.7 T, and 1.4 T at t = 0.52 ns. As shown in Fig. 4.15, both ion andelectron cyclotron radii exceed the plasma expansion scale length: L ~ 1 mm, so the plasma is unmagnetized in B0 = 0.014 T. The region of cavity of plasma in the center as explained in thesection 4.2 and the plasma expansion front were not analyzed because the number of superparticle is small in these regions. As shown in Figs. 4.1618, ion and electron cyclotron radii
exceed the plasma expansion scale length: L ~ 1 mm in downstream, in contrast, an electron
cyclotron radius is below L in near the coil current region: x = -0.45 mm, -0.15 mm < y < 0.15
mm. Here, electrons are magnetized as shown in Figs. 4.1618 at 12 ns. In addition, an electron
cyclotron radius is comparable to L in -0.5 mm < x < 0 mm as shown in Fig.4.16, -0.2 mm < x < 0mm as shown in Fig.4.17, and -0.1 mm < x < 0.1 mm as shown in Fig. 4.18. Electrons are weakly
magnetized there. Hence, these three conditions show that ions are unmagnetized except the coilregion, electrons change the magnetized condition as they move to downstream, therefore, three
regions are defined: magnetized (x < -0.2 mm), weakly magnetized, unmagnetized (x > 0.4 mm).
Fig. 4.19 shows an almost uniform distribution.

37

Chapter 4

Ion 0.5 ns

Electron 0.5 ns

1 ns

1 ns

1.5 ns1.5 ns

2 ns

2 ns

Figure 4.15: Ion and electron cyclotron radii with B0 = 0.014 T (b = 6).

38

Chapter 4
Ion 0.5 ns

Electron 0.5 ns

1 ns

1 ns

1.5 ns1.5 ns

2 ns

2 ns

Figure 4.16: Ion and electron cyclotron radii with B0 = 0.14 T (b = 2.8).

39

Chapter 4
Ion 0.5 ns

Electron 0.5 ns

1 ns

1 ns

1.5 ns1.5 ns

2 ns

2 ns

Figure 4.17: Ion and electron cyclotron radii with B0 = 0.35 T (b = 2).

40

Chapter 4

Ion 0.5 ns

Electron 0.5 ns

1 ns

1 ns

1.5 ns1.5 ns

2 ns

2 ns

Figure 4.18: Ion and electron cyclotron radii with B0 = 0.7 T (b = 1.6).
41

Chapter 4
Ion 0.5 ns

Electron 0.5 ns

1 ns

1 ns

1.5 ns1.5 ns

2 ns

2 ns

Figure 4.19: Ion and electron cyclotron radii with B0 = 1.4 T (b = 1.3).

42Chapter 4

4-5 Discussion
The expected detachment condition as stated in the section 1-4 and 1-5 was expressed asfollows:
|∇𝐵|

(4-5)

𝑟ce 𝐵 ≫ 1.

When Eq. (4-5) is satisfied, the electron may detach from the magnetic field. Figures 4.204.22
show the calculation results of rce|B|/B, here we used averaged B and |B| in 50  50 cells.
rce|B|/B exceeds the unity everywhere in the case of the low magnetic field of 0.014 T as shownin Fig. 4.20. rce|B|/B changes 0.3 located on x = 0.4 mm at t = 0.5 ns to 1 located on x = 0.8 mm
at t = 1 ns as shown in Fig. 4.21. In contrast, rce|B|/B is less than the unity in x < 0 mm in the caseof the large magnetic field of 1.4 T as shown in Fig. 4.22. This means that electrons detach frommagnetic field under the high b = 6 and attach to magnetic field under low b = 1.3.
Eq. (4-5) showed the possibility to explain the electron attachment and detachment. The
analysis shown here suggested that electron detachment can be discussed by using a simple Eq. (45). However, the accuracy of analysis became worse in the downstream because we could notanalyze the edge of downstream due to smaller number of particles in specific cells. Hence, theincrease of the total number of super particles is required to solve this.
The ions freely expands in the low magnetic field strength condition: b =6, but the plasma
structure changes in the high magnetic field strength condition: b =1.3, as shown in Figs. 4.14.5.
The ions detach from the magnetic field when b is low, and the electrons may follow that. The
ions weakly attach to the magnetic field in the downstream when b =1.3, and also the electron
strongly attach to the magnetic field. To confirm this supposition, particle trajectories would beanalyzed in future.

43

Chapter 4
B0 = 0.014 T
0.5 nsB0 = 0.14 T
0.5 ns

1 ns

1 ns

1.5 ns1.5 ns

2 ns

2 ns

Figure 4.20: Calculated rce|B|/B under B0 = 0.014 T and b = 6 (left panels), and B0 = 0.14 T andb = 2.8(right panels).
44

Chapter 4
B0 = 0.35 T
0.5 nsB0 = 0.7 T
0.5 ns

1 ns

1 ns

1.5 ns1.5 ns

2 ns

2 ns

Figure 4.21: Calculated rce|B|/B under B0 = 0.35 T and b = 2 (left panels), and B0 = 0.7 T and b= 1.6(right panels).
45

Chapter 4
B0 = 1.4 T
0.5 ns

1 ns

1.5 ns2 ns

Figure 4.22: Calculated rce|B|/B under B0 = 1.4 T and b = 1.3

46Chapter 5

Chapter 5: Conclusion
In this research, a full Particle-In-Cell (PIC) simulation was conducted to investigate theplasma behavior in the magnetic nozzle. A PIC method simulated the motion of charged particles
in plasma: ion and electron, and the equation of motion solved these. Also this method handled fieldvalues on grid points: magnetic field, electric field, and plasma current. the Maxwell’s equations
were solved by FDTD method.
The problem was the high calculation cost of particle calculation. Therefore, we have used aparallel calculation method and confirmed the reduction of the calculation cost. The method ofparallelization was explained in chapter 2, and the effect of parallelization was stated in chapter 3.
We have confirmed that the maximum number of parallelization for the calculation of magnetic
nozzle.
We have analyzed plasma flow and confirmed a difference caused by the conditions of adifferent ratio of b = rci / rb where rci is cyclotron radius of ion, and rb is length derived from theratio of the plasma kinetic energy to the density of the magnetic field energy in the section 4-1 inthe chapter 4., We have analyzed the cyclotron radius of ion and electron statistically, the result
showed the magnetized parameters of ions and electrons in the section 4-4.
Finally, we have discussed about the evaluation method of electron detachment to use therelationship between the cyclotron radius and the magnetic scale length: B/|B|. rce|B|/B
exceeded the order of unity in the case of the low magnetic field. In contrast, rce|B|/B was lessthan the order of unity in the case of the high magnetic field density. We have found that electronsdetach from magnetic field in the downstream and attach to magnetic field in near the coil.
However, the analysis is insufficient to evaluate the region where |B| changes locally such as aplasma cavity. Therefore, other detachment condition should be find to analyze this region. In
addition, the accuracy of analysis became worse in the downstream because we could not analyze
the edge of downstream due to the less of number of particles. Hence, the increase of the total
number of super particles is required to solve this.

47

References

References[1] A.S. Yen et al. An Integrated View of the Chemistry and Mineralogy of Martian Soils [J].
Nature 2005, 436, 49-54.
[2] Hayabusa Project [EB/OL]. [2010-6-8]
http://www.isas.jaxa.jp/en/missions/spacecraft/past/hayabusa.html[3] S. Abe et al. Mass and local topography measurements of Itokawa by Hayabusa [J]. Science,
2006, Vol. 312, no. 5778, pp. 1344-1347.
[4] NASA Science MARS EXPORATION PROGRAM [EB/OL]. [2018-4-6]
https://mars.nasa.gov/#red_planet/0[5] JUICE DEFINITON STUDY REPORT [RED BOOK] [EB/OL]. [2014-9-30]
http://sci.esa.int/juice/54994-juice-definition-study-report/#[6] 栗木恭一, 荒川義博. 電気推進ロケット入門 [M]. 日本:東京大学出版会, 2003: p. 20-21[7] H. Nabuchi et al., Thrust Enhanced by a Magnetic Laval Nozzle in an Applied-Field MagnetoPlasma-DynamicThruster [J]. Plasma and Fusion Research: Regular Articles 2016, 11,2406033[8] M. Colleti. Simple thrust formula for an MPD thruster with applied-magnetic field frommagnetic stress tensor [C]. in 43rd AIAA/ASME/SAE/ASEE Joint Propulsion Conference &Exhibit. 2007. Cincinnati, Ohio: AIAA.
[9] F. N. Gesto, et al., Ion Detachment in the Helicon Double-Layer Thruster Exhaust Beam [J].
Journal of Propulsion and Power, 2006. 22(1): p. 24-30.
[10] Hyde R, Jr L W and Nuckolls J [C]1972 AIAA Paper 72-1063[11] Hideki NAKASHIMA, et al., Parametric design study of laser fusion rocket [C] Proceedingsof 18th International Symposium on Space Technology and Science,1992, Vols. 1&2, 19-24[12] Yoshihiko Nagamine, and Hideki Nakashima, Analysis of Plasma Behavior in a Magnetic
Thrust Chamber of a Laser Fusion Rocket [J]. Fusion Science and Technology, 1999, 35, 62-70[13] Konstantin V. Vchivkov, et al., Optimization of thrust efficiency in laser fusion rocket by using
three-dimensional hybrid particle-in-cell code [J], Vacuum, 2004, 73, 427-432[14] Akihiro Maeno, et al., Effect of Magnetic Field Strength on a Magnetic Thrust Chamber System[J]. Journal of Propulsion and Power, 2004, 30, 54-61[15] Christopher S. Olsen, Experimental Characterization of Plasma Detachment from Magnetic
Nozzles [D]. Houston, Texas: Rice University, 2013[16] R. Kawabuchi, et al., Numerical simulation of plasma detachment from a magnetic nozzle byusing fully Particle-In-Cell code[J]. Journal of Physics: Conference Series, 2008. 112(4): p. 4[17] M. D. Carter, et al. Radio Frequency Plasma Applications for Space Propulsion [C].
International Conference on Electromagnetics in Advanced Applications. 1999. Torino, Italy.
[18] Christopher S. Olsen, et al. Investigation of Plasma Detachment From a Magnetic Nozzle in
the Plume of the VX-200 Magnetoplasma Thruster [J]. IEEE Transactions on Plasma Science ,2015 ,Volume: 43 , Issue: 1, 252-268[19] K. Terasaka, et al., Observation of Ion Streamline Detachment and Onset of Azimuthal Rotation
in a Diverging Magnetic Field [J]. IEEE Transactions on Plasma Science, 2011. 39(11): p. 24702471.
[20] K. Terasaka, et al., Experimental studies on ion acceleration and stream line detachment in adiverging magnetic field [J]. Physics of Plasmas, 2010. 17(7): p. 072106-6.
[21] A. I. Podgorny and I. M. Podgorny, MHD Simulation of Solar Flare Current Sheet Position andComparison with X-ray Observations in active region NOAA 10365 [J], Sun and Geosphere, 2013;
ISSN 1819-0839, 8(2):71-76[22] J. Buechner, Vlasov-code simulation [M]. Advanced Methods for Space Simulation, edited by
H. Usui and Y. Omura, Tokyo, TERRAPUB, 2007, pp.23-46.
[23] Y. Omura, One-dimensional Electromagnetic Particle Code: KEMPO1 [M]. Advanced48

References

Methods for Space Simulation, edited by H. Usui and Y. Omura, Tokyo, TERRAPUB, 2007, pp.121.
[24]
T. D. Arber, et al, Contemporary particle-in-cell approach to laser-plasma modeling [J],Plasma Physics and Controlled Fusion 2015, 57, 113001[25] K.S. Yee, Numerical solution of initial boundary value problems involving Maxwell’sequations in isotropic media [J]. IEEE Transactions on Antennas and Propagation. 1966, Vol.14,
No.3, 302-307[26] Hartmut Ruhl, Classical Particle Simulations with the PSC code [EB/OL]. https://www.plasmasimulation-code.net/documentation/psc-docu.pdf[27] C. K. Birdsall and A. B. Langdon, Plasma Physics via Computer Simulation [J]. Institute ofPhysics Publishing Bristol and Philadelphia, 1995, ISBN: 0-7503-0117-1, 356.
[28] Villasenor and Buneman, Rigorous charge conservation for local electromagnetic field solvers[J]. Computer Physics Communications 1992, 69, 306-316[29] User Manual for the EPOCH PIC codes [EB/OL]. [2018-4-13] University of Warwick[30] G. M. Amdahl, Validity of the Single Processor Approach to Achieving Large-Scale Computing
Capabilities [C]. AFIPS Conference, 1967, Proceedings (30): 483–485.
[31] K. V. VCHIVKOV, et al., Laser-Produced Plasma Experiments and Particle in Cell Simulation
to Study Thrust Conversion Process in a Laser Fusion Rocket [J]. Japanese Journal of Applied
Physics, 2003, 42, 6590
[32] NRL Plasma Formulary 2018 [EB/OL]. https://www.nrl.navy.mil/ppd/content/nrl-plasmaformulary , Naval Research Laboratory49

Acknowledgements

Acknowledgements
I would like to express my deepest appreciation to my supervisor Prof. Naoji Yamamoto, Assist. Prof. Taichi Morita,and Hideki Nakashima. They gave me constructive comments and warm encouragement to carry out my research andreach the completion of my Master’s thesis.
I am very thankful to Mrs. Eiko Ouchida who takes care of my paperwork. I would like to thank also the LFR
research group members; Naoya Saito, Masafumi Edamoto, Mariko Takagi, Yuki Nishioka, Toru Yamamura, and Satoru
Kawabe. In addition, I thank my labmates; Koichi Ushio, Yusuke Egawa, Ali Yousefian, Naruya Hiroike, Ryo Ikeda,
Masatoshi, Chono, Ippei Takesue, Masakazu Otsubo, Yuki Kuroe, Misaki Adachi, Jun Saiki, Nao Nishimuta, Atsumu
Inoue, Masato Kaswazu, and Naoya Kuwabara.
The computation was mainly carried out using the computer facilities at Research Institute for Information
Technology, Kyushu University.
Finally, I would like to express the deepest gratitude to my family for their kind support to help.

50