九州大学学術情報リポジトリ
Kyushu University Institutional Repository

Study on Experimental Demonstration andNumerical Simulation of Magnetic Thrust Chamber
System前野, 旭弘https://doi.org/10.15017/1485068
出版情報：九州大学, 2014, 博士（工学）, 課程博士バージョン：published
権利関係：全文ファイル公表済

Doctoral Dissertation

Study on Experimental Demonstration and Numerical Simulation ofMagnetic Thrust Chamber System
磁気スラストチャンバーシステムの実験的実証と数値シミュレーションに関
する研究

指導教員

中島

秀紀

教授

指導教員

山本

直嗣

准教授

Kyushu University
Interdisciplinary Graduate School of Engineering Science Department ofAdvanced Energy Engineering School of Engineering Sciences

Akihiro Maeno

August 6, 2014

Contents
1 Introduction

1

1.1

Interplanetary spaceflight . . . . . . . . . . . . . . . . . . . . .

1

1.2

Nuclear fusion . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.3

Inertial Confinement Fusion . . . . . . . . . . . . . . . . . . .

4

1.4

The Sun’s energy . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.5

Laser Fusion Rocket . . . . . . . . . . . . . . . . . . . . . . .

7

1.6

Magnetic thrust chamber system of LFR . . . . . . . . . . . .

7

1.7

Past works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

1.8

Aim of this study . . . . . . . . . . . . . . . . . . . . . . . . . 11

1.9

Table of contents . . . . . . . . . . . . . . . . . . . . . . . . . 12

Bibliography. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2 Laser-target interaction

18

2.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.2

Laser parameters . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.3

Elements of plasma physics
2.3.1. . . . . . . . . . . . . . . . . . . 19

Dielectric functions for transverse and longitudinal plasma
waves . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.3.2
2.4

Transverse electromagnetic waves . . . . . . . . . . . . 22

Collisional absorption of laser light in laser-produced plasma . 24
2.4.1

Scaling with wavelength and incident intensity . . . . . 24

i

CONTENTS

ii

2.4.2

Absorption coeﬃcient . . . . . . . . . . . . . . . . . . . 24

2.4.3

Model for collisional absorption . . . . . . . . . . . . . 25

2.4.4

Dependence of absorption on wavelength . . . . . . . . 25

2.5

2.6

1D Lagrangian radiation hydrodynamic code . . . . . . . . . . 26
2.5.1

Hydrodynamic equations . . . . . . . . . . . . . . . . . 26

2.5.2

Radiative transport equations . . . . . . . . . . . . . . 29

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

Bibliography. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3 3D Hybrid PIC Code

34

3.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

3.2

Particles in electromagnetic field . . . . . . . . . . . . . . . . . 35
3.2.1

PIC algorithm . . . . . . . . . . . . . . . . . . . . . . . 35

3.2.2

Digital filter . . . . . . . . . . . . . . . . . . . . . . . . 37

3.3

Basic equations . . . . . . . . . . . . . . . . . . . . . . . . . . 39

3.4

Leapfrog algorithm . . . . . . . . . . . . . . . . . . . . . . . . 41

3.5

3.4.1

Initial steps . . . . . . . . . . . . . . . . . . . . . . . . 42

3.4.2

Cycle steps . . . . . . . . . . . . . . . . . . . . . . . . 44

Calculation of electromagnetic field . . . . . . . . . . . . . . . 44
3.5.1

Plasma-vacuum interface . . . . . . . . . . . . . . . . . 44

3.5.2

Predictor corrector algorithm for plasma region . . . . 46
Initial step . . . . . . . . . . . . . . . . . . . . . . . . . 46
Cycle steps . . . . . . . . . . . . . . . . . . . . . . . . 47

3.5.3
3.6

SOR algorithm for vacuum region . . . . . . . . . . . . 48

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

Bibliography. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

4 Direct measurement of impulse

53

4.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

4.2

Dependence of impulse on laser pulse duration . . . . . . . . . 54

CONTENTS

iii

4.2.1

Experimental setup for EUV facility . . . . . . . . . . 54

4.2.2

Pendulum type thrust stand for EUV facility . . . . . . 57

4.2.3

Results and discussion for EUV facility . . . . . . . . . 59

4.3

Dependence of impulse on laser pulse wavelength and incident
laser energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

4.4

4.3.1

Experimental setup for GXII facility . . . . . . . . . . 62

4.3.2

Pendulum type thrust stand for GXII facility . . . . . 62

4.3.3

Results and discussion for GXII facility . . . . . . . . . 66

Comparison of measured and calculated impulse . . . . . . . . 68
4.4.1

Numerical simulation with and without considering creation process of laser-produced plasma . . . . . . . . . 68

4.4.2
4.5

Results and discussion for numerical simulations . . . . 72

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

Bibliography. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

5 Comparison of ion expansion velocity

78

5.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

5.2

Experiment for measuring ion expansion velocity . . . . . . . . 79
5.2.1

Experimental setup . . . . . . . . . . . . . . . . . . . . 79

5.2.2

Charge collector with a Faraday cup . . . . . . . . . . 81

5.2.3

Pulsed power supply for generating high pulsed current 83
Pulse Forming Network . . . . . . . . . . . . . . . . . . 83
Rogowski coil . . . . . . . . . . . . . . . . . . . . . . . 86

5.3

Numerical simulation with and without considering creation
process of laser-produced plasma . . . . . . . . . . . . . . . . 90

5.4

5.3.1

Calculation model of Star-1D . . . . . . . . . . . . . . 90

5.3.2

Results and discussion of Star-1D . . . . . . . . . . . . 92

5.3.3

Calculation model of 3D hybrid PIC code . . . . . . . 94

Results and discussion . . . . . . . . . . . . . . . . . . . . . . 96

CONTENTS
5.5

iv

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

Bibliography. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

6 Eﬀect of magnetic field

100

6.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

6.2

Dependence on ion expansion velocity distribution in initiallaser-produced plasma . . . . . . . . . . . . . . . . . . . . . . 101
6.2.1

Calculation model of 3D hybrid PIC code . . . . . . . 101

6.2.2

Results and discussion . . . . . . . . . . . . . . . . . . 104
Diamagnetic cavity . . . . . . . . . . . . . . . . . . . . 104
Thrust energy and impulse . . . . . . . . . . . . . . . . 108
Discussion . . . . . . . . . . . . . . . . . . . . . . . . . 110

6.3

6.4

Thrust performances for various magnetic field strengths . . . 112
6.3.1

Numerical simulation conditions . . . . . . . . . . . . . 112

6.3.2

Results and discussion . . . . . . . . . . . . . . . . . . 112

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

7 Summary118

A Ampere’s law with Darwin approximation

123

B Components of matrix equations

127

n+ 12

B.1 x component of vi. . . . . . . . . . . . . . . . . . . . . . . 127

B.2 x component of E n . . . . . . . . . . . . . . . . . . . . . . . . 130
C Neumann boundary condition

131

D Normalization

135

E CFL condition

136

List of Figures
1.1

Four stages of a direct-drive implosion. . . . . . . . . . . . . .

1.2

Mechanism for generating impulse in the magnetic thrust chamber system of LFR. . . . . . . . . . . . . . . . . . . . . . . . .

5

9

2.1

Target irradiated and focused by a single laser beam. . . . . . 20

2.2

Normal incidence of the laser light on planar ablation plasma.

2.3

Ion and electron temperature, ablation plasma velocity and27electron density distributions in a laser-produced plasma. . . . 27
3.1

Particle and eight nearest grid points. • and ◦ shows particle
and grid point, respectively. . . . . . . . . . . . . . . . . . . . 36

3.2

Eighteen nearest grid points. The points indicate middle, side
and corner points. . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.3

Update steps in time of leapfrog algorithm. . . . . . . . . . . . 43

3.4

A schematic of the grid points used in the plasma and vacuum
region. The black area is the outside of calculation region. . . 45

4.1

Vacuum chamber in the EUV facility. . . . . . . . . . . . . . . 55

4.2

Experimental setup for EUV facility. . . . . . . . . . . . . . . 55

4.3

Calculated magnetic field strength produced by the permanent
magnets at the z-axis at x = y = 0 and in the z − x plane at
y = 0. The magnetic field strength at the origin is 0.12 T. . . 56v

LIST OF FIGURES
4.4vi

Schematic of the pendulum type thrust stand for EUV facility.
(vertical sectional view) . . . . . . . . . . . . . . . . . . . . . 58

4.5

Applied calibration force as a function of the signals of LED
sensors and the liner curve fitting to the calibration data. . . . 58

4.6

Impulses measured on the pendulum type thrust stand as afunction of incident laser energies for laser pulse duration of 3
ns and 9 ns. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

4.7

Impulse model of the product of plasma ablation mass and ion
acoustic speed. . . . . . . . . . . . . . . . . . . . . . . . . . . 61

4.8

Vacuum chamber in the GXII facility. . . . . . . . . . . . . . . 64

4.9

Experimental setup & Schematic of the pendulum type thruststand for GXII facility. . . . . . . . . . . . . . . . . . . . . . . 64

4.10 Calculated magnetic field strength produced by the permanent
magnet at the z-axis at x = y = 0 and in the z − x plane at
y = 0. The magnetic field strength at the origin is 0.1 T. . . . 65
4.11 Applied calibration force as a function of the signals of LED
sensors and the liner curve fitting to the calibration data. . . . 65
4.12 Impulses measured on the pendulum type thrust stand as afunction of incident laser energies for laser pulse wavelength. . 67
4.13 Maxwell super particle velocity distribution. The characteristic velocity is assumed to be vp =2.4 km/s. . . . . . . . . . . . 70
4.14 Cubic calculation model of 3D hybrid PIC code. . . . . . . . . 71
4.15 Unit vector distributions of diamagnetic current calculated bymodified and original numerical simulation. . . . . . . . . . . . 73
5.1

Experimental setup & Cubic calculation model of 3D hybrid
PIC code. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

5.2

Edgewise coil for generating a pulsed magnetic filed. . . . . . . 80

LIST OF FIGURES
5.3

vii

Charge collector with a Faraday cup and an I-V converter.
The 2-mm-diameter hole is covered with a mesh to prevent
electrons in the plasma from entering the Faraday cup. . . . . 82

5.4

Circuit diagram of pulsed power supply with a 3-step PFN. . . 85

5.5

Pulsed power supply made by myself. . . . . . . . . . . . . . . 85

5.6

Rogowski coil for measuring a pulsed current. . . . . . . . . . 87

5.7

Equivalent circuit diagram of integral and non-inverting amplifier circuits. . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

5.8

Experimental and theoretical output currents of pulsed power
supply at the peak current of 450 A. Current starts flowing 31µs before laser ignition. . . . . . . . . . . . . . . . . . . . . . . 89

5.9

Lagrangian calculation model for ray-tracing. . . . . . . . . . . 91

5.10 Gaussian time-variable power profile of a single laser beam. . . 91
5.11 1D distributions of plasma for ion density, ion expansion velocity, and charge state as a function of plasma radius at 75
ns. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
5.12 Time-variation of plasma kinetic, internal and radiation energies. 93
5.13 Calculated magnetic field strength produced by the edgewise
coil at the peak current of 450 A. The magnetic field strength
at the origin is 0.0321 T. . . . . . . . . . . . . . . . . . . . . . 95
6.1

Cubic calculation model of 3D hybrid PIC code. . . . . . . . . 102

6.2

Calculated magnetic field strength produced by the edgewise
coil when the peak current is 5.2 kA. The magnetic field strength
at the center is 0.51 T. . . . . . . . . . . . . . . . . . . . . . . 102

6.3

Magnetic field strength distribution for the original case inz − x plane at y = 0 for calculation times of 325, 575, 825 ns,
and 1.075 µs. A diamagnetic cavity is formed. . . . . . . . . . 105

LIST OF FIGURES
6.4

viii

Magnetic field strength distribution for the original case inz − x plane at y = 0 for calculation times of 1.325, 1.575,
3.075, and 3.575 µs. A diamagnetic cavity is formed. . . . . . 106

6.5

Magnetic field strength distribution for the modified case inz −x plane at y = 0 for calculation times of 1.325, 2.575, 3.825,
and 5.075 µs. A diamagnetic cavity is formed. . . . . . . . . . 107

6.6Time evolutions of thrust energy and impulse. For the originalcase, thrust energy and impulse saturate at 0.36 J (31%) and151 µNs (41%), respectively. For the modified case, thrustenergy and impulse saturate at 0.37 J (12%) and 31.6 µNs(41%), respectively. . . . . . . . . . . . . . . . . . . . . . . . . 109

6.7Time evolutions of thrust energy and impulse for (a) slow velocity group and (b) fast velocity group. For group A, thrustenergy and impulse saturate at 0.0337 J (125%) and 24.1 µNs(102%). For group B, thrust energy and impulse saturate at
0.34 J (11%) and 7.5 µNs (14%), respectively. . . . . . . . . . 111

6.8

Ion super particle distributions at B=0.0321 T calculated bymodified and original numerical simulations. . . . . . . . . . . 114

6.9

Ion super particle distributions at B=0.961 T calculated bymodified and original numerical simulations. . . . . . . . . . . 114

6.10 Ion super particle distributions at B=1.92 T calculated bymodified and original numerical simulations. . . . . . . . . . . 115
6.11 Ion super particle distributions at B=2.88 T calculated bymodified and original numerical simulations. . . . . . . . . . . 115
6.12 Impulse for various magnetic field strengths calculated by modified and original numerical simulations. . . . . . . . . . . . . 116
C.1 Orthogonal and nonorthogonal cells. . . . . . . . . . . . . . . 132
C.2 Two examples of Neumann boundary condition. . . . . . . . . 132

List of Tables
4.1

Plasma conditions for original and modified cases. . . . . . . . 71

5.1

Plasma conditions for original and modified cases. . . . . . . . 95

6.1

Plasma conditions for original and modified cases. . . . . . . . 103

ix

Chapter 1
Introduction
1.1

Interplanetary spaceflight

Great progress has already been made in exploring the earth and themoon with manned and unmanned missions and surveying other planetsin the solar system by using unmanned probes. In the last 18 years, Jet
Propulsion Laboratory (JPL) of National Aeronautics Space Administration(NASA) has strongly carried out some robotic space missions to Mars as partof JPL’s Mars Exploration Program [1].
Firstly, “Mars Pathfinder” [2] was launched out in December 1996 fromthe earth and landed on Mars in July 1997. This spaceflight spent 212 days
going to the planet. In this mission, from the investigations by scientificinstruments on both the lander and the rover, it was suggested that Mars
had been warm and existed liquid water in the atmosphere. Eight years
later, in January 2004, two robotic geologists (Spirit and Opportunity) [3]landed on Mars after 202 days from they were launched out. This mission
is named “Mars Exploration Rover.” These robotic explorers have trekked
across the Martian surface with surveying the field geology and observing
the atmosphere. As a result, they have found that there were wet and habit-

1CHAPTER 1. INTRODUCTION

2

able conditions on ancient Mars. NASA’s next rover mission, “Mars Science
Laboratory” [4], carried “Curiosity” to Mars. The rover arrived on Marsin August 2012 after a long trip of 254 days. Mars Science Laboratory has
been collecting Martian soil and rock samples. These samples have been
analyzed whether organic compounds and environmental conditions which
could have supported microbial life are existing now or in the past. This
rover also carries a radioisotope power system generating electricity from theheat of plutonium’s radioactive decay. This power source enables Curiosity
to operate a full Martian year (687 days) on Mars’ surface.
The aim of JPL’s Mars Exploration Program is to investigate the presenceof liquid water on Mars, either existed in its ancient past or preserved in thesubsurface today. Water is essential for lives. If Mars once had liquid wateror still does today, it has the potential to support life. From these factors,
the focus for manned transport in a century will undoubtedly be Mars.
However, it takes long time to Mars using electric propulsion (i.e., lowthrust and high specific impulse). (Specific impulse is the ratio of the effective exhaust velocity and the acceleration due to gravity.) This meansthat spaceship crew stay in the harmful radiation environment of the interplanetary space for a long time. It is too dangerous for human. Thus, space
propulsion systems which is capable of generating large thrust with high specific impulse will be required for manned interplanetary flights. Large thrust
permits high accelerations up to a speed. In the systems, fusion-based nuclear propulsion has the potential to enable fast interplanetary transaction
by the huge energy generated by fusion reaction.

1.2

Nuclear fusion

Einstein’s Special relativity, probably the most famous equation in physics,
is commonly written as E = mc2 . This equation represents the amount ofCHAPTER 1. INTRODUCTION

3energy associated with the forces holding nucleus together. When light elements combine into a heavier element, the mass of the heavier element isslightly less than the sum of the mass of the light elements. The diﬀerencein mass per nucleus is called mass defect. The mass defect multiplied by the
speed of light squared is the release energy by nuclear fusion [5]. Nuclear
fusion is about 40 billion times more eﬃcient in converting mass into energy than chemical reactions, and about 4 times more eﬃcient than nuclear
fission.
One promising nuclear fusion reaction is between deuterium and tritium.
Deuterium-tritium (DT) reaction releases 17.6 MeV and produces helium-4(4 He) and neutron (n). Deuterium is naturally present in water and available
soon. Tritium is not available naturally and has to be produced from lithium
being available widely.

D + T → 4 He (3.5MeV) + n (14.1MeV)(1.1)

Another reaction is deuterium-helium-3 (D3 He) reaction. This reaction
releases 18.4 MeV and produces only charged particle as helium-4 and afree proton (p), not neutron. However, there is a major problem to obtainhelium-3 because it is not present on the earth. Although the moon’s surface
contains helium-3 brought by solar wind, the concentration is extremely low
and it would require enormous mining eﬀort to extract helium-3.

D + 3 He → 4 He (3.7MeV) + p (14.7MeV)(1.2)

To initiate these reactions, a lot of energy is needed to bring two nuclei
into an encounter close enough. This is because the electric force trying
to keep them apart dominates stronger than the attractive nuclear force
between two nuclei. One typical approach is magnetic confinement fusion

CHAPTER 1. INTRODUCTION

4(MCF). Microwaves or electricity heat and ionize a fuel gas in a vacuum
chamber into a plasma state. Magnetic field guides the charged particles
along the direction of the field and heats the fuel up. Another is inertial
confinement fusion (ICF). Laser beams compress and heat so quickly before
the fuel escapes.

1.3

Inertial Confinement Fusion

The ICF [6, 7] releases large amounts of energy and creates high energy
density plasma in about 10 ps by focusing homogeneously laser beams onto amillimeter-sized capsule filled with nuclear fusion fuel. When the density ofDT fuel is 0.3 g/cm3 , 425 MJ is released from a 2-mm-diameter capsule filledwith DT fuel at 100% burn up. This energy equals to the energy releasedby the combustion of about 10 kg of gasoline. The main sequences of adirect-drive implosion are as follows and sketched in Fig. 1.1.
(a) Strong laser beams are homogeneously focused onto a millimeter-sized
spherical fuel target.
(b) The laser energy is absorbed on the critical surface and generates ablation plasma on the outer layer. High pressure of the ablation plasma
implodes inwards via a rocket eﬀect.
(c) The capsule is compressed to a density several thousand times greaterthan that of solid materials. In the center of capsule, the kinetic energy of the compression is converted into the internal energy. When
suﬃcient density and temperature is reached, a central hot spot regionis created and ignites via fusion reactions.
(d) Then, a burn wave starting from the hot spot propagates and ignites
the whole fuel. The nuclear fusion energy and productions are released.

CHAPTER 1. INTRODUCTION

5

Ablation plasma

Laser

Thermal energy(a) Laser heating of the capsule outer(b) Ablation plasma of the outer layer

layer.

compresses the capsule.

Burning(c) Core reaches the density and tem-(d) Fusion reaction spreads rapidly

perature enough for ignition.

through the compressed fuel.

Figure 1.1: Four stages of a direct-drive implosion.

CHAPTER 1. INTRODUCTION

1.4

6

The Sun’s energy

ICF is the same principle as the Sun [8]. The Sun has a plenty of hydrogenfor nuclear fusion. The nuclear fusion reaction in the Sun is the conversionof four protons into a helium-4. This reaction does not happen in a single
step. Firstly, two protons combine to form a deuterium as shown in (1.3).
The deuterium combines with another proton to form helium-3 as shown in(1.4). Finally, helium-4 and two protons are released from two helium-3 asshown in (1.5). The total amount of energy released from this reaction isabout 10 million times more than the chemical reaction.

3

p + p → D + positron(1.3)

D + p → 3 He

(1.4)

He + 3 He → 4 He + 2p

(1.5)

The Sun has high dense core and low dense outer layers. The core exists
from the center to 0.3 of the Sun’s radius (about 0.7 million km). The hydrogen fuel is compressed to very high density (about 150 g/cm3 ), roughly 150
times denser than water. The temperature in the core is about 15 million
K. Nuclear fusion reactions take place only in the core. The combinationof high density and high temperature makes an enormous pressure for expanding outward. To prevent this pressure, the Sun confines the expansion
fusion plasma by the presence of a massive gravitational field of each layer.
Gravity can confine the hot plasma despite a very weak force compared to
the forces of nuclear physics. This is called hydrostatic equilibrium. In theother hand, the energy is transported outward from the core as the radiation
heat through the radiative zone to 0.7 of the radius. As the distance increases from the core, the zone becomes more opaque and radiation becomesCHAPTER 1. INTRODUCTION

7

less eﬃcient. The temperature in this zone cools to 2 million K. Then, theenergy begins to move through the convective zone outside 0.3 of the Sun’sradius. Finally, the energy arrives at the photosphere (the temperature ofabout 6,000 K) emitting the sunlight. A series of these events is quite a slow
process and it takes about a million years for getting out the energy.

1.5

Laser Fusion Rocket

The first concept of fusion-based nuclear propulsion [9] was proposed byUlam of Los Alamos National Laboratory (LANL) in 1946 and extended
further by Ulam and Everett in 1955. Then, Taylor and Dyson [10] conducted Project Orion at the LANL, in the 1950s. Since this project wasplanning to use hydrogen bomb in space, it was prohibited by the Partial
Test Ban Treaty in 1963. Then, Project Daedalus [11] originated by Winterberg [12] was motivated by the British Interplanetary Society. At the sametime, the application of ICF to space propulsion was firstly proposed in the
1970s by Balcomb [13] of Los Alamos Scientific Laboratory and Hyde, Woodand Nuckolls [14] of the Lawrence Livermore National Laboratory (LLNL).
The propulsion system named as laser fusion rocket (LFR) [15] redirects theemitted charged particles in one direction by a magnetic field. Orth [16, 17]of LLNL estimated that an LFR employing DT reaction have an eﬀective
specific impulse of 15,500 s and a thrust of 230 kN. The peak magnetic field
is 12 T. The period round-trip to Mars was also estimated to be close to 145
days. This is much shorter than the existing propulsion systems.

1.6

Magnetic thrust chamber system of LFRLFR employs a magnetic thrust chamber system to generate impulse.
The system produces impulse by inducing a diamagnetic current in a laser-CHAPTER 1. INTRODUCTION

8

produced plasma as the following processes. Figure 1.2 depicts the mechanism for producing impulse in the system of LFR.
(a) ICF plasma expands in an applied magnetic field with a suitably designed geometry. The magnetic field is typically produced by a solenoid
super-conducting magnet (SCM).
(b) Since plasmas are good conductors, the charged particles in the plasma
undergo Larmor motion when a magnetic field is applied. This Larmor
motion induces the diamagnetic current to sweep aside and compress
the magnetic field.
(c) The compressed magnetic field pushes the plasma back in a manner
similar to a spring. Then, LFR generates impulse through the interaction between the diamagnetic current in the plasma and the magnetic
field generated by the SCM.

CHAPTER 1. INTRODUCTION

Coil

9

Magnetic field

Expanding plasma(a)(a) Expanding ICF plasma in an appliedmagnetic field that has a suitably designed
geometry.

Compressed magnetic field

Diamagnetic current(b)(b) Compressed magnetic field and diamagnetic current induced by Larmor motion.

Impulse(c)(c) Compressed magnetic field pushes
plasma back.

LFR produces impulse

through the interaction between the diamagnetic current and the magnetic field.

Figure 1.2: Mechanism for generating impulse in the magnetic thrust chamber system of LFR.

CHAPTER 1. INTRODUCTION

10

The magnetic thrust chamber system of LFR has the following threeadvantages.
1. The galactic cosmic radiation is the greatest problem associated withmanned spaceflight. LFR reduces the period of space missions and theexposure of the spaceflight crew to cosmic radiation because LFR cangenerate large thrusts with high specific impulses.
2. The thrust can be easily varied during the mission by varying the system parameters. For example, the thrust can be changed by varying theoperation frequency of the system because the thrust is proportional
to the production frequency of impulses.
3. This system has a long lifetime. The direct plasma impingement on
the solid wall of the vehicle is minimized as the plasma is confined by
the magnetic field produced by the SCM. Consequently, little charged
plasma energy is wasted and the damage of the spacecraft structure isminimized, resulting in long lifetime.

1.7

Past works

The most important thing for realizing LFR is the magnetic thrust chamber system. Several numerical simulations and experiments have been performed to investigate the behavior of charged particles in a magnetic field
for the magnetic thrust chamber system of LFR. In experiments, Vchivokv
et al. [18] observed the temporal plasma behavior in an axially symmetric
dipole magnetic field. By using a magnetic probe, they discovered that a diamagnetic cavity forms in the applied magnetic field. Pioneer works related
to the plasma cloud expansion were conducted by a Japanese group. Sato etal. [19] studied experimentally and numerically the physical mechanism ofreconverting the kinetic expansion energy of plasma into thermal energy asCHAPTER 1. INTRODUCTION

11

the result of plasma-magnetic field interaction and the occurrence of plasma
instability. In numerical simulations [20, 21, 22], Rayleigh-Taylor instability
and electrons detachment from magnetic field were suﬀered. Rayleigh-Taylor
instability is a serious problem for plasma explosion and would lead to theimpulse reduction. The detachment problem is caused by the Larmor radiusdiﬀerence between ions and electrons. When electrons would not escape fromthe magnetic field of the system, ions would be pulled back to the space vehicle by the charge-separation electric field generated by the Larmor motionsof the ions and electrons. As a result, the space vehicle cannot get enoughimpulse. Nagamine et al. [20] found that Rayleigh-Taylor instability would
not be serious problems for the system by using a three-dimensional (3D)
hybrid Particle-In-Cell (PIC) code which simulates the behavior of charged
particles in a magnetic field. Kawabuchi et al. [21] found that electrons detachment would occur and generate impulse by using a two-dimensional (2D)
coordinate space and a 3D velocity space electromagnetic full PIC code.

1.8

Aim of this study

Past works have evaluated the behavior of charged particles in a magnetic field for the magnetic thrust chamber system of LFR by using mainly
numerical simulation. However, there have already left same challenges inthe study of the magnetic thrust chamber system for LFR. To study threechallenges of them, the followings are conducted.
1. The direct measurement of impulse due to the interaction between adiamagnetic current and a magnetic field has not been demonstrated
experimentally yet. To produce the impulse by using laser-produced
plasma and magnet, we attempt to conduct experiments in the actual
research laboratory instead of ICF plasma and SCM. The experiments
are conducted for various laser pulse energies at Extreme Ultra-Violet

CHAPTER 1. INTRODUCTION

12(EUV) database and Gekko XII (GXII) laser facilities of Institute ofLaser Engineering (ILE) in the collaboration with Osaka University.
2. The creation process of a laser-produced plasma has not suﬃciently
been considered in the numerical simulation yet. This process is needed
to compare numerical simulations with experiments in the actual research laboratory. By combining a 3D hybrid PIC code with a onedimensional (1D) Lagrangian radiation hydrodynamic code which cancompute the hydrodynamic evolution of radiating plasma heated bylaser pulses, we attempt to compute high energy laser injection onto aspherical target and the ablation plasma dynamics in a magnetic field.
3. To find the scaling law of the magnetic thrust chamber system, it isnecessary to confirm the validity of the numerical simulation. In order to validate the modified numerical simulation with considering thecreation process, the experimental results are compared with these of
the numerical simulations which are with and without considering thecreation process of a laser-produced plasma.

1.9

Table of contents

The structure of this thesis is as follows. Chapter 2 reviews the energy
deposition of high power beams and the 1D Lagrangian radiation hydrodynamic code via the deposition. Chapter 3 describes the basic equations andthe numerical algorithms used in the 3D hybrid PIC code. Chapter 4 is devoted to the direct measurement of impulse and the dependences of impulseon the laser parameters. Numerical simulations with and without considering the creation process of a laser-produced plasma are also performed toconfirm the validity of the modified numerical simulation. In Chapter 5, theion expansion velocity of a laser-produced plasma is described. NumericalCHAPTER 1. INTRODUCTION

13

simulations are also performed to confirm the validity of the modified numerical simulation. In addition to these, the interaction between the initialplasma condition and magnetic field strength is investigated in Chapter 6.
Finally, Chapter 7 contains conclusions and suggestions for the future work.

Bibliography[1] JPL, “Mars Exploration Program,” URL: http://mars.jpl.nasa.gov[cited 1 January 2013].
[2] JPL, “NASA Facts:

Mars Pathfinder,” NASA Fact Sheet, URL:

http://www.jpl.nasa.gov/news/fact_sheets/mpf.pdf [cited 1 January 2013].
[3] JPL, “NASA Facts: Mars Exploration Rover,” NASA Fact Sheet, URL:
http://marsrovers.jpl.nasa.gov/newsroom/factsheets/pdfs/
Mars03Rover041020.pdf [cited 1 January 2013].
[4] JPL, “NASA Facts: Mars Science Laboratory/Curiosity,” NASA Fact
Sheet, URL: http://mars.jpl.nasa.gov/msl/news/pdfs/MSL_Fact_
Sheet.pdf [cited 1 January 2013].
[5] McCracken, G. M., and Stott, P. E., Fusion: The Energy of the Universe,
1st ed., Academic Press, Massachusetts, 2005, Chaps. 1 and 2.
[6] Atzeni, S., and Meyer-ter-Vehn, J., The Physics of Inertial Fusion:
Beam Plasma Interaction, Hydrodynamics, Hot Dense Matter, Oxford
University Press, London, 2004, Chap. 3.
[7] McCracken, G. M., and Stott, P. E., Fusion: The Energy of the Universe,
1st ed., Academic Press, Massachusetts, 2005, Chap. 7.

14

BIBLIOGRAPHY

15[8] McCracken, G. M., and Stott, P. E., Fusion: The Energy of the Universe,
1st ed., Academic Press, Massachusetts, 2005, Chap. 3.
[9] Long, K. F., Deep Space Propulsion: A Roadmap to Interstellar Flight,
Springer, Berlin, 2012, Chap. 12.
[10] Dyson, F. J., “Interstellar Transport,” Physics Today, October 1968, pp.
41-45. doi: 10.1063/1.3034534[11] Project Daedalus Study Group and British Interplanetary Society,
Project Daedalus: The Final Report on the BIS Starship Study, Space
Educational Aids, London, 1978.
[12] Winterberg, F., “Rocket Propulsion by Thermonuclear Microbombs Ignited with Intense Relativistic Electron Beams,” Raumfahrtforschung,Vol. 15, 1971, pp. 208-217.
[13] Balcomb, J. D., Booth, L. A., Cotter, T. P., Hedstrom, J. C., Robinson,
C. P., Springer, T. E., Watson, C. W., “Nuclear Pulsed Space Propulsion
Systems,” Los Alamos Scientific Laboratory LA-4541-MS, 1970.
[14] Hyde, R. A., Wood, L. L., and Nuckolls, J. H., “Prospects for Rocket
Propulsion with Laser-Induced Fusion Microexplosions,” AIAA Paper
72-1063, 1972.
[15] Hyde, R. A., “A Laser Fusion Rocket for Interplanetary Propulsion,”
Lawrence Livermore National Laboratory UCRL-88857, 1983.
[16] Orth, C. D., Klein, G., Sercel, J., Hoﬀman, N., Murray, K., and ChangDiaz, F., “Transport Vehicle for Manned Mars Missions Powered by
Inertial Confinement Fusion,” Lawrence Livermore National Laboratory
UCRL-96832, 1987.

BIBLIOGRAPHY

16[17] Orth, C. D., “VISTA - A Vehicle for Interplanetary Space Transport Application Powered by Inertial Confinement Fusion,” Lawrence Livermore
National Laboratory UCRL-LR-110500, 2003.
[18] Vchivokv, K. V., Nakashima, H., Zakharov, Y. P., Esaki, T., Kawano,T., and Muranaka, T., “Laser-Produced Plasma Experiments and Particle in Cell Simulation to Study Thrust Conversion Processes in a Laser
Fusion Rocket,” Jpn. J. Appl. Phys., Vol. 42, Issue 10, 2003, pp. 65906597. doi: 10.1143/JJAP.42.6590[19] Sato, K., Okada, S., Kogoshi, S., Sudo, S., Tsuji, H., Ohwadano,Y.,

Sekiguchi,

T.,

Katsurai,

M. and Tanaka,

H.,“HIGH-TEMPERATURE HIGH-QUALITY DEUTERIUM PLASMA PRODUCTION BY LASER BEAM(S) AND INTERACTIONS WITHMAGNETIC FIELDS,” Proceedings of the 6th International Conferenceon Plasma Physics and Controlled Nuclear Fusion Research, held by theInternational Atomic Energy Agency, at Berchtesgaden, 6-13 October
1976, Nuclear Fusion Supplement, Vol. 2, 1977, pp. 567-575.
[20] Nagamine, Y., and Nakashima, H., “Analysis of Plasma Behavior in a
Magnetic Thrust Chamber of a Laser Fusion Rocket,” Fusion Scienceand Technology, Vol. 35, Issue 1, 1999, pp. 62-70.
[21] Kawabuchi, R., Matsuda, N., Kajimura, Y., Nakashima, H., and Zakharov, Y. P., “Numerical Simulation of Plasma Detachment froma Magnetic Nozzle by using Fully Particle-In-Cell Code,” J. Phys.
Conf. Ser., Vol. 112, Issue 4, 2008, pp. 1-4. doi:

10.1088/1742-

6596/112/4/042082[22] Kajimura, Y., Kawabuchi, R., and Nakashima, H., “Control Techniquesof Thrust Vector for Magnetic Nozzle in Laser Fusion Rocket,” Fusion

BIBLIOGRAPHY

17

Engineering and Design, Vol. 81, Issues 23-24, 2006, pp. 2871-2875. doi:
10.1016/j.fusengdes.2006.07.046

Chapter 2
Laser-target interaction
2.1

Introduction

For creation of a laser-produced plasma, a fuel target is driven by externallaser beams. The beams heat the target rapidly to high temperatures. The
beam energy is transferred predominantly to electrons and electron plasma
waves. Then, thermalization occurs due to Coulomb collisions between theelectrons and the ions and a laser-produced plasma is created. Here, we
describe the energy deposition of high power beams and a 1D Lagrangianradiation hydrodynamic code via the deposition. This code can compute thehydrodynamic evolution of radiating plasma heated by laser pulses.

2.2

Laser parameters

The important parameter that dominates the creation process of a laserproduced plasma is the critical density nc ,nc =4π 2 meµe2 λ2L18(2.1)

CHAPTER 2. LASER-TARGET INTERACTION

λL =

2πc
ω

19

(2.2)

where me is the electron mass, µ is the magnetic permeability and e is theelementary charge. λL is the laser wavelength as follows, where c is the speedof light (i.e., c2 = 1/µε) and ω is the laser frequency. When the target isdriven by a laser beam, the beam is reflected on the critical surface andcannot penetrate inside. The incident laser intensity, IL , is defined as thelaser power per unit area per one shot as shown in Fig. 2.1, and is given by

IL =

EL
2
πrspot · τL

where EL is the incident laser energy, rspot is the radius of the laser beam
spot and τL is the laser pulse duration, respectively.

2.3

Elements of plasma physics

2.3.1

Dielectric functions for transverse and longitudinal plasma waves

The absorption of external laser beams involves collective plasma eﬀects.
The theoretical study is based on Maxwell’s equations and the linearized
plasma equations, describing small perturbations in an uniform plasma. Equations (2.3) and (2.4) are Maxwell-Faraday equation and the equation forAmpere’s law, respectively.

∂B
∂t
(
)
∂E
∇×B =µ j+ε∂t∇×E =−

(2.3)
(2.4)

CHAPTER 2. LASER-TARGET INTERACTION

Laser

rspot
TargetFigure 2.1: Target irradiated and focused by a single laser beam.

20

CHAPTER 2. LASER-TARGET INTERACTION

21

where E is the electric field, and B is the magnetic field, t is time and ε isthe permeability. By substituting Eq. (2.3) into the time partial diﬀerential
equation of Eq. (2.4), we obtain the following equation.

∇ × (∇ × E) = −µ∂j1 ∂ 2E
− 2 2
∂t
c ∂t

(2.5)

Equation (2.5) describes the propagation of electromagnetic waves in aplasma. For the sake of conveniences, all perturbations are expressed in termsof their Fourier modes as shown below,

A = ℜ[Â exp(−iωt + ik · r)]
where Â is the complex magnitude and r is the position. Time and space
derivatives are replaced by the frequency −iω and the wave vector ik, respectively. Then, Eq. (2.5) has the form by Fourier transformation (i.e.,∇ = ik, ∂/∂t = −iω), and the wave equation is−k × (k × E) = iµωj +
=

ω2
E
c2

ω2
εE
c2

(2.6)

Here, from the Ohm’s law j = σE, the definition of dielectric tensor ε is

ε=1+i

µc2
σ
ω

(2.7)

where σ is the electrical conductivity tensor. For non-magnetized plasmas,
the dielectric tensor has the only diagonal elements, εT for the two transverse
directions (E ⊥ k) and εL for the longitudinal direction (E ∥ k; k × E = 0).
Consequently, Eq. (2.6) separates into a transverse part shown in Eq. (2.8)
and a longitudinal part shown in Eq. (2.9) [1, 2].

CHAPTER 2. LASER-TARGET INTERACTION|k|2 =

ω2εT
c2εL = 0

22

(2.8)
(2.9)

In the case of k × (k × E) = 0, the plasma and displacement currents
compensate each other. The dielectric functions εT and εL have to be determined from the plasma equations and they depend on the wave number |k|and ω for an isotropic plasma.

2.3.2

Transverse electromagnetic waves

In the lowest linear order, when the ion is assumed to be static, theelectromagnetic waves induce a transverse perturbation in the plasma electrons. This is described by Newton’s equation for the electron velocity ve
and Fourier transformation.

medve= −iωme ve = −eE − me νei ve
dt(2.10)

Equation (2.10) includes electron-ion collisions, where νei is the electronion collision frequency. The friction term (∝ νei ) is essential for collisional
absorption. The electron velocity is given by Eq. (2.10).

ve = −eEme (νei − iω)(2.11)

From Eq. (2.11), the transverse current density is obtained.

j = −ne eve =ne e2
E = σEme (νei − iω)(2.12)

and from Eqs. (2.7) and (2.12), the transverse dielectric function is given by

CHAPTER 2. LASER-TARGET INTERACTION

23

2
iωpe
ε
σ=
ω + iνei

ne e2
2
ωpe =
εme
2
ωpe
εT = 1 −ω(ω + iνei )(2.13)

where ωpe is the electron plasma frequency and ne is the electron density.
The imaginary term in Eq. (2.13) leads to damping of the electromagnetic
wave. Replacing k by k + iKa /2, the incident laser intensity becomes IL ∝|E|2 ∝ exp(−Ka r), where Ka is the absorption coeﬃcient. By inserting Eq.
(2.13) into Eq. (2.8), the dispersion equation is given as Eq. (2.14). From
the real part of Eq. (2.14), the dispersion relation for the electromagnetic
waves, Eq. (2.15), is obtained. The coeﬃcient for collisional absorption, Eq.
(2.16) is obtained from the imaginary part of the dispersion equation [1, 2].

|k|2( c )2ω=1−( ω )2peωω
ω + iνei(2.14)

Here, we assumed ω ≫ νei . The group velocity vg is the speed of theoverall shape of the wave’s amplitude propagation through the space.

ω
k=
c

Ka =

√
1−

( ω )2peωνei ( ωpe )2vg ω√( ω )2
∂ω
pe
vg =
=c 1−
∂k
ω

(2.15)(2.16)(2.17)

CHAPTER 2. LASER-TARGET INTERACTION

2.4

24

Collisional absorption of laser light in laserproduced plasma

2.4.1

Scaling with wavelength and incident intensity

When ωpe is larger than ω, Eq. (2.15) is not valid. This means thatthe laser beam is reflected on the critical surface defined from the relation2
ωpe/ω 2 = ne /nc . The beam is not propagating inside, as shown in Fig. 2.2.

Here, ω 2 = nc e2 /εme is obtained by substituting Eq. (2.2) into Eq. (2.1).
For driving matters with laser beams, the collisional absorption of the light in
the ablation plasma is the dominant absorption mechanism, and it is called
inverse-bremsstrahlung [3]. Other mechanisms such as resonance and parametric absorptions are not considered because they mostly convert the lightenergy to the energy of the plasma waves and produce superthermal electrons.
In general, such mechanisms do not cause local thermalization on the criticalsurface. However, the deposited energy by the inverse-bremsstrahlung does
local thermalization. The light can travel through the ablation plasma up
to the critical surface. For suﬃciently small wavelengths, this collisional absorption is considered to be the dominant process and assumed to be about
1, since the critical density strongly depends on the laser wavelength as seenin Eq. (2.1).

2.4.2

Absorption coeﬃcient

The coeﬃcient of the collisional absorption [4, 5], given by Eq. (2.16),
2/ω 2 = ne /nc ,may be rewritten in the following form by using ωpeCHAPTER 2. LASER-TARGET INTERACTION

νei
Ka =
vg

(

ne
nc

25

)2

√
( )2
1 4 2me π 3/2 c2 e2 Z ∗ ln Λei ne=vg
3λ2L (kB Te )3/2
nc

(2.18)

where Z is the charge state of ion, ln Λei is the Coulomb logarithm [6] forelectron-ion collisions, kB is the Boltzmann constant and Te is the electron
temperature, respectively. Z ∗ is defined as < Z 2 > / < Z >, where <>indicates averaging over the ion species. The collisional frequency, νei , on thecritical surface is taken from the Spitzer relaxation time [7],√8 2πme e4 ni < Z 2 > ln Λeiνei =3mi (kB Te )3/2

(2.19)

where ni is the ion density and mi is the ion mass.

2.4.3

Model for collisional absorption

The total fraction of the collisional absorption of the laser light on its pathto the critical surface and the backward path after reflection is determinedfrom the Lambert-Beer law dIL /dx = −Ka IL as shown in Fig. 2.2.
( ∫ L)
ηa = 1 − exp −2
Ka dx(2.20)

02.4.4

Dependence of absorption on wavelength

The most outstanding feature of the collisional absorption is the strong
dependence on wavelength, since the laser beam is reflected on the criticalsurface. When the wavelength becomes a half, the critical density becomes
four times larger as seen from Eq. (2.1). Thus, the short wavelength light
can penetrate deeper, and enhance the collisional absorption given in Eq.

CHAPTER 2. LASER-TARGET INTERACTION

26(2.18). Finally, the deposited energy in the low density region is transported
to the high density plasma and continues creating ablation plasma on theablation surface, as shown in Fig. 2.3.

2.5

1D Lagrangian radiation hydrodynamic
code

2.5.1

Hydrodynamic equations

A 1D Lagrangian radiation hydrodynamic code [8, 9, 10, 11] named Star1D was developed by Sunahara [12, 13]. Star-1D is able to compute thehydrodynamic evolution of radiating plasma heated by laser pulses. It employs a one-fluid, two-temperature (Ti ̸= Te ) model, as shown in Fig. 2.3,
where Ti is the ion temperature. The plasma is assumed to be quasi-neutral
due to co-moving electrons (i.e., Zni = ne ). We solve the following equations.
Equation (2.21) is the equation of continuity, where ρ is the mass densityand v is the fluid velocity. Equation (2.22) is the equation of motion inplanar and spherical geometries in which the artificial viscosity q smoothens
the shock front. ps is the pressure. Here, the subscript s indicates either anion (i) or an electron (e).

dρ+ ρ∇ · v = 0
dt(2.21)dv
1= − ∇(pi + pe + q)dt
ρ(2.22)

Equations (2.23) and (2.24) are the energy equations for ions and electrons, where cVs is the specific heat, κs the thermal conductivity and α thethermal relaxation coeﬃcient, respectively.

CHAPTER 2. LASER-TARGET INTERACTION

27

ne
nc
x
Critical front
L
Laser
Figure 2.2: Normal incidence of the laser light on planar ablation plasma.

Low density
plasmaHigh density
plasma

Te
v
Laser

Solid
Ti
ne
Critical front Ablation front

Figure 2.3: Ion and electron temperature, ablation plasma velocity and electron density distributions in a laser-produced plasma.

CHAPTER 2. LASER-TARGET INTERACTION

28

ρcVi

dTi
= −pTHi ∇ · v + ∇ · (κi ∇Ti ) + α(Te − Ti ) + Qi
dt

(2.23)

ρcVe

dTe= −pTHe ∇ · v + ∇ · (κe ∇Te ) − α(Te − Ti ) + Qe
dt

(2.24)(
pTHs = Ts

∂ps
∂Ts

)
(2.25)
ρIn Eqs. (2.22)-(2.25), the variables ps , cVs , and pTHs are functions of ni
and Ts . They are obtained from the SESAME equation of state tables [14]based on the local thermodynamic equilibrium (LTE) model. The second
terms on the right-hand sides of Eqs. (2.23) and (2.24) represent the thermal
conduction due to ions [15] and electrons [16, 17], as shown in Eq. (2.26).

κi = √

2.925(kB Ti )5/2
πmi e4 < Z 2 >2 ln Λii√
16 2(kB Te )5/2
κe =,√ϕπ 5/2 me e4 Z ∗ ln Λei(2.26)
Z ∗ + 4.2
ϕ= ∗Z + 0.24

We used the flux-limited Spitzer-Härm model [18]（with a flux limit factorof 0.1 [19]）. The third terms in Eqs. (2.23) and (2.24) indicate the electronion temperature relaxation, and α(= 1.5kB ne νei /f ) is determined. f is the
reduction factor of the laser absorption due to the Langdon eﬀect [20], which
is important for long wavelength laser irradiation. Qi (= −q∇ · v) is theviscous heat, generated by conversion of the kinetic energy to the ion thermalenergy, due to the viscosity. Qe represents a laser energy absorption forelectrons and radiation. The laser energy deposited between the vacuum
plasma boundary (the expansion surface) and the critical surface is calculatedby ray tracing using 100 rays. We assumed the inverse bremsstrahlung tooccur in the laser energy absorption process. The calculated deposition ateach calculation mesh is P ∝ exp(−vg Ka ∆t), where P is the local laser

CHAPTER 2. LASER-TARGET INTERACTION

29

power, and ∆t is the time required for the laser ray to pass through eachmesh size at its group velocity.

2.5.2

Radiative transport equations

Equation (2.27) is the radiative transport equation, where Eν is the photon energy density for ν, Dν the diﬀusion coeﬃcient for ν, ην the emissivity forν, χν the absorption coeﬃcient for ν, and ν the frequency group, respectively.
Here, the emissivity and the opacity have [W/m3 /second/dν/steradian] and[1/m] in SI-unit, respectively. To study the radiation transport process, themulti-group diﬀusion approximation model [21] using a flux limit as shown in
Eq. (2.28) is applied. The photon energy hν is in the range of 0-1 keV, andthe range is divided into 100 bins for ν to solve Eq. (2.27) numerically. h is
the Planck constant. hν and the opacity κν (= χν /ρ) are evaluated using thecollisional radiative equilibrium (CRE) model. In CRE model, the plasma isassumed to be optically thin: plasma does not reabsorb its radiation energy.
d
ρ
dt

(

Eν
ρ

)
+ ∇ · (Dν ∇Eν ) = 4πην − cχν Eν

Dν = −

2.6c3χν + |∇Eν |/Eν

(2.27)(2.28)

Summary

In this Chapter, we described the energy deposition of high power beams
and a 1D Lagrangian radiation hydrodynamic code (named Star-1D) which
can compute the hydrodynamic evolution of radiating plasma heated by laser
pulses.
In the modified numerical simulation computing high energy laser injection onto a spherical target and the ablation plasma dynamics in a magnetic

CHAPTER 2. LASER-TARGET INTERACTION

30

field, the results of Star-1D are used as 3D distributions of ion super particle
for the particle position, expansion velocity, charge state and mass. They are
inputted into a 3D hybrid PIC code. To convert from a 1D Lagrangian meshmodel to a 3D Cartesian coordinate particle model, the plasma is assumed
to expand into the vacuum isotropically as the 1D meshes are elements along
the radius of multi-layer sphere. Since the number of ions in each shell is calculated by multiplying the shell volume by the ion density, it is easy that todivide the ion super particles into each shell. Finally, the physical quantitiesof the 1D meshes are distributed among the ion super particles.

Bibliography[1] Dawson, J. M., and Oberman, C. R., “High-Frequency Conductivity andthe Emission and Absorption Coeﬃcients of a Fully Ionized Plasma,”
Phys. Fluids, Vol. 5, Issue 5, 1962, pp. 517-524. doi: 10.1063/1.1706652[2] Atzeni, S., and Meyer-ter-Vehn, J., The Physics of Inertial Fusion:
Beam Plasma Interaction, Hydrodynamics, Hot Dense Matter, Oxford
University Press, London, 2004, Chap. 11.
[3] Kruer, W. L., The Physics of Laser Plasma Interactions, Westview
Press, Colorado, 2003, pp. 45-56.
[4] Shearer, J. W., “Eﬀect of Oblique Incidence on Optical Absorption of
Laser Light by a Plasma,” Phys. Fluids, Vol. 14, Issue 1, 1970, pp.
183-185. doi: 10.1063/1.1693272[5] Johnston, T. W., and Dawson, J. M., “Correct Values for HighFrequency Power Absorption by Inverse Bremsstrahlung in Plasmas”
Phys. Fluids, Vol. 16, Issue 5, 1973, pp. 722. doi: 10.1063/1.1694419[6] Skupsky, S., “Coulomb Logarithm for Inverse-Bremsstrahlung Laser Absorption,” Phys. Rev. A, Vol. 36, Issue 12, 1987, pp. 5701-5712. doi:
10.1103/PhysRevA.36.5701[7] Spitzer Jr., L., Physics of Fully Ionized Gases, 2nd ed., Dover Publications, New York, 2006, pp. 131-136.
31

BIBLIOGRAPHY

32[8] Sizyuk, V., Hassanein, A., Morozov, V., Tolkach, V., Sizyuk, T., andRice, B., “Numerical Simulation of Laser-Produced Plasma Devices
for EUV Lithography Using the Heights Integrated Model,” Numer.
Heat Transfer Part A Appl., Vol. 49, Issue 3, 2006, pp. 215-236. doi:
10.1080/10407780500324996[9] MacFarlane, J. J., Golovkin, I. E., and Woodruﬀ, P. R., “HELIOSCR: A 1-D Radiation-Magnetohydrodynamics Code with Inline Atomic
Kinetics Modeling,” J. Quant. Spectrosc. Radiat. Transfer, Vol. 99, Issue
1-3, 2006, pp. 381-397. doi: 10.1016/j.jqsrt.2005.05.031[10] Nagatomo, H., “Computational Simulation for Laser Driven Implosion,”
J. Plasma Fusion Res., Vol. 81, 2005, pp. 59-63. [In Japanese][11] Ohnishi, N., “Fundamentals of Radiation Hydrodynamics Simulation,”
J. Plasma Fusion Res., Vol. 88, Issue 9, 2012, pp. 467-474. [In Japanese][12] Sunahara, A., Nishihara, K., and Sasaki, A., “Optimization of Extreme
Ultraviolet Emission from Laser-Produced Tin Plasmas Based on Radiation Hydrodynamics Simulations,” Plasma and Fusion Res., Vol. 3,
Issue 43, 2008, pp. 1-5. doi: 10.1585/pfr.3.043[13] Sunahara, A., Sasaki A., Tanuma, H., Nishihara, K., Nishikawa,T., Koike, F., Fujioka, S., Aota, T., Yamamura, M., Shimada, Y.,
Nishimura, H., Izawa, Y., Miyanaga, N., and Mima, K., “Radiation Hydrodynamic Simulation of Extreme Ultra-Violet Emission from LaserProduced Tin Plasmas,” J. Plasma Fusion Res., Vol. 83, Issue 11, 2007,
pp. 920-926. [In Japanese][14] Lyon, S. P., and Johnson, J. D., “SESAME: The Los Alamos National
Laboratory Equation of State. Database,” Los Alamos National Laboratory LA-UR-92-3407, 1992.

BIBLIOGRAPHY

33[15] Braginskii, S. I., “Transport Processes in a Plasma,” Rev. Plasma Phys.,Vol. 1, Consultants Bureau, New York, 1965, pp. 250.
[16] Spitzer Jr., L., and Härm, R., “Transport Phenomena in a Completely
Ionized Gas,” Phys. Rev., Vol. 89, Issue 5, 1953, pp. 977-981. doi:
10.1103/PhysRev.89.977[17] Lee, Y. T., and More, R. M., “An Electron Conductivity Model for
Dense Plasmas,” Phys. Fluids, Vol. 27, Issue 5, 1984, pp. 1273-1286.
doi: 10.1063/1.864744[18] Malone, R. C., McCrory, R. L., and Morse, R. L., “Indications ofStrongly Flux-Limited Electron Thermal Conduction in Laser-Target
Experiments,” Phys. Rev. Lett., Vol. 34, 1975, pp. 721-724. doi:
10.1103/PhysRevLett.34.721[19] Azechi, H., “Concept of Laser Fusion.,” J. Plasma Fusion Res., Vol. 81,
2005, pp. 2-10. [in Japanese].
[20] Langdon, A. B., “Nonlinear Inverse Bremsstrahlung and HeatedElectron Distributions,” Phys. Rev. Lett., Vol. 44, Issue 9, 1980, pp.
575-579. doi: 10.1103/PhysRevLett.44.575[21] Minerbo, G. N., “Maximum Entropy Eddington Factors,” J. Quant.
Spectrosc. Radiat. Transfer, Vol. 20, Issue 16, 1978, pp. 541-545. doi:
10.1016/0022-4073(78)90024-9

Chapter 3
3D Hybrid PIC Code
3.1

Introduction

The 3D hybrid PIC code used in the present study for simulating the
behavior of charged particles in a magnetic field is based on the model proposed by Horowitz [1]. Ions and electrons are treated as individual superparticles and a magnetic hydrodynamic fluid, respectively. This approach isvalid when the system is dominated by ion dynamics [2, 3, 4, 5]. The motionof an ion is solved self-consistently, that is, ions aﬀect electromagnetic field in
a 3D Cartesian coordinates and move in the field. In addition, this code hassome assumptions about ions, electrons and the field. The mass of electron
is assumed to be zero (i.e., me =0) to eliminate the kinetic electron eﬀect
such as electron plasma oscillation and cyclotron motion. The plasma is also
assumed to be quasi-neutrality due to co-moving electrons (i.e., Zni = ne ).
Z, ni and ne is charge state of ion, ion density and electron density, respectively. In order to combine this code with Star-1D, this code was improved
to treat 3D distributions of ion super particle for the charge state and massin addition to the particle position and expansion velocity. The results ofStar-1D are inputted into the 3D hybrid PIC code. Here, we describe the34CHAPTER 3. 3D HYBRID PIC CODE

35

basic equations and the numerical algorithms used in this code.

3.2

Particles in electromagnetic field

3.2.1

PIC algorithm

In this code, the electromagnetic field is only calculated on the discrete
rectangular grid points as shown in Fig. 3.1. • and ◦ indicate a particle andgrid points, respectively. The field is linearly interpolated to the continuous
particle positions for giving the correct ion acceleration. After the ions are
moved, the ion density and ion current density redeposit to the eight nearest
grid points. The algorithm to avoid a self-force (i.e., a particle accelerates
itself.) is PIC algorithm [6, 7] by using the same weighting in both densityand field calculations. The first-order weighting of Eq. (3.1) reduces highfrequency noise and smooths the density and field variation, where V is thecell volume and Vn is the assignment volume. Equation (3.2) solves physicalcomponent ϕ at the particle position by using the triangle weighting.

δVn =

V =

Vn
V
8
∑

(3.1)

Vn

1

ϕ =ϕi+1,j,k × δV1 + ϕi,j,k × δV2+ϕi+1,j+1,k × δV3 + ϕi,j+1,k × δV4+ϕi+1,j,k+1 × δV5 + ϕi,j,k+1 × δV6+ϕi+1,j+1,k+1 × δV7 + ϕi,j+1,k+1 × δV8(3.2)

CHAPTER 3. 3D HYBRID PIC CODE

36

φi,j+1,k+1

φi,j,k+1

φi+1,j+1,k+1

φi+1,j,k+1

φi+1,j+1,kφi,j,kφi+1,j,kFigure 3.1: Particle and eight nearest grid points. • and ◦ shows particle
and grid point, respectively.

CHAPTER 3. 3D HYBRID PIC CODE

3.2.2

37

Digital filter

The digital filter [8] improves overall accuracy and reduces the noise by ignoring source terms at short wavelengths. The phenomena of high frequency
plasma are also neglected. The filter is limited to operate the eighteen nearest grid points as shown in Fig. 3.2. This filter is designated by the values(M, S, K), where M =middle, S=side and K=corner weights, as shown in
Eqs. (3.3), (3.4) and (3.5). New physical component ϕnew
i,j,k is obtained Eq.
(3.6) given by substituting a binomial eighteen-points filter (12, 4, 1) into
Eq. (3.3).

ϕnewi,j,k ←

cornerM ϕi,j,k + Sϕside
i,j,k + Kϕi,j,k
M + 6S + 12K

ϕside
i,j,k = ϕi−1,j,k + ϕi+1,j,k + ϕi,j−1,k + ϕi,j+1,k + ϕi,j,k−1 + ϕi,j,k+1

(3.3)

(3.4)

ϕcorneri,j,k = ϕi,j−1,k−1 + ϕi,j+1,k−1 + ϕi,j−1,k+1 + ϕi,j+1,k+1 + ϕi−1,j,k−1 + ϕi+1,j,k−1+ ϕi−1,j,k+1 + ϕi+1,j,k+1 + ϕi−1,j−1,k + ϕi+1,j−1,k + ϕi−1,j+1,k + ϕi+1,j+1,k(3.5)

1 side1
1
ϕi,j,k + ϕcorner
ϕnewi,j,k = ϕi,j,k +4
12
48 i,j,k(3.6)

CHAPTER 3. 3D HYBRID PIC CODE38

φi,j+1,k+1
φi-1,j,k+1
φi,j-1,k+1
φi-1,j,kφi+1,j,k+1

φi,j,k+1
φi-1,j+1,kφi,j+1,k

φi,j,kφi-1,j-1,kφi+1,j-1,kφi,j-1,kφi+1,j,kφi,j+1,k-1

φi-1,j,k+1
φi,j-1,k-1

φi+1,j+1,kφi,j,k-1

φi+1,j,k-1

Figure 3.2: Eighteen nearest grid points. The points indicate middle, side
and corner points.

CHAPTER 3. 3D HYBRID PIC CODE

3.3

39

Basic equations

The basic equations of this code are given below, where mi is the ion
mass, vi is the ion velocity, t is time, e is the elementary charge, E is theelectric field, B is the magnetic field, ve is the electron velocity and pe is theelectron pressure, respectively.

n i midvi= Zeni (E + vi × B) − R
dt

(3.7)

dx
= vi
dt
n e medve= 0 = −ene (E + ve × B) + R − ∇pe
dt

(3.8)

Equations (3.7) and (3.8) are the equations for particle and inertia-lessfluid motion, respectively. The term R on the equations is friction term
for the momentum of electron colliding with ion [9, 10]. This term solves
magnetic field diﬀusion on the boundary between plasma and magnetic field
and is given by

R = ne me (vi − ve )νei
= ne ηe(ji − ene ve )(3.9)

where νei is the electron-ion collision frequency and ji (= Zeni vi ) is the ion
current density. The resistivity η is usually taken as a scalar with a constantcoeﬃcient me νei /ne e2 and is a parameter depending on the electron-ion collision frequency. The third term on Eq. (3.8) is the relationship between thepressure and the temperature. Because this relationship is assumed to be an
isothermal or adiabatic while plasma expands, the electron temperature isassumed to be uniform for simplicity. ∇pe is given by

CHAPTER 3. 3D HYBRID PIC CODE40

∇pe = ∇ (ne kB Te )= ne ∇kB Te + kB Te ∇ne= kB Te ∇ne(3.10)

where kB is the Boltzmann constant and Te is the electron temperature. By
substituting Eqs. (3.9) and (3.10) into Eq. (3.8), we obtain the following
equation, where j(= ji + je ) is the current density and je (= −ene ve ) is theelectron current density.

0 = −ene (E + ve × B − ηj) − kB Te ∇ne(3.11)

Ampere’s law with Darwin approximation (See Appendix A) of Eq. (3.12)
neglects the eﬀects of high frequency phenomena such as electromagnetic
radiation and waves associated with electron inertia, that is, the electromagnetic field is treated in the low frequency approximation. The electron
velocity is expressed as follows, where µ0 is the permeability of vacuum.

∇ × B = µ0 j

(3.12)= µ0 (ji − ene ve )
1
ve =
ene

(
)
∇×B
ji −
µ0

(3.13)The electric field in the plasma region is computed by substituting Eq.
(3.13) into Eq. (3.11).

CHAPTER 3. 3D HYBRID PIC CODE

41

{
(
)
}
1
∇×B
0 = ene E +ji −× B − ηj + kB Te ∇ne
ene
µ0
{(
)
}
1
∇×B
E=−
ji −
× B + kB Te Z∇ni + ηj
Zeni
µ0
{(
)
}
1
∇×B
η
=−
ji −
× B + kB Te Z∇ni + ∇ × B
Zeni
µ0
µ0

(3.14)The partial diﬀerential equation for the time development of the magnetic
field is advanced by Eq. (3.14), Maxwell-Faraday’s law of Eq. (3.15) and thevector triple product of magnetic field (i.e., ∇×(∇×B) = −∇2 B+∇(∇·B)).

∇×E =−

∂B
1
=
∇×
∂t
Zeni

{(

∇×B
ji −
µ0

∂B
∂t

)

(3.15)}

× B + kB Te Z∇ni

+

η 2
∇B
µ0

(3.16)

After the advance of magnetic field is solved, the electric field can besolved. Then, the motion of ion is simulated by using leapfrog algorithm.
Here, η and Te is assumed to be zero for simplicity.

3.4

Leapfrog algorithm

The time integration scheme for solving partial diﬀerential equation isfollowed by the second-order leapfrog algorithm. The updates in time are
relatively straightforward. The calculation steps at the time step n (denotedby the superscript) are as follows. Figure 3.3 indicates the initial step (1, 2,
3, 4) and the cycle steps (5, 6, 7) of leapfrog algorithm. (The components ofion velocity and electric field are indicated in Appendix B.) The time step ∆t
is according to the Courant Friedrichs Lewy (CFL) condition as the speed oflight is replaced by the Alfven speed, since hybrid PIC code uses the implicit
method not the explicit method.

CHAPTER 3. 3D HYBRID PIC CODE

3.4.1

42

Initial steps

1. x0 , vi0 , Z and mi of each particle are inputted, where x is the particle
position. B 0 of each grid points is also inputted to this 3D hybrid PIC
code. PIC algorithm interpolates n0i , ji0 , Z 0 and m0i for each grid pointfrom x0 and vi0 . For calculating electromagnetic field, the grid points
in the calculation region are divided to the plasma and vacuum region.
in Fig. 3.3.)(

2. Since n0i , ji0 , Z 0 and B 0 are known, E 0 is obtained by using Eq. (3.17)
in Fig. 3.3.)

given from Eq. (3.14). (1
E =− n n
Z eni
n

{(
)
}
∇ × Bn
η
n
n
n
n
ji −
× B + kB Te Z ∇ni + ∇ × B nµ0µ0(3.17)

3. By substituting Eqs. (3.9) and (3.13) into Eq. (3.7), Eq. (3.18) as
n+ 12

implicit formula is obtained for solving vi

1

analytically. Thus, vi2 isobtained from vi0 , Z 0 , m0i , E 0 and B 0 by usinfg Eq. (3.19). (in Fig.

3.3.)

n
n dvi
mi
= Z ne

dt

(
)
η
n
n
n
n
E + vi × B − ∇ × Bµ0(3.18)

1v 2 − vi0
dvi0
= i
dt
∆t/2
n+ 12

4. The central diﬀerence of vi
12(3.19)

as Eq. (3.20) is divided into Eqs. (3.21)
1

and (3.22). Since vi is known, x 2 and x1 are obtained from Eqs. (3.21)
1

1

1

and (3.22), respectively. Then, PIC algorithm interpolates ni2 , ji2 , Z 2
1

and mi2 for each grid point. (in Fig. 3.3.)

CHAPTER 3. 3D HYBRID PIC CODE

n
0
xi
vi
ni, ji, Z, miE
B

n+1/2 n+1
1/2
2
1
4
4
1
3
1
4
2
5
5
1
5
5

43

n+3/2 n+2
6

6

6
6
7
7

Figure 3.3: Update steps in time of leapfrog algorithm.

CHAPTER 3. 3D HYBRID PIC CODE

441dxn+ 2
xn+1 − xn
n+ 1== vi 2dt∆t
1

xn+ 2 = xn +∆t n+ 12v
2 i

1

xn+1 = xn+ 2 +3.4.2(3.20)(3.21)∆t n+ 12v
2 i

(3.22)

Cycle steps
n+ 21Firstly, xn , B n , E n , vi

n+ 12

, ni

n+ 21

, ji

n+ 12

1

, Z n+ 2 , miand xn+1 are alreadyknown. The predictor corrector algorithm and the successive over-relaxation(SOR) algorithm are used to solve equations iteratively for E n+1 and B n+1 .
By using the new electromagnetic field, the particle positions can be advancedn+ 32

to n + 3/2. Finally, xn+1 , B n+1 , E n+1 , vi

n+ 32

, ni

n+ 32

, ji

3

n+ 32

, Z n+ 2 , miandxn+2 are obtained.

3.5

Calculation of electromagnetic field

3.5.1

Plasma-vacuum interface

In order to calculate the field quantities at the grid points, this code
determines by whether that the grid points is in the plasma or vacuum region
as shown in Fig. 3.4. The black area indicates the outside of calculation
region. The moving plasma-vacuum interface is not used to treat vacuum
region (See Appendix C). This method adequately treats variation in low
density plasma region to depart the two regions. The threshold of ion density
is empirically determined as 1016 cm−3 . Thus, the grid points which ion
density is larger than the threshold are in the plasma region. The other points
are determined to be the vacuum region. Then, in each region, the electricfield is calculated by predictor corrector and SOR algorithms, respectively.

CHAPTER 3. 3D HYBRID PIC CODE

45

y

x
vacuumplasmavacuumFigure 3.4: A schematic of the grid points used in the plasma and vacuum
region. The black area is the outside of calculation region.

CHAPTER 3. 3D HYBRID PIC CODE

3.5.2

46

Predictor corrector algorithm for plasma region

The time advance of the electromagnetic field in the plasma region isaccomplished by the predictor corrector algorithm. This algorithm is numerically stability, gives good energy conservation and is also used the arbitrary
plasma-vacuum interface. The calculation processes are as follows. (Figure
3.3 indicates the each step.)
Initial step
5. Equation (3.23) shows Maxwell-Faraday’s law. Equation (3.24) indicates the central diﬀerence of B n . By substituting Eq. (3.24) into
(3.23), Eq. (3.25) is obtained. The straightforward way to advance themagnetic field to level n + 1/2 is to use Eq. (3.25) and E n .

1

1

B n+ 2 − B n− 2dB n== −∇ × E n
dt
∆t
1

1

B n+ 2 + B n− 2
B =
2
n

1

B n+ 2 = B n −1(3.23)∆t∇ × En2(3.24)(3.25)

1

n+1
Since B n+ 2 is known, E n+ 2 is calculated from Eq. (3.17). Thus, Epredis predicted from the following equation Eq. (3.26).

E

n+ 12

n+1
+ En
Epred
=
2
1

n+1
= 2E n+ 2 − E n
Epred

n+1
is predicted from Eq. (3.27). (Bpred

in Fig. 3.3.)(3.26)

CHAPTER 3. 3D HYBRID PIC CODE

3

47

1

B n+ 2 − B n+ 2
= −∇ × E n+1
∆t
3

B

n+1

1

B n+ 2 + B n+ 2
=
2
1

n+1
Bpred
= B n+ 2 −∆tn+1∇ × Epred2(3.27)n+ 3

2
6. Because dvin+1 /dt is given as following equation Eq. (3.28), vi,predn+ 3

is predicted by using Eq. (3.18). xpred2 is also predicted by using
n+ 3

n+ 3

n+ 3

n+ 3

2
2
2
Eq. (3.21). PIC algorithm calculates ni,pred, ji,pred, Zpred2 and mi,pred.

Finally, xn+2
pred is obtained by using Eq. (3.22). (n+ 3

in Fig. 3.3.)

n+ 12

2
vi,pred− vi
dvin+1=dt∆t(3.28)

7. Some iteration is needed to give convergence because the magnetic field
has to been improved in advancing after the electric field advanced. (in Fig. 3.3.)

Cycle steps(a) The predicted fields at the time step n + 3/2 is computed with the1

newly E n+ 2 evaluated by using Eq. (3.17). Maxwell-Faraday’s
law as shown in Eq. (3.27) can be used to push magnetic field
n+ 3

n+1
. Equation (3.29) indicates Bpred2 by using Eq.
to the newly Bpred(3.25).

n+ 3

n+1
−
Bpred2 = Bpred∆tn+1∇ × Epred2(3.29)

CHAPTER 3. 3D HYBRID PIC CODE

48

n+ 3

Epred2 is predicted by using Eq. (3.17). The corrected electronic
field at the time step n + 1 is determined from Eq. (3.30).

n+ 3

E

n+1

1

Epred2 + E n+ 2=2

(3.30)(b) The corrected magnetic field at the time step n + 1 is determinedfrom Eq. (3.31).

1

B n+1 = B n+ 2 −∆t∇ × E n+12

(3.31)(c) The processes could be repeated to improve the accuracy as shownin Eq. (3.32).

|B n+1,new − B n+1,old |≤ 0.01|B n+1,old |3.5.3(3.32)

SOR algorithm for vacuum region

The rhs of Eq. (3.14) includes the term inversely proportional to the ion
density. In the vacuum region, this term becomes infinity and the calculation
diverges. To avoid this problem, the Laplace equation ∇2 E = 0 is solved toobtain the electric field in the vacuum region. Equation (3.33) is the central
diﬀerence of the Laplace equation, where ∆p is the mesh width and Ep is
the electric field for each direction. The subscription p indicates x or y or z
direction.

CHAPTER 3. 3D HYBRID PIC CODE49

Epi−1,j,k − 2Epi,j,k + Epi+1,j,k∆x2+

Epi,j−1,k − 2Epi,j,k + Epi,j+1,k∆y 2+

Epi,j,k−1 − 2Epi,j,k + Epi,j,k+1=0∆z 2(3.33)

When ∆x is equal to ∆y, the newly electric field is taken from Eq. (3.33).

Ep∗oldi,j,k =
+

∆z 2
old
old
old
{Epoldi−1,j,k + Epi+1,j,k + Epi,j−1,k + Epi,j+1,k }4∆z 2 + 2∆x2
∆x2
old
{Epoldi,j,k−1 + Epi,j,k+1 }4∆z 2 + 2∆x2

(3.34)The new electric field is predicted by using Eq. (3.35) and SOR method.

old
∗old
old
Epnew
i,j,k = Epi,j,k + α(Epi,j,k − Epi,j,k )(3.35)

Equation (3.36) indicates the acceleration coeﬃcient for the rectangulararea which is usually set between 0 and 2, where L is the represent length.
Here, α is determined to be 1. This is called Gauss-Seidel method.

α=2
1 + sin(π/L)(3.36)The processes could be repeated to improve the accuracy as shown in Eq.
(3.37).

new
Ep∗new
i,j,k − Epi,j,k≤ 0.01old
Ep∗oldi,j,k − Epi,j,k(3.37)

CHAPTER 3. 3D HYBRID PIC CODE

3.6

50

Summary

In this Chapter, we described the basic equations and the numerical algorithms of the 3D hybrid PIC code used in the present study for simulating
the behavior of charged particles in a magnetic field.

Bibliography[1] Horowitz, E. J., Shumaker, D. E., and Anderson, D. V., “QN3D: A
Three-Dimensional Quasi-neutral Hybrid Particle-in-Cell Code with Applications to the Tilt Mode Instability in Field Reversed Configurations,” J. Comput. Phys., Vol. 84, Issue 2, 1989, pp. 279-310. doi:
10.1016/0021-9991(89)90234-9[2] Harned, D. S., “Quasineutral Hybrid Simulation of Macroscopic Plasma
Phenomena,” J. Comput. Phys., Vol. 47, Issue 3, 1982, pp. 452-462. doi:
10.1016/0021-9991(82)90094-8[3] Winske, D., Yin, L., Omidi, N., Karimabadi, H., and Quest, K. B.,“Hybrid Simulation Codes: Past, Present, Future-A Tutorial,” Space
Plasma Simulations, Vol. 615, 2003, pp. 136-165.
[4] Kajimura, Y., and Nakashima, H., “Verification of Hybrid Particlein-Cell Simulation Model for Advanced Plasma Propulsions: Magneto
Plasma Sail and Magnetic Nozzle for Laser Fusion Rocket,” J. SpaceTechnol. Sci., Vol. 25, Issue 2, 2011, pp. 34-54.
[5] Matsumoto, M., Kajimura, Y., Funaki, I., and Shinohara, I., “Development of Hybrid Particle-In-Cell Simulation Code for Research of Magnet
Plasma Sail,” JAXA Research and Development Memorandum RM-11019, 2012. [In Japanese]51

BIBLIOGRAPHY

52[6] Birdsall, C. K., and Langdon, A. B., “Plasma Physics via Computer
Simulation,” Francis & Taylor, New York, 2004, Chapter 2.
[7] Naitou, Y., “Basic Theory of Particle Simulation,” J. Plasma FusionRes., Vol. 74, Issue 5, 1998, pp. 470-478. [In Japanese][8] Birdsall, C. K., and Langdon, A. B., “Plasma Physics via Computer
Simulation,” Francis & Taylor, New York, 2004, Appendix C.
[9] Mankofsky, A., Sudan, R. N., and Denavit, J., “Hybrid Simulation of
Ion Beams in Background Plasma,” J. Comput. Phys., Vol. 70, Issue 1,
1987, pp. 89-116. doi: 10.1016/0021-9991(87)90003-9[10] Bagdonat, T., and Motschmann, U., “3D Hybrid Simulation Code Using
Curvilinear Coordinates,” J. Comp. Phys., Vol. 183, Issue 2, pp. 470485. doi: 10.1006/jcph.2002.7203

Chapter 4
Direct measurement of impulse
4.1

Introduction

Direct measurement of impulse due to the interaction between a diamagnetic current in a laser-produced plasma and a magnetic field produced bycylindrical neodymium permanent magnets were performed. The impulse
was directly measured by using a pendulum type thrust stand for various
laser pulse energies at EUV database and GXII laser facilities of ILE in
the collaboration with Osaka University. [1, 2, 3, 4]. In addition to theexperiments, to confirm the validity of the modified numerical simulation
computing high energy laser injection onto a spherical target and the ablation plasma dynamics in a magnetic field, the numerical simulations were
performed for comparing the measured impulse with the calculated impulse.

53

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

4.2

54

Dependence of impulse on laser pulse duration

4.2.1

Experimental setup for EUV facility

Figure 4.1 is a photo of the vacuum chamber in the EUV database laser
facility. The experimental setup for EUV facility is shown in Fig. 4.2. A
laser-produced plasma is created by a 1064 nm single-beam neodymium:yttrium aluminum garnet (Nd:YAG) laser focusing into a 500 µm wide indiameter polystyrene ([-CH2 -CH(C6 H6 )-]n ) spherical target. The target issuspended at the end of a carbon fiber attached to a glass rod to reduce theinfluence of ablation plasma generated by the glass rod. The laser beam spot
diameter is the same as the target diameter by using a lens having F-numberof 20. Six cylindrical neodymium permanent magnets connected in the z-axis
are used for producing the magnetic field. The permanent magnet itself is
16 mm in diameter and 10 mm long. The surface magnetic flux density ofeach individual permanent magnet is 0.47 T. The magnetic field strength at
the target is calculated to be 0.12 T, as shown in Fig. 4.3. This value wasconfirmed by a gauss meter. The permanent magnets are arranged so as not
to disturb the incoming laser beam. The distance between the surface of thepermanent magnets and the target is 11 mm. The angle between the z-axisof the permanent magnets and the laser beam is also 45 ◦ .

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

Chamber

Laser

Figure 4.1: Vacuum chamber in the EUV facility.

ChamberThrust stand11 mm z-axis of magnets
45 degrees
Nd:YAG laser

Target
235 mmFigure 4.2: Experimental setup for EUV facility.

55

Magnetic field strength, T

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

1.20

1.4

Magnetic field strength on the z-x plane, T

1.2

1.16
1.11

1.0

1.06
1.01

0.8

0.955

0.6

0.905
0.854

0.4

0.804

0.12 T

0.2
0.0
-50

56

0.754
-40

-30

-20

-10

0

10

20

30

40

50

z, mm

50

0.704
0.653
0.603
0.553

Magnetic field strength on the z-x plane, T

0.503

40

0.452

x, mm

0.402
30

0.352
0.302
0.251

20

0.201
0.151

10
0
-50

0.101
0.0503
-40

-30

-20

-10

0

10

20

30

40

50

0

z, mm

Figure 4.3: Calculated magnetic field strength produced by the permanent
magnets at the z-axis at x = y = 0 and in the z − x plane at y = 0. The
magnetic field strength at the origin is 0.12 T.

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

4.2.2

57

Pendulum type thrust stand for EUV facility

The impulse is measured by means of a pendulum type thrust stand for
EUV facility, as shown in Fig. 4.4. The stand is composed of the permanent
magnets, knife edges, counter weights, light emitting diode (LED) displacement sensors, and an oil damper. The magnets are installed on the topof the pendulum. The pendulum is supported by the knife edges mounted
on side posts. The counterweight is placed on the pendulum such that thecenter of gravity lies in the pendulum’s axis of rotation. This ensures thatthe pendulum swings easily even for a small impulse. The natural periodof this pendulum is 1.2 second. The LED sensors detect the deflection of
the pendulum. To decrease the error in the impulse measurement, i.e., tooﬀset the mechanical noise experienced by the sensors, the LED sensors areplaced face-to-face and supplied by a single power source. The oil damper
is installed to suppress mechanical noise induced by background vibrations,which are critical problems for small impulse measurements. To prevent thestand from momentum transfer due to the direct collisions between the expanding plasma and the permanent magnets, a ceramic plate is placed at a
distance of 1 mm from the magnets to shield from the plasma. This plateis connected to a grounded vacuum chamber. The stand is calibrated by using small hammer which takes impulse to the stand. The impact is directly
measured by a force transducer attached to the counterweight. The force
transducer records the impulsive force as a function of time; integrating thisyields the applied impulse. It is assumed that the impulse on the force transducer is proportional to the amplitude of the actual response of the LED
sensors. This assumption allows for a linear curve fit to the calibration data,
as shown in Fig. 4.5. The correlation coeﬃcient (R2 ) is determined to be
0.99922. The experimental data are converted into an impulse force using
the conversion factor 70 µNs/V found from the calibration.

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

58

Figure 4.4: Schematic of the pendulum type thrust stand for EUV facility.
(vertical sectional view)100

y = 70 x
2

Impulse,Ns

R= 0.99922

50

0
0.0

0.5

1.0

1.5

Actual response of LED sensors, V

Figure 4.5: Applied calibration force as a function of the signals of LED
sensors and the liner curve fitting to the calibration data.

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

4.2.3

59

Results and discussion for EUV facility

Figure 4.6 shows the impulses measured on the pendulum type thruststand for EUV facility as a function of incident laser energies. The two data
sets show the measured impulses for single laser shot with laser pulse durationof 9 ns and 3 ns, respectively. Every data point plotted is the average result
over more than three shots; the error is estimated for the associated standard
deviation. In every shot, the incident laser energy and laser pulse duration
have some error from the setting value. The eﬀects of them are considered
to be the main cause of the error bars on the measured impulse.
When a non-magnetized iron cylinder is installed on the pendulum typethrust stand instead of the permanent magnets, the measured impulse is zero.
The result shows the interaction between the diamagnetic current in the laserproduced plasma and the magnetic field of the permanent magnets. Thus,
the magnetic thrust chamber system for LFR is demonstrated to produce
impulse experimentally.
The impulse is 2.0 µNs at the laser pulse duration of 9 ns, and it is 1.5µNs at the laser pulse duration of 3 ns, when the incident laser energy is 0.7
to 0.8 J. The diﬀerence between 9 ns and 3 ns can be considered as follows.
When charged particles of the laser-produced plasma are pushed back by themagnetic field before they collide to the shield plate, the estimated impulse isassumed to be the product of the plasma ablation mass and the ion acoustic
speed Cs , i.e., SτL ṁCs . S is the plasma ablation area of the target by focusing
the incident laser beam. τL is the laser pulse duration. The plasma expansion
velocity into a vacuum is approximately equal to the ion acoustic speed. The
plasma ablation mass per unit area per unit time ṁ and the ion acoustic
speed are expressed as follows [5], where me is the electron mass and mi is
the ion mass. The critical density nc is 9.8 ×1026 m−3 .

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE)1/3 √
2/3 (nc
3Ka IL
1/3ṁ =(1 + Z)me mi
Z
4f
√
(
)1/3
1/3
(1 + Z)me
3Ka IL
Cs =
4f nc
mi

60

(4.1)
(4.2)

Here, the charge number of ion Z is assumed to be 4 (C4+ is assumed tobe the primary ion in the plasma), the laser absorption coeﬃcient Ka and theflux limit factor f is empirically assumed to be 1 and 0.08 [5], respectively. IL
is the incident laser intensity, which has the relationship IL = EL /SτL , where
EL is the incident laser energy. In the experiments, the plasma ablation area
is assumed to be the surface of target where the solid angle is 0.4π steradian,as shown in Fig. 4.7. Thus, the plasma ablation area is fixed to be 7.9 ×10−8
m2 .
The curves in Fig. 4.6 represent the estimated impulse, i.e., SτL ṁCs .
The solid and dashed lines are shown for laser pulse duration of 9 ns and 3
ns, respectively. According to Eqs. (4.1) and (4.2), the plasma ablation massand expansion velocity increase with increasing incident laser intensity, that
1/3

is, ṁ ∝ IL

and Cs ∝ IL . When EL is constant (τL ∝ IL−1 ), the impulse is
1/3

−1/3

proportional to SτL ṁCs ∝ IL

1/3

, that is, τL . Thus, the impulse increases

theoretically with the laser pulse duration. The experimental dependence of
impulse on laser pulse duration shown in Fig. 4.6 is in good agreement withthe estimated values.

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

61

4.5

4.0

9 ns
3 ns

3.5

Impulse,

Ns

3.0

2.5

2.0

1.5

1.0

0.5

0.0
0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1.8

2.0

Incident laser energy, J

Figure 4.6: Impulses measured on the pendulum type thrust stand as afunction of incident laser energies for laser pulse duration of 3 ns and 9 ns.

0.4π sr
8 mm
Magnets

11 mm

TargetS
Plasma motion

Figure 4.7: Impulse model of the product of plasma ablation mass and ion
acoustic speed.

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

4.3

62

Dependence of impulse on laser pulse wavelength and incident laser energy

4.3.1

Experimental setup for GXII facility

A plasma is generated by focusing a 1.3 ns laser beam onto a 1.5 mm
diameter spherical polyacetal ([-CH2 O-]n ) target at GXII facility. The target
is suspended at the end of a carbon fiber attached to a glass rod to reducethe influence of ablation plasma generated by the glass rod. Figure 4.8 isa picture of the vacuum chamber in the GXII facility. The experimental
setup is shown in Fig. 4.9. The glass laser beam is focused using a lens (Fnumber: 3) so that its diameter is equal to that of the target. A cylindricalneodymium permanent magnet is used to generate the magnetic field. The
maximum energy product of the magnet is 0.382 mJ/mm3 . The permanent
magnet is 50 mm diameter and 40 mm long. It is located at where it does
not obstruct the laser beam. The distance between the permanent magnet
surface and the target center is 33 mm and the angle between the permanent
magnet axis and the laser beam is 66.5◦ . The magnetic field strength at
the target is calculated to be 0.1 T, as shown in Fig. 4.10. This value wasconfirmed by a gauss meter.

4.3.2

Pendulum type thrust stand for GXII facility

The impulse is measured using a pendulum type thrust stand for GXII
facility (See Fig. 4.9). The pendulum type thrust stand for GXII facility islargely similar to the stand for EUV facility and consists of the permanent
magnet, bearings, a counterweight and LED displacement sensors. The pendulum is supported by the bearings, which are mounted on side posts. The
pendulum has a natural period of 1.3 s. This is short enough to reduce themechanical noise induced by background vibrations, which are problematic

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

63

when measuring small impulses. A ceramic plate connected to a groundedvacuum chamber is placed 2 mm from the magnet to shield it from theplasma. Figure 4.11 indicates the calibration data. The R2 coeﬃcient is
0.99953. The experimental data are converted into forces using a conversion
factor of 240 µNs/V, which is obtained from the calibration.

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

64

Laser beams

Chamber

Figure 4.8: Vacuum chamber in the GXII facility.

Figure 4.9: Experimental setup & Schematic of the pendulum type thruststand for GXII facility.

Magnetic field strength, T

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE
0.8

0.960
0.910

0.6

0.859
0.809

0.4

0.759
0.709

0.2

0.658

0.1 T

0.608
0.0
-50

-40

-30

-20

-10

0

10

20

30

40

0.558

50

0.508

z, mm

50

0.457

Magnetic field strength on the z-x plane, T

0.407

40

x, mm

65

0.357
0.307

30

0.256
0.206

20

0.156
0.106

10

0.0553
0
-50

-40

-30

-20

-10

0

10

20

30

40

0.00500

50

z, mm

Figure 4.10: Calculated magnetic field strength produced by the permanent
magnet at the z-axis at x = y = 0 and in the z − x plane at y = 0. The
magnetic field strength at the origin is 0.1 T.

800

700

y = 240 x

Impulse,Ns

R2= 0.99953

600

500

400

300
1.5

2.0

2.5

3.0

Actual response of LED sensors, V

Figure 4.11: Applied calibration force as a function of the signals of LED
sensors and the liner curve fitting to the calibration data.

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

4.3.3

66

Results and discussion for GXII facility

Figure 4.12 shows the impulses measured on the pendulum type thruststand as a function of the incident laser energy. The three data sets show
impulses for single laser shots with laser frequencies of ω (λL =1,053 nm), 2ω
(λL =527 nm) and 3ω (λL =351 nm), where λL is the laser pulse wavelength.
Impulses of 1.3, 5.0, and 6.4 mNs are generated at laser frequencies of ω,2ω and 3ω for incident laser energies of 548, 568 and 550 J, respectively.
Figure 4.12 shows that impulses 3ω is the largest, followed in order by thoseof 2ω and ω. In addition, it shows that the impulse increases with increasing
incident laser energy.
The impulse for a specific laser wavelength as shown in Fig. 4.12 isestimated using the estimation model similar to that of EUV facility. Thelines in Fig. 4.12 are the estimated impulse curves obtained by assuming thatthe impulse is equal to the product of the plasma ablation mass and the ion
acoustic speed. All the plasma is pushed back by the magnetic field before
it collides with the plate of the stand. In this case, the plasma ablation areaof the target is constant at 3.5 × 10−6 m2 , when the area is assumed to beequal to the target surface for a solid angle of 2π since the target is irradiatedby a single laser beam through a lens with an F-number of 3.
The critical density for each laser frequency is 1.0 × 1027 m−3 for ω, 4.0
× 1027 m−3 for 2ω, and 9.0 × 1027 m−3 for 3ω. According to Eqs. (2.1), (4.1)
and (4.2), the product of the plasma ablation mass and the expansion velocity
increases with increasing critical density and incident laser energy; that is,−2/3ṁCs ∝ λL

2/3

and IL . Since τL is constant, the impulse is proportional−2/3

to SτL ṁCs ∝ λL

2/3

and EL . Thus, the estimated impulse increases withdecreasing laser pulse wavelength and increasing incident laser energy. The
measured impulses, which are indicated by the plots in Fig. 4.12. The results
show that the same trends as the estimation model.

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

67

10

2

Impulse, mNs

3
1

0.1

0.01
10

100

1000

Incident laser energy, J

Figure 4.12: Impulses measured on the pendulum type thrust stand as afunction of incident laser energies for laser pulse wavelength.

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

4.4

68

Comparison of measured and calculated
impulse4.4.1

Numerical simulation with and without considering creation process of laser-produced plasma

In the original numerical simulation modeling the experiment of EUV
facility, Maxwellian parameters is used as the initial plasma parameters inputted into the 3D hybrid PIC code. The laser-produced plasma is assumed
to be uniform. The plasma radius is set to 9 mm because the plasma expands isotropically. The plasma mass N mi is equal to the target mass (6.9
× 10−8 kg) under the assumption that all the target is ionized. The numberof ion super particles N is 1,000,000. The charge state of the ions < Z >is assumed to be 1.25 (=C4+ +H+ /4), since laser-produced plasma is weakly
ionized plasma. <> indicates averaging over the ion species. The atomic
mass mi is 6.5 amu. The ion density is 2.1 × 1018 cm−3 since the plasma
radius is 9 mm and the mass density of polystyrene is 1.06 g/cm3 . The ion
super particle expansion velocity distribution in the plasma is assumed to be
Maxwellian, as shown in Eq. (4.3) and Fig. 4.13, where v is the expansion
velocity, kB is the Boltzmann constant and Ti is the ion temperature. The
plasma kinetic energy is 0.3 J for incident laser energy of 1.0 J since theabsorption rate is assumed to be 0.3. This is our empirical prediction fromthe target and laser conditions. The characteristic velocity vp is 2.4 km/s.
(
f (v) = 4πv

2

3
N kB Ti = 0.3 J
21kB Ti = mi vp22

mi
2πkB Ti

) 32

(
exp

−mi v 2
2kB Ti)(4.3)

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

69

Furthermore, to improve the numerical simulation of the magnetic thrustchamber system for LFR, the modified numerical simulation employing a 3D
hybrid PIC code and Star-1D, which accounts for the creation process of alaser-produced plasma, was performed. Table 4.1 lists the characteristic initial plasma conditions inputted into the 3D hybrid PIC code for the original
and modified cases.
Numerical simulation of a 3D hybrid PIC code is performed for two initial conditions with and without considering the creation process of a laserproduced plasma. The calculation model of the 3D hybrid PIC code is illustrated in Fig. 4.14. This model is an 100∆x × 100∆x × 100∆x cube,
where ∆x=1 mm is the mesh width. The calculation time is 100,000∆t (5µs), where ∆t=50 ps.

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

3.5x104<
3.0x10

Number

2.5x102.0x101.5x101.0x105.0x1070

4v >=2.4 km/s
2

4

4

4

4

3

0.0
0.0

2.5

5.0

7.5

Ion super particle velocity, km/s

Figure 4.13: Maxwell super particle velocity distribution. The characteristic
velocity is assumed to be vp =2.4 km/s.

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

71

Table 4.1: Plasma conditions for original and modified cases.
Original case

Modified case

Plasma radius, mm

96Plasma mass, kg

6.9 × 10−8

2.5 × 10−8

Atomic mass, amu

6.5

6.5

<Atomic number>

1.25

0.36

Ion density, cm−3

2.1 × 1018

non-uniform

Kinetic energy, J

0.30

0.58

Ion expansion velocity

Maxwell distribution

non-uniform

Initial max. velocity, km/s

7.5

117

x

Magnets

11 mm

PlasmazyRegion : 100 x 100 x 100 mm
Meshes: 100 x 100 x 100

Figure 4.14: Cubic calculation model of 3D hybrid PIC code.

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

4.4.2

72

Results and discussion for numerical simulations

The modified numerical simulation with considering the creation processof a laser-produced plasma combines Star-1D and the 3D hybrid PIC code
described in Chaps. 2 and 3, respectively. Figure 4.15 shows the unit vector
distributions of the diamagnetic current calculated by the modified numericalsimulation at calculation times of 75 ns, 1.325, 2.575 and 3.825 µs. At 75 ns,
the plasma is expanding isotropically. Then, with updating calculation time,
the current rotates clockwise. It is indicated that the plasma produced the
diamagnetic cavity as shown in Fig. 1.2.
The measured impulse is 2.42 ± 0.23 µNs at a laser pulse energy of 1.05
J and a laser pulse duration of 9 ns. In the numerical simulations, theoriginal numerical simulation used only a 3D hybrid PIC code, which employs
the initial plasma conditions without considering the creation process of alaser-produced plasma; it computed the impulse to be 72.7 µNs. There isa large diﬀerence between these results. On the other hand, the calculated
impulse of the modified numerical simulation is 10.3 µNs. The impulse of themodified numerical simulation does not agree with the value of the previous
experiment. This is considered that the ablation plasma on the opposite
surface of laser injection is not calculated well. However, the value is closer
to that of the experiment than that of the original numerical simulation. This
indicates that the creation process of a laser-produced plasma (initial plasma
condition) is important when investigating plasma behavior in a magnetic
field by numerical simulations. This is considered to be because the diﬀerent
ion expansion velocities calculated by the two numerical simulations aﬀect theimpulse. Consequently, the modified numerical simulation can describe theablation plasma dynamics in the magnetic field but the original numericalsimulation does not. This confirms the validity of the modified numericalsimulation.

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

0.01

73

0.02

0.015

0.005

0.01

y, m

y, m

0.005

0

0

-0.005-0.005-0.01-0.015

-0.01
-0.01

-0.005

0
x, m

0.005

-0.02
-0.02

0.01

-0.015

0.02

0.02

0.015

0.015

0.01

0.01

0.005

0.005

0

-0.005

-0.01

-0.01

-0.015-0.015-0.015-0.01-0.005

0
x, m

0
x, m

0.005

0.01

0.015

0.02

0.01

0.015

0.02

0

-0.005

-0.02
-0.02

-0.005(b) 1.325 µs

y, m

y, m

(a) 75 ns

-0.01

0.005

0.01

0.015

0.02

-0.02
-0.02

(c) 2.575 µs-0.015-0.01-0.005

0
x, m

0.005

(d) 3.825 µs

Figure 4.15: Unit vector distributions of diamagnetic current calculated bymodified and original numerical simulation.

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

4.5

74

Summary

In this Chapter, we described the direct measurement of impulse due tothe interaction between a diamagnetic current in the laser-produced plasma
and a magnetic field produced by the cylindrical neodymium permanent magnet.
In the experiments, the impulse was measured directly by means of a pendulum type thrust stand for various incident laser energies at EUV databaseand GXII laser facilities of ILE in the collaboration with Osaka University.
The pendulum type thrust stand was fabricated for measuring a small impulse. Without a magnetic field, the measured impulse was zero. On theother hand, with a magnetic field, impulse was generated. These results
show that the system for generating impulse works as theory tells. Thus, thesystem was demonstrated experimentally for eventual use in a space propulsion system. The experiments also found that the dependences of the impulseon the laser pulse duration, laser wavelength and incident laser energy in thesystem. Thus, for producing large impulse, it is revealed that longer pulse
durations, shorter wavelengths and higher energies are better.
In the numerical simulations, in order to compare the measured impulseof EUV facility with the calculated ones, the original and modified numericalsimulations were performed. The modified numerical simulation combining
Star-1D and the 3D hybrid PIC code computes high energy laser injection
onto a spherical target and the ablation plasma dynamics in a magnetic field.
The original numerical simulation does not consider the creation process ofa laser-produced plasma. The measured impulse was 2.42 ± 0.23 µNs at alaser pulse energy of 1.05 J and a laser pulse duration of 9 ns. The calculated
impulse of the original case was 72.7 µNs. There is a large diﬀerence between
these results. The calculated impulse for the modified case was 10.3 µNs.
From these values, the calculated impulse of the modified case is closer toCHAPTER 4. DIRECT MEASUREMENT OF IMPULSE

75

the measured impulse than the original case. This indicates that the creation
process of a laser-produced plasma is important factor for the numericalsimulation.
The achievements in this Chapter are as follows.
Experiments• By using cylindrical neodymium permanent magnet and non-magnetizediron cylinder, impulse due to the interaction between a diamagnetic
current in a laser-produced plasma and a magnetic field was generated
and the magnetic thrust chamber system for LFR was demonstrated
experimentally to work as theory tells.
• The impulse for the duration of 9 ns was larger than that of 3 ns. The
impulse increases, as the laser pulse duration becomes longer.
• The impulse for the wavelength of 3ω is the largest, followed in orderby that of 2ω and ω. The impulse increases, as the laser wavelength
becomes shorter.
• The impulse increases with increasing the incident laser energy.
• In order to produce large impulse, longer pulse durations, shorter wavelengths and higher energies enlarge the impulse.
Numerical simulations• The measured impulse of EUV facility was compared with the calculated impulse.
• The modified numerical simulation with considering the creation process of a laser-produced plasma as the initial plasma conditions of the
3D hybrid PIC code was performed.

CHAPTER 4. DIRECT MEASUREMENT OF IMPULSE76• The modified numerical simulation can describe the ablation plasma
dynamics in the magnetic field but the original numerical simulation
does not. This confirms the validity of the modified numerical simulation.

Bibliography[1] Akira Kakami, “Study of pulsed plasma thrusters by using liquid propellants,” Ph.D. Dissertation, University of Tokyo Graduate School ofEngineering, 2002. [In Japanese][2] Koizumi, H., Kakami. A., Komurasaki, K., and Arakawa, Y., “A ThrustStand for Low-Power Electric Propulsion,” Journal of the Japan Society
for Aeronautical and Space Sciences, Vol. 51, Issue 593, 2003, pp. 270275. [In Japanese][3] Koizumi, H., Komurasaki, K., and Arakawa, Y., “Development of thruststand for low impulse measurement from microthrusters,” Rev. Sci. Instrum., Vol. 75, Issue 10, 2004, pp. 3185-3190. doi: 10.1063/1.1790568[4] Tahara, H., “Dependence of Solid Propellant Species on Performance ofSmall Laser Ablative Thrusters,” J. IAPS, Vol. 16, Issue 1, 2008, pp.
27-32. [In Japanese][5] Azechi, H., “Concept of Laser Fusion.,” J. Plasma Fusion Res., Vol. 81,
2005, pp. 2-10. [in Japanese].

77

Chapter 5
Comparison of ion expansion
velocity
5.1

Introduction

In addition to the direct measurement of the impulse, we measured theion expansion velocity of a laser-produced plasma in a magnetic field. The
magnetic field is produced by an edgewise coil at EUV database laser facilityof ILE in the collaboration with Osaka University. To confirm the validity ofthe modified numerical simulation with considering the creation process of alaser-produced plasma, the numerical simulations with and without considering the creation process were performed for comparing with the experimental
result.

78

CHAPTER 5. COMPARISON OF ION EXPANSION VELOCITY

5.2

79

Experiment for measuring ion expansion
velocity

5.2.1

Experimental setup

Figure 5.1 shows the experimental setup. A laser-produced plasma iscreated by focusing a 1064 nm Nd:YAG laser beam (laser pulse duration: 9ns; incident laser energy: 4 J) onto a 500-µm-diameter spherical polystyrene
target. The laser beam is focused on the central axis of the coil by a lens
(F-number: 20) so that its diameter is equal to that of the target. To reducethe measurement error caused by the ablation plasma generated by the glassrod, the target is suspended on the end of a carbon fiber attached to a glassrod. A pulsed power supply and an edgewise coil are used to generate apulsed magnetic field. The coil (outer diameter: 60 mm; inner diameter: 50
mm; length: 10 mm) that generates the magnetic field consists of 10 turnsof a rectangular copper wire (width: 5 mm; thickness perpendicular to thecoil axis: 1 mm) as shown in Fig. 5.2. The wire is coated with a polyimide
thin film, as an insulator, whose breakdown voltage is 10 kV. The edgewise
coil has a high lamination factor, a small temperature diﬀerence between theinside and outside of the coil, and a high-frequency response. The coil is set
25 mm from the origin as shown in Fig. 5.1.

CHAPTER 5. COMPARISON OF ION EXPANSION VELOCITYx80

330 mmCoil

Chargecollector
14

Plasma 60 mm
Laser

z

y
50 mm

155 mm
10 mm 25 mm

3 mm

Figure 5.1: Experimental setup & Cubic calculation model of 3D hybrid PIC
code.

Figure 5.2: Edgewise coil for generating a pulsed magnetic filed.

CHAPTER 5. COMPARISON OF ION EXPANSION VELOCITY

5.2.2

81

Charge collector with a Faraday cup

The ion beam current of a laser-produced plasma is measured using acharge collector with a Faraday cup. The distance between the charge collector and the origin is 155 mm. The angle between the charge collector
and a single laser beam is 14◦ as shown in Fig. 5.1. Figure 5.3 shows the
charge collector that consists of a Faraday cup [1, 2, 3] and the currentto-voltage (I-V ) converter. The charge collector in the vacuum chamber isbiased to -50 V to repel electrons. Four resistors with a resistance of 50 Ωeach are used for impedance matching and signal denoising through the highpass filter. The Faraday cup reabsorbs secondary electrons, making it very
suitable for absolute measurement. This is because that almost none of thesecondary electrons generated in the Faraday cup escape from the chargecollector through the 2-mm-diameter hole, as shown in Fig. 5.3. To prevent
electrons in the plasma from entering the Faraday cup, the hole is covered
with a mesh which has 500 lines/inch with a line diameter of 0.025 mm andan ion penetration rate of 25%.

CHAPTER 5. COMPARISON OF ION EXPANSION VELOCITY

1 nF

Hole
Mesh
Faraday cup

1 MΩ
50 V

50 Ω
Oscilloscope

Charge collector

82

Chamber
Figure 5.3: Charge collector with a Faraday cup and an I-V converter. The
2-mm-diameter hole is covered with a mesh to prevent electrons in the plasma
from entering the Faraday cup.

CHAPTER 5. COMPARISON OF ION EXPANSION VELOCITY

5.2.3

83

Pulsed power supply for generating high pulsed
current

Pulse Forming Network
It is diﬃcult for neodymium permanent magnet to generate strong magnetic field flexibility, since the surface magnetic flux density of the neodymium
permanent magnet is within 0.5 T. On the other hand, coil and high-current
can make strong magnetic field than the permanent magnets. For producing
high pulsed current, pulsed power supply is superior to steady power supply
under the thermal limitation of circuit diagram. Thus, we made a pulsed
power supply using a 3-step pulse forming network (PFN) [4, 5, 6]. PFN
consists of a series of energy storage capacitors for high-voltage and inductors which are interconnected as a ladder network. The circuit diagram ofthe supply is shown in Fig. 5.4. The edgewise coil in the vacuum chamber
has a self-inductance of 7 µH.
The oil capacitors made by General Atomics corporation (Model 39745)
are rapidly charged by means of a high-voltage direct-current power supply
via a high-voltage switch; the charging voltage is VC . The switch is GH2(withstand voltage; 10 kV and allowable current; 25 A) made by Gigavac
corporation. Another assemblies for storing electrical energy consists of acharging resistor (100 kΩ), and a potential divider (1 GΩ and 1 MΩ).
The section for releasing high pulsed current consists of three capacitors(withstand voltage; 5 kV, capacitance C; 32 µF), two self-made inductors(inductance L; 2 µH), a metal clad wire wound resistor, a fast high-voltage
solid-state diode and a thyristor switch as shown in Fig. 5.5. The inductanceof the inductors is determined by the measurement with probe and calculation [7]. The diode and the thyristor switch are FDA-100-300 (withstandvoltage; 16 kV and peak current; 3 kA) and HTS 120-1600-SCR (withstandvoltage; 12 kV and peak current; 16 kA) made by Behlke corporation, respec-

CHAPTER 5. COMPARISON OF ION EXPANSION VELOCITY

84

tively. They have much faster switching capability because of their unipolar
conduction. The functional drawback of the thyristor switch, like a diode,
only conducts in one direction.
To prevent an oscillating current, a matching impedance ZPFN is needed.
The impedance is the following equation.
√
ZPFN =LC(5.1)

The impedance is determined to be 250 mΩ from Eq. (5.1). The resistor
is IRH300 (power consumption; 300 W) made by PCN corporation. PFN accumulates electrical energy over a comparatively long time, and then releases
the stored energy in the form of a relatively square pulse of comparatively
short duration. The duration τPFN is expressed as follows, where n is thenumber of step of PFN. Here, n is three.

√τPFN = 2n LC(5.2)

The peak current IPFN of PFN is also expressed as follows. The supply
can theoretically generate a maximum current of 10 kA with 48 µs when the
charging voltage is 5 kV.

IPFN =

VC
2ZPFN

(5.3)

CHAPTER 5. COMPARISON OF ION EXPANSION VELOCITY

85

SCR
SW 100 kΩ
Power
1 GΩ

2 µΗ

2 µH

Chamber
250 mΩ

32 µF

32 µF

32 µF

Coil
7 µΗ

1 MΩ

Figure 5.4: Circuit diagram of pulsed power supply with a 3-step PFN.

Figure 5.5: Pulsed power supply made by myself.

CHAPTER 5. COMPARISON OF ION EXPANSION VELOCITY

86

Rogowski coil
Rogowski coil [8] is an electrical device for measuring high-frequency
pulsed current as shown in Fig. 5.6, where N is the number of turns and lis the circumference of the coil. The coil is a helical coil. The both terminals of lead wire are at the same end of the coil. This allows the coil to be
open-ended and flexible, and be installed and wrapped easily around a live
conductor of the pulsed power supply without disturbing the conductor. It
is highly linear for the current because the coil has no iron core to saturate
and a low inductance. Equation (5.4) indicates the voltage produced by arogowski coil, where µ0 is the permeability of vacuum, A is the area of oneof the poloidal loops, and dIPFN /dt is the derivation of the circuit current.
This formula assumes that the turns are evenly spaced.

Ein = −dIPFN
µ0 AN dIPFN
= Kr
l
dt
dt

(5.4)

The output of the rogowski coil is connected to the integral and noninverting amplifier circuits. Figure 5.7 indicates equivalent circuit diagramof these circuit. The output of the external circuit is proportional to thecurrent as shown in Eq. (5.5), where Af is 1 + 100kΩ/4.7kΩ.

Eout =Kr Af
IPFN
1kΩ × 1µF

(5.5)

Since 0.001Kr Af in Eq. (5.5) calibrated and determined, the current of
this pulsed power supply is measured from the output of this rogowski coil.

CHAPTER 5. COMPARISON OF ION EXPANSION VELOCITY

87

Figure 5.6: Rogowski coil for measuring a pulsed current.

1 kΩ
Ein

1 µF

+
100 kΩ

Eout

4.7 kΩ

Figure 5.7: Equivalent circuit diagram of integral and non-inverting amplifier
circuits.

CHAPTER 5. COMPARISON OF ION EXPANSION VELOCITY

88

Figure 5.8 shows the experimental and simulated output currents of this
supply when measuring a rogowski coil and a potential drop across a 0.002mΩ current shunt resistor for a charging voltage of 230 V. The noises at
the time of 0 s, -31 and 59 µs are caused by a laser injection and a flip of
the thyristor switch. The monitored currents agree well with the simulated
value. The peak current is 450 A at the time of 10 µs. The current starts
to flow 31 µs before laser ignition induced by a digital delay/pulse generator
for generating the strongest peak current at the laser ignition.

CHAPTER 5. COMPARISON OF ION EXPANSION VELOCITY

89

500

400

Current, A

300

200

100
Shunt value
Rogowski value

0

-100
-40

Simulated value

-30

-20

-10

0

10

Time,

20

30

40

50

60

s

Figure 5.8: Experimental and theoretical output currents of pulsed power
supply at the peak current of 450 A. Current starts flowing 31 µs before laser
ignition.

CHAPTER 5. COMPARISON OF ION EXPANSION VELOCITY

5.3

90

Numerical simulation with and without
considering creation process of laser-produced
plasma5.3.1

Calculation model of Star-1D

To compute this experiment using a CH spherical target with a diameterof 500 µm, the target is divided into six parts. Each part has a solid angle of
0.636π steradian, as shown in Fig. 5.9. In the experiment, a single laser beam
is focused onto the target. In contrast, a single laser beam is assumed to befocused onto each part in Star-1D. The surface area of the part in Star-1D
is equal to the laser beam spot size (i.e., the full width at half maximum of
the laser beam). Each laser beam has energy of 0.636 J; i.e., the total energyof the laser beams focused onto the target is about 4 J. The laser beam
is generated by a 1064 nm Nd:YAG laser with a Gaussian intensity profile.
Figure 5.10 shows the Gaussian time-variable power profile of the laser beam.
The rise and fall times of the peak power are 10 and 20 ns, respectively. The
peak power of the laser beam is 2.51 × 108 W at the calculation time of 25
ns. The distance between the laser focal point and the target is set to 0.5 cm.
The number of radial meshes with widths in the range 0.01 to 1 µm is 500.
They move with the hydrodynamic evolution of the plasma heated by thelaser beam. For a CH material, the atomic mass averaged over carbon andhydrogen and the mass density are assumed to be 6.5 amu and 1.06 g/cm3 ,respectively. The calculation start and finish times are set to zero and 75 ns,
respectively. The calculation time step ∆t is initially set to 0.1 ps and it issubsequently varied according to the CFL condition [9].

CHAPTER 5. COMPARISON OF ION EXPANSION VELOCITY

91

Spherical target
500 µmNd:YAG laser
0.5 cmOrigin 0.636π sr

Cone region
Figure 5.9: Lagrangian calculation model for ray-tracing.

3.0x102.5x10Power, W

2.0x1082.51 x 1081.0x105.0x10W

8

10 ns
1.5x10

8

20 ns

8

8

7

0.0
0

5

10 15 20 25 30 35 40 45 50 55 60 65 70 75
Time, ns

Figure 5.10: Gaussian time-variable power profile of a single laser beam.

CHAPTER 5. COMPARISON OF ION EXPANSION VELOCITY

5.3.2

92

Results and discussion of Star-1D

Figure 5.11 shows the results of Star-1D at the calculation time of 75 ns.
The horizontal axis represents the radius of the laser-produced plasma fromthe center of plasma. The ablation and expansion surfaces are located atradii of 0.25 and 9 mm. The vertical axes represent the 1D ion density ni ,the ion expansion velocity v, and the charge state of ion Z in the plasma.
The material between the center and a radius of 0.25 mm is not considered
to be ionized because the ion expansion velocity in this region is zero.
Figure 5.12 indicates the time-variation of plasma kinetic, internal andradiation energies. The absorption rate is 93% of the incident laser energy(4.0 J). The plasma has kinetic and inertial energies at 75 ns of 3.1 J (78%of incident laser energy) and 0.26 J (7% of the incident laser energy), respectively. The remaining absorption energy is the heat loss from the expansion
surface. The ablation plasma mass is 2.1×10−8 kg (30%) for a target mass of
6.9 × 10−8 kg. The ablation plasma expanding into vacuum expands nearly
isothermally. Consequently, the ion expansion velocity increases to 175 km/s
and the ion density decreases exponentially, as shown in Fig. 5.11. The
charge state of ion near a radius of 5 mm is rapidly altered by shock waves.
The results of Star-1D are used as the initial plasma conditions, which are
inputted into the 3D hybrid PIC code. The conditions are 3D distributionsof ion super particle for the particle position, expansion velocity, charge stateand mass.

CHAPTER 5. COMPARISON OF ION EXPANSION VELOCITY

1

0

10

24

10

23

10

22

10

21

10

20

10

19

10

18

10

17

10

16

40

10

15

20

10

14

180

Z

160
140
120

v

100
80

n

0

1

2

3

60
i

4

5

6

7

8

9

Ion expansion velocity, km/s

-3

Ion density, cm

Charge state of ion

2

93

0

Laser-produced plasma radius, mm

Figure 5.11: 1D distributions of plasma for ion density, ion expansion velocity,and charge state as a function of plasma radius at 75 ns.
4Plasma kinetic energy
Plasma internal energy

Energy, J

3

Radiation energy

3.06

2

1
0.37

0

0.26

0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75
Time, ns

Figure 5.12: Time-variation of plasma kinetic, internal and radiation energies.

CHAPTER 5. COMPARISON OF ION EXPANSION VELOCITY

5.3.3

94

Calculation model of 3D hybrid PIC code

The initial laser-produced plasma is located at the origin of the model asshown in Fig. 5.1. The model is an 110∆x × 110∆x × 110∆x cube, where∆x=3 mm is the mesh width. Figure 5.13 shows the magnetic field strength
along the central axis of z at x = y = 0 and in the z − x plane at y = 0
produced by the edgewise coil at the peak current of 450 A. The magnetic
field strength at the origin is calculated to be 0.0321 T. The calculation start
time is set to zero. At this time, numerical simulations are changed from
Star-1D to the 3D hybrid PIC code. The number of ion super particles is
1,000,000. The calculation time is 100,000∆t (5 µs), where ∆t=50 ps (ioncyclotron frequency ωci =1/2047∆t) is the time step. It is according to CFL
condition. Table 5.1 lists the characteristic plasma conditions for the original
and modified cases.

CHAPTER 5. COMPARISON OF ION EXPANSION VELOCITY

95

Magnetic field strength, T

0.240
0.227

0.12

0.215

0.10

0.202

0.08

0.189

0.06

0.177

0.04

0.0321 T

0.164
0.152

0.02
0.00
-100
75

0.139
-75

-50

-25

0

25

50

0.126
0.114

Magnetic field strength on the z-x plane, T

0.101
0.0884x, mm

50

0.07580.06320.050525

0.03790.02530.0126

0
-100

-75

-50

-25

0

25

50

0

z, mm

Figure 5.13: Calculated magnetic field strength produced by the edgewise
coil at the peak current of 450 A. The magnetic field strength at the origin
is 0.0321 T.

Table 5.1: Plasma conditions for original and modified cases.
Original case

Modified case

Plasma radius, mm

99Plasma mass, kg

6.9 × 10−8

2.1 × 10−8

Atomic mass, amu

6.5

6.5

<Atomic number>

1.25

0.5

Ion density, cm−3

2.1 × 1018

non-uniform

Kinetic energy, J

1.2

3.1

Ion expansion velocity Maxwell distribution

non-uniform

CHAPTER 5. COMPARISON OF ION EXPANSION VELOCITY

5.4

96

Results and discussion

The ion expansion velocity was measured by the charge collector with a
Faraday cup in the experiment and calculated by the modified and original
numerical simulations. The values were derived by dividing the distance
between the target and the probe (155 mm) by the arrival time of ion beam
current.
In the experiment, the measured ion expansion velocity at B=0.0321 T
is 199 km/s. The uncertainty in the time range at the oscilloscope is in
the order of nanoseconds, so that the uncertainty in the measured velocity
is negligible. In the original numerical simulation without considering thecreation process of a laser-produced plasma, the velocity is calculated to be
15 km/s. The velocity obtained using the modified numerical simulation
with considering the creation process is 194 km/s. Thus, the experimentallymeasured velocity and that calculated by the modified numerical simulation
at B=0.0321 T are almost the same; that is, the initial plasma velocity forthe modified numerical simulations is adequate for the experiment.

5.5

Summary

In this Chapter, we described ion expansion velocity of a laser-produced
plasma in a magnetic field produced by a edgewise coil.
In the experiments at EUV database laser facility of ILE in the collaboration with Osaka University, to generate a magnetic field, a pulsed power
supply and an edgewise coil were made. The ion expansion velocity was measured by using a charge collector with a Faraday cup. To confirm the validityof the modified numerical simulation, the results obtained using the modified
numerical simulation were compared with those of experiments and the original numerical simulation. The modified numerical simulation compute high

CHAPTER 5. COMPARISON OF ION EXPANSION VELOCITY

97

energy laser injection onto a spherical target and the ablation plasma dynamics in a magnetic field. The measured ion expansion velocity at B=0.0321
T was 199 km/s. In the numerical simulations, the calculated velocity forthe modified and original numerical simulations were 194 km/s and 15 km/s,
respectively. Comparison of the ion expansion velocities revealed that theexperiment and the modified numerical simulation at B=0.0321 T give almost the same ion expansion velocities. These results demonstrate that themodified numerical simulation can describe the ablation plasma dynamics
in the magnetic field but the original numerical simulation does not. Thisconfirms the validity of the modified numerical simulation.
The achievements in this Chapter are as follows.
Experiments• In order to measure the ion beam current of a laser-produced plasma ina pulsed magnetic field, charge collector with a Faraday cup and pulsed
power supply using a three-step PFN were made.
• By measuring the ion beam current was measured, the ion expansion
velocity was estimated.
Numerical simulations• The measured ion expansion velocity was compared with the calculatedone.
• The experimentally measured velocity and that calculated by the modified numerical simulation are almost the same; that is, the initialplasma velocity for the modified numerical simulation is adequate forthe experiment.

Bibliography[1] Pearlman, J. S., “Faraday Cups for Laser Plasmas,” Rev. Sci. Instrum.,Vol. 48, Issue 8, 1977, pp. 1064-1067. doi: 10.1063/1.1135184[2] Sakawa, Y., Watanabe, D., Shibahara, D., Sugiyama K., and Tanabe, T., “Pulsed-laser Ablation of Co-deposits on JT-60 Graphite
Tile,” J. Nuclear Materials, Vol. 367370, 2007, pp. 1516-1521. doi:
10.1016/j.jnucmat.2007.04.014[3] Ikeda, S., Nakajima, M., and Horioka, K., “Control of Laser Ablation
Plasma with Longitudinal Magnetic Field,” Plasma Fusion Res., Vol. 7,
2012, 1201015. doi: 10.1585/pfr.7.1201015[4] Nojiri, H., “パルス強磁場用小型コンデンサ電源の自作,” Solid State
Physics, Vol. 37, Issue 7, 2002, pp. 465-472. [In Japanese][5] Akiyama, H., 高電圧パルスパワー工学, EE Text Series, Ohmsha, Tokyo,
2003, Chaps. 10 and 11. [In Japanese][6] Minami, T., Funaki, I., Yamakawa, H., and Nakayama, Y., “Experimental Study of Magnetic Sails,” ISAS Research Note 790, 2005. [InJapanese][7] Hayashi, I., “Formulas for Use in Calculating Electromagnetic Forceson the Cylindrical Air Core Coils and its Applications for the Power
Reactors,” J. IEEJ, Vol. 81, Issue 9, 1961, pp. 60-65. [In Japanese]98BIBLIOGRAPHY

99[8] 川崎昌二, “ロゴスキーコイルによる交流大電流測定,” 九州大学応用力学研究所技術レポート, Vol. 1, 2000, pp. 14-23. [In Japanese][9] Ohsawa, Y., “Relativistic Electromagnetic Particle Simulation,” J.
Plasma Fusion Res., Vol. 74, Issue 7, 1998, pp. 746-752. [In Japanese]Chapter 6Eﬀect of magnetic fieldstrength on impulse
6.1

Introduction

In the previous Chapters, it was brought up the diﬀerences between themodified and original numerical simulations, which are with and without
considering the creation process of a laser-produced plasma for the impulse
and ion expansion velocity under the same conditions of the magnetic fieldstrength and geometry. It is considered that the initial plasma condition inputted into the 3D hybrid PIC code is important. In this Chapter, the detailsof the initial plasma condition aﬀecting the impulse is analyzed by using themodified and original numerical simulations. Finally, the interaction between
initial plasma condition and magnetic field strength is investigated.

100

CHAPTER 6. EFFECT OF MAGNETIC FIELD

6.2

101

Dependence on ion expansion velocity distribution in initial laser-produced plasma

6.2.1

Calculation model of 3D hybrid PIC code

The calculation model of 3D hybrid PIC code is shown in Fig. 6.1. The
initial plasma is located at the center of model. The model is a 100∆x×
100∆x× 100∆x cube, where ∆x=1 mm is the mesh width. The edgewise
coil that generates the pulsed magnetic field is formed by winding five turnsof a rectangular copper wire (width: 5 mm; thickness perpendicular to thecoil axis: 1 mm; external diameter: 18 mm; internal diameter: 8 mm; length:
5 mm). The coil is set at a distance of 6.5 mm from the center. Figure 6.2
shows the magnetic field strength along the z-axis at x = y = 0 and in thez − x plane at y = 0 produced by the coil at the peak current of 5.2 kA.
The magnetic field strength at the center is 0.51 T. The calculation start
time is set to 75 ns. The calculation time is 100,000∆t, where ∆t=50 ps(ion cyclotron frequency ωci =1/2047∆t) is the time step. The numericalsimulation finishes at 5.075 µs. The number of particles is 1,000,000.
The plasma condition in Chap. 5 are used as the initial conditions forthe calculations. The conditions are 3D distributions of ion super particle forthe particle position, expansion velocity, charge state and mass. Table 6.1
lists the characteristic plasma conditions for the original and modified cases.

CHAPTER 6. EFFECT OF MAGNETIC FIELD

x

Coil
18 mm

102

6.5 mmPlasma

8 mm

z

y
5 mm
100 mm1 mm

Figure 6.1: Cubic calculation model of 3D hybrid PIC code.

Magnetic field strength, T

2.85
2.70

3.0

2.55

2.5

2.40

2.0

2.25

1.5

2.10

1.0

1.95
0.51 T

0.5
0.0
-30
20

1.80
1.65

-20

-10

0

10

1.50
1.35

Magnetic field strength on z-x plane, T

1.20
1.05

x, mm

0.900
0.750

10

0.600
0.450
0.300
0.150

0
-30

-20

-10

0

10

0

z, mm

Figure 6.2: Calculated magnetic field strength produced by the edgewise coilwhen the peak current is 5.2 kA. The magnetic field strength at the center
is 0.51 T.

CHAPTER 6. EFFECT OF MAGNETIC FIELD

103

Table 6.1: Plasma conditions for original and modified cases.
Original case

Modified case

Plasma radius, mm

99Plasma mass, kg

6.9 × 10−8

2.1 × 10−8

Atomic mass, amu

6.5

6.5

<Atomic number>

1.25

0.5

Ion density, cm−3

2.1 × 1018

non-uniform

Kinetic energy, J

1.2

3.1

Ion expansion velocity Maxwell distribution

non-uniform

CHAPTER 6. EFFECT OF MAGNETIC FIELD

6.2.2

104

Results and discussion

Diamagnetic cavity
Figures 6.3 and 6.4 show the distribution of the magnetic field strength in
the z − x plane at y = 0 at calculation times of 325, 575, 825 ns, 1.075, 1.325,
1.575, 3.075 and 3.575 µs for the original case. Comparison with the initial
magnetic field shown in Fig. 6.2 reveals that the magnetic field disappears
inside the area indicated by the red circle in Fig. 6.3(a). The diamagnetic
current in a laser-produced plasma forms a diamagnetic cavity. Figures 6.3
and 6.4 show that the magnetic field returns to the initial magnetic field up
to 3.575 µs and the cavity moves in the +z direction for pushing plasmaback.
Figure 6.5 shows the distribution of the magnetic field strength in thez − x plane at y = 0 at calculation times of 1.325, 2.575, 3.825 and 5.075 µs
for the modified case. The distribution does not change from 75 ns to 1.325µs. Comparison with the initial magnetic field shown in Fig. 6.2 reveals thatthe magnetic field disappears inside the area indicated by the red circle inFig. 6.5(a). The diamagnetic current in a laser-produced plasma forms adiamagnetic cavity. Figure 6.5 shows that the magnetic field returns to theinitial magnetic field up to 5.075 µs and the cavity moves in the +z direction
for pushing the plasma back.

CHAPTER 6. EFFECT OF MAGNETIC FIELD

105(a) 325 ns & 575 ns

20

2.85

Magnetic field strength on z-x plane at thecalculation time of 825 ns, T

2.70
2.55

x, mm

2.40
2.25

10

2.10
1.95
1.80
0
-30
20

1.65
-20

-10

z, mm

0

Magnetic field strength on z-x plane at thecalculation time of 1.075

10

1.50
1.35

s, T

1.20
1.05

x, mm

0.900
0.750

10

0.600
0.450
0.300
0.150

0
-30

-20

-10
z, mm

0

10

0.00

(b) 825 ns & 1.075 µs

Figure 6.3: Magnetic field strength distribution for the original case in z − xplane at y = 0 for calculation times of 325, 575, 825 ns, and 1.075 µs. A
diamagnetic cavity is formed.

CHAPTER 6. EFFECT OF MAGNETIC FIELD

20

Magnetic field strength on z-x plane at thecalculation time of 1.325

106

s, T

2.85
2.70
2.55

x, mm

2.40
2.25

10

2.10
1.95
1.80
0
-30
20

1.65
-20

-10
z, mm

0

Magnetic field strength on z-x plane at thecalculation time of 1.575

10

1.50
1.35

s, T

1.20
1.05

x, mm

0.900
0.750

10

0.600
0.450
0.300
0.150

0
-30

-20

-10
z, mm

0

10

0.00

(a) 1.325 µs & 1.575 µs

20Magnetic field strength on z-x plane at thecalculation time of 3.075

s, T

2.85
2.70
2.55

x, mm

2.40
2.25

10

2.10
1.95
1.80
0
-30
20

1.65
-20

-10

z, mm

0

Magnetic field strength on z-x plane at thecalculation time of 3.575

10

1.50
1.35

s, T

1.20
1.05

x, mm

0.900
0.750

10

0.600
0.450
0.300
0.150

0
-30

-20

-10
z, mm

0

10

0.00

(b) 3.075 µs & 3.575 µs

Figure 6.4: Magnetic field strength distribution for the original case in z − xplane at y = 0 for calculation times of 1.325, 1.575, 3.075, and 3.575 µs. A
diamagnetic cavity is formed.

CHAPTER 6. EFFECT OF MAGNETIC FIELD

107(a) 1.325 µs & 2.575 µs

20Magnetic field strength on z-x plane at thecalculation time of 3.825

s, T

2.85
2.70
2.55

x, mm

2.40
2.25

10

2.10
1.95
1.80
0
-30
20

1.65
-20

-10

z, mm

0

Magnetic field strength on z-x plane at thecalculation time of 5.075

10

1.50
1.35

s, T

1.20
1.05

x, mm

0.900
0.750

10

0.600
0.450
0.300
0.150

0
-30

-20

-10
z, mm

0

10

0.00

(b) 3.825 ns & 5.075 µs

Figure 6.5: Magnetic field strength distribution for the modified case in z − xplane at y = 0 for calculation times of 1.325, 2.575, 3.825, and 5.075 µs. A
diamagnetic cavity is formed.

CHAPTER 6. EFFECT OF MAGNETIC FIELD

108

Thrust energy and impulse
The eﬃciencies in terms of the kinetic energy ηkin and momentum ηIbit
are calculated using Eqs. (6.1) and (6.2), where v±z is the velocity of asuper particle which moves in the ±z direction. Ibit stands for impulse bit(impulse per pulse). The sums are performed over all super particles. Figure
6.6 shows the time evolutions of the thrust energy and impulse for the original
and modified cases.
∑1
2
m
v
−
m v2
i
+z
2
2 i −z
∑1
m v2
2 i 0

(6.1)

∑
mi vz − mi vz0∑
mi |vz0 |

(6.2)

∑1
ηkin =∑
ηIbit =For the original case, the thrust energy increases gradually up to 3 µs,
after which it saturates at 0.36 J; the final eﬃciency ηkin is 31%. In a similar
fashion, the impulse increases and saturates at 151 µNs; the final eﬃciency
ηIbit is 41%. For the modified case, the thrust energy increases rapidly up
to 200 ns, after which it saturates at 0.37 J (12%). In a similar fashion,
the impulse increases and saturates at 31.6 µNs (41%). Figure 6.6 indicatesthat the direction of movement of many super particles is changed in the +z
direction by the magnetic field of the coil, as shown in Figs. 6.3, 6.4 and6.5. The time-variable magnetic field, thrust energy, and impulse for the twocases exhibit a time diﬀerence.

CHAPTER 6. EFFECT OF MAGNETIC FIELD

109

0.5

0.5

E =0.37 J &0.4
Thrust energy, J

kin

kin

0.4

=12%

0.3

0.3

Original case
Modified case

0.2

0.2

E =0.36 J &
kin

0.1

kin

=31%

0.1

0.0

0.0
0

1

2

3
Time,

4

5

s

(a) Thrust energy

200

200

180

180
I

Impulse bit,

Ns

160

bit

=151

Ns &

Ibit

=41%

160

140

140

120

120

Original case

100

100

Modified case

80

80

60

60
I

40

bit

=31.6

Ns &

Ibit

=41%

40

20

20
0

1

2

3
Time,

4

5

s

(b) Impulse

Figure 6.6: Time evolutions of thrust energy and impulse. For the originalcase, thrust energy and impulse saturate at 0.36 J (31%) and 151 µNs (41%),
respectively. For the modified case, thrust energy and impulse saturate at
0.37 J (12%) and 31.6 µNs (41%), respectively.

CHAPTER 6. EFFECT OF MAGNETIC FIELD

110

Discussion
When the Larmor radius is shorter than the coil radius, it is consideredthat the charged particles move in the +z direction eﬀectively. The thresholdof the Larmor radius is set as 5 mm which is a bit larger than the internal
radius of the edgewise coil. The region of slow expansion velocity in therange 0-20 km/s is focused using the threshold. In this region, the charge
state and ion density are almost same for the two cases.
For the original case, most initial plasma is in the region of slow expansion
velocity. For the modified case, the initial plasma can be divided into twogroups: group A that has a slow expansion velocity (kinetic energy: 0.027
J), and group B that has a fast expansion velocity in the range 20-175 km/s
(kinetic energy: 3.0 J). Figure 6.7 shows the thrust energy and impulse of
these groups. For the group B as shown in Fig. 6.7(b), the thrust energy andimpulse increase rapidly up to 200 ns, after which they saturate at 0.34 J(11%) and 7.5 µNs (14%), respectively. The group B forms the diamagnetic
cavity but it does not contribute the impulse. For group A as shown in
Fig. 6.7(a), the thrust energy and impulse increase slowly up to 1.5 µs, after
which they saturate at 0.0337 J (125%) and 24.1 µNs (102%), respectively.
Thus, the group A contributes the impulse eﬃciently. The reason for thehigh eﬃciency is that the group A moves in the diamagnetic cavity formedby the group B and pushed back by the diamagnetic cavity.

CHAPTER 6. EFFECT OF MAGNETIC FIELD

111

0.05

0.05
Thrust energy, J

0.04
0.03

E =0.0337 J &

kin

kin

0.04
0.03

=125%

0.02

0.02

0.01

0.01

0.00

0.00
30

30
Impulse bit,

Ns

20

20
I

bit

=24.1

Ns &

Ibit

=102%

10

10

0
0.0

0
0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

5.0

s

Time,

(a) Slow velocity group

0.5

0.5
Thrust energy, J

0.4

0.4

0.3

0.3

E =0.34 J &

0.2

kin

kin

=11%

0.2

0.1

0.1

0.0

0.0

10

10
Impulse bit,

8

Ns

8

6

6

4

I

bit

2

=7.5

Ns &

0

Ibit

=14%

4
2
0

0

50

100

150

200

Time, ns(b) Fast velocity group

Figure 6.7: Time evolutions of thrust energy and impulse for (a) slow velocity
group and (b) fast velocity group. For group A, thrust energy and impulse
saturate at 0.0337 J (125%) and 24.1 µNs (102%). For group B, thrust energy
and impulse saturate at 0.34 J (11%) and 7.5 µNs (14%), respectively.

CHAPTER 6. EFFECT OF MAGNETIC FIELD

6.3

112

Thrust performances for various magnetic
field strengths

6.3.1

Numerical simulation conditions

In the experiment of Chap. 5, an impulse due to the interaction between
a diamagnetic current in a laser-produced plasma and a magnetic field generated by an edgewise coil has been measured by using a torsional thrust stand.
However, the level of the impulse has been lower than the detection limit ofa torsional thrust stand. This is considered that the magnetic field strengthof the experiment (B=0.0321 T) has been too weak. When the magnetic
field strength is suﬃciently high, the impulse should increase. Calculations
were conducted at the origin of stronger magnetic field, which strengths are
B=0.0321, 0.961, 1.92, and 2.88 T.

6.3.2

Results and discussion

Figures 6.8-6.11 show the ion super particle distributions at a calculation
time of 5.075 µs in the z − x plane at y = 0 for various magnetic fields
calculated by the two numerical simulations. They show the eﬀect of themagnetic field strength in the magnetic thrust chamber system for LFR. The
squares in the figures indicate the edgewise coil.
For the modified case, most particles of B=0.0321 T are expanding isotropically, while the other particles are excluded by colliding the edgewise coil asshown in Fig. 6.8(a). This indicates that the interaction between the particles and the magnetic field is too weak to generate the impulse. Then, Figs.
6.9(a), 6.10(a) and 6.11(a) show that the directions of particles movement
are limited along z-axis with increasing magnetic field strength. The particles move in the exhaust direction with anisotropic expansion. Thus, theimpulse is predicted to be improved by increasing the magnetic field strength.

CHAPTER 6. EFFECT OF MAGNETIC FIELD

113

The reason for this is considered that the particles of fast expansion velocity
are interacting with the magnetic field and the Larmor radius is decreasing
with increasing the strength. This indicates that the interaction between theparticles and the magnetic field is stronger with increasing the strength.
For the original case, Figs. 6.8(b), 6.9(b), 6.10(b) and 6.11(b) show thatthe particles are guided by the magnetic field through the edgewise coil. The
particle distribution at B=0.961 T shown in Fig. 6.9(b) is similar to theresults for B=1.92 T shown in Fig. 6.10(b) and for B=2.88 T shown in Fig.
6.11(b). The directions of particles movement are limited enough along zaxis. This suggests that the impulse will saturate at magnetic field strengths
greater than B=0.961 T.
Figure 6.12 shows the impulses calculated by the modified and original
numerical simulations. As the ion super particle distributions in Figs. 6.86.11 show, the impulse in the modified numerical simulation increases with increasing magnetic field strength from B=0.961 T, whereas the impulse in theoriginal numerical simulation saturates at high magnetic field strengths. The
reason for this proportion relation between the magnetic field strength andthe impulse is that the diamagnetic cavity which increases with increasing themagnetic field strength interacts more with the particles of slow expansion
velocity. Therefore, it is revealed that the impulse increases with increasing
magnetic field strength and saturates at higher magnetic field strength.

CHAPTER 6. EFFECT OF MAGNETIC FIELD(a) Modified case.

114(b) Original case.

Figure 6.8: Ion super particle distributions at B=0.0321 T calculated bymodified and original numerical simulations.

(a) Modified case.

(b) Original case.

Figure 6.9: Ion super particle distributions at B=0.961 T calculated by modified and original numerical simulations.

CHAPTER 6. EFFECT OF MAGNETIC FIELD(a) Modified case.

115(b) Original case.

Figure 6.10: Ion super particle distributions at B=1.92 T calculated by modified and original numerical simulations.

(a) Modified case.

(b) Original case.

Figure 6.11: Ion super particle distributions at B=2.88 T calculated by modified and original numerical simulations.

CHAPTER 6. EFFECT OF MAGNETIC FIELD

116

350

300

Modified numerical simulation
Original numerical simulation

Impulse bit,

Ns

250

200

150

100

50

0
0.0

0.5

1.0

1.5

2.0

2.5

3.0

Magnetic field strength, T

Figure 6.12: Impulse for various magnetic field strengths calculated by modified and original numerical simulations.

CHAPTER 6. EFFECT OF MAGNETIC FIELD

6.4

117

Summary

In this Chapter, we described the dependences of the interaction between
initial laser-produced plasma property and magnetic field strength by comparing the two numerical simulations. The following results were obtained.
• The particles of slow expansion velocity contributes the impulse eﬃciently. The reason for the high eﬃciency is that the diamagnetic cavity
produced by the particles of fast expansion velocity interact with theparticles of slow expansion velocity.
• Calculations of the impulse reveal that the impulse increases with increasing magnetic field strength and saturates at higher magnetic fieldstrength.

Chapter 7
Summary
In this study, the followings are conducted to study the magnetic thrustchamber system for LFR.
1. The direct measurement of the impulse due to the interaction between
a diamagnetic current and a magnetic field has not been demonstrated
experimentally yet. To produce the impulse by using laser-produced
plasma and magnet, we attempt to conduct experiments in the actual
research laboratory instead of ICF plasma and SCM. The experiments
are conducted for various laser pulse energies at Extreme Ultra-Violet(EUV) database and Gekko XII (GXII) laser facilities of Institute ofLaser Engineering (ILE) in the collaboration with Osaka University.
2. The creation process of a laser-produced plasma has not suﬃciently
been considered in the numerical simulation yet. This process is needed
to compare numerical simulations with experiments in the actual research laboratory. By combining a 3D hybrid PIC code with a onedimensional (1D) Lagrangian radiation hydrodynamic code which cancompute the hydrodynamic evolution of radiating plasma heated bylaser pulses, we attempt to compute high energy laser injection onto aspherical target and the ablation plasma dynamics in a magnetic field.
118

CHAPTER 7. SUMMARY

119

3. To find the scaling law of the magnetic thrust chamber system, it isnecessary to confirm the validity of the numerical simulation. In order to validate the modified numerical simulation with considering thecreation process, the experimental results are compared with these of
the numerical simulations which are with and without considering thecreation process of a laser-produced plasma.
As a result, the achievements in this study are obtained as follows.
• By using cylindrical neodymium permanent magnet and non-magnetizediron cylinder, impulse due to the interaction between a diamagnetic
current in a laser-produced plasma and a magnetic field was generated
and the magnetic thrust chamber system for LFR was demonstrated
experimentally to work as theory tells.
• The modified numerical simulation with considering the creation process of a laser-produced plasma as the initial plasma conditions of the
3D hybrid PIC code was performed. The modified numerical simulation
can describe the ablation plasma dynamics in the magnetic field but
the original numerical simulation does not. This confirms the validityof the modified numerical simulation.
• The experimentally measured velocity and that calculated by the modified numerical simulation are almost the same; that is, the initialplasma velocity for the modified numerical simulation is adequate forthe experiment.
• In order to produce large impulse, longer pulse durations, shorter wavelengths and higher energies are better.
• Calculations of the impulse reveal that the impulse increases with increasing magnetic field strength and saturates at higher magnetic fieldstrength.

CHAPTER 7. SUMMARY

120• The particles of slow expansion velocity contributes the impulse eﬃciently. The reason for the high eﬃciency is that the diamagnetic cavity
produced by the particles of fast expansion velocity interact with theparticles of slow expansion velocity.
Then, the magnetic thrust chamber system for LFR is demonstrated experimentally to work as theory tells. In addition to that, the modified numerical simulation with considering the creation process of a laser-produced
plasma as the initial plasma conditions of the 3D hybrid PIC code can determine the ablation plasma dynamics in the magnetic field more accurately
than the original numerical simulation. Finally, the dependences of the impulse on the laser parameters and initial conditions are obtained in order toimprove the thrust performance of the magnetic thrust chamber system for
LFR.
The suggestions for the further work are as follows.
• Under the guidances obtained in this study for improving the thrust
performance of the magnetic thrust chamber system for LFR, we areplanning to conduct experiments and numerical simulation of the system. In addition to that, it is also necessary to consider how to create
the optimal laser-produced plasma for the system in the experiments.
• The modified numerical simulation can describe the ablation plasma
dynamics in the magnetic field but the original numerical simulation
does not. This confirms the validity of the modified numerical simulation. However, there is the discrepancy between measured and calculated values. The reason for the discrepancy is considered to be due tothe calculation assumptions in Star-1D and the 3D hybrid PIC code.
Star-1D assumes that a laser-produced plasma expands into the vacuum isotropically. However, in the experiments, the plasma created bya single laser beam is ablating in one direction. This is considered that

CHAPTER 7. SUMMARY

121

the ablation plasma on the opposite surface of laser injection is not
calculated well in the modified numerical simulation. This diﬀerenceof plasma ablation may aﬀect initial plasma property. Then, we mustuse a 2D radiation hydrodynamic code which can compute the axial
symmetry hydrodynamic evolution of radiating plasma heated by laser
pulses.
• Because it is impossible to conduct experiment using ICF plasma and
SCM in several decades, to investigate the magnetic thrust chambersystem of LFR is proceeding by computing high energy laser injection
onto a ICF target and the ablation plasma dynamics in a magnetic field
produced by SCM.

Acknowledgments
I would like to express the deepest appreciation to my doctoral advisors,
Professor Hideki Nakashima, Professor Kazuo Nakamura, Professor Kiichiro
Uchino, and Associate Professor Naoji Yamamoto.
I would like to express the appreciation to my advisers, Dr. Atsushi
Sunahara of Institute for Laser Technology, Dr. Tomoyuki Johzaki of Hiroshima University, Dr. Yoshihiro Kajimura of Akashi National College ofTechnology, Dr. Shinsuke Fujioka of Osaka University, and Yoshitaka Mori of
The Graduate School for the Creation of New Photonics Industries. I would
like to appreciate Megumi Oogami, Kazuyo Mine, Mieko Inoue, and Takako
Mawatari who take care of my various paperwork.
I would like to thank my junior labmates whose work about my research
theme together, Masato Yasunaga, Tomoyuki Hinaga, Yasuma Hidehira, Hiroshi Tominaga and Yudai Miwa. In addition, I thank my junior labmates of
Kyushu University, Yusuke Nakamura, Kenji Hirano, Naruya Hiroike, Yutaro
Sasagawa, Yuji Toyoda, Yuta Tajiri, Hiroshi Umishio, Shouhei Araki, Kohei
Nagaoka, Shouta Maeda, Yuki Sakamoto and Keisuke Kataoka.
This work was supported by the Japan Society for the Promotion of
Science (JSPS) KAKENHI (Grant Number 21360418) and Grant-in-Aid for
JSPS Fellows (Grant Number 22-1672) and Collaboration Research with theInstitute of Laser Engineering at Osaka University since 2008.
Finally, I would like to express the deepest gratitude to my friends, myfamily and my precious one for their kind support and help.
122

Appendix A
Ampere’s law with Darwin
approximation
Here, we derive Ampere’s law with Darwin approximation. Equation(A.1) is Gauss’s law for magnetism, where B is magnetic field.

∇·B =0(A.1)

Because magnetic field is solenoidal field, magnetic field is expressed byvector potential A (i.e., ∇ · A = 0).

B =∇×A(A.2)

Equation (A.3) is Maxwell-Faraday equation, where E is electric field.

∇×E =−

∂B
∂t

(A.3)

By substituting Eq. (A.2) into Eq. (A.3), the following equation is obtained.

123

APPENDIX A. AMPERE’S LAW WITH DARWIN APPROXIMATION124

)
(
∂A
=0
∇× E+
∂t

(A.4)

Electric field is expressed as follows by using Eq. (A.4) and scalar potential ϕ (i.e., ∇ × ∇ϕ = 0).

E = −∇ϕ −

∂A
∂t

(A.5)

Equation (A.6) is Ampere’s law, where j is current density, µ0 is thepermeability of vacuum and ε0 is the permittivity of vacuum.
)(∂E∇ × B = µ0 j + ε0
∂t

(A.6)

By substituting Eq. (A.5) into Eq. (A.6), the following equation is obtained.
)(1 ∂ 2A∂ϕ∇ × B = µ0 j − ε0 ∇
− 2 2
∂t
c ∂t

(A.7)

The first and second terms on r.h.s of Eq. (A.7) have the order of (v/c)2and (v/c)3 , respectively. When the speed of light c (= 1/µ0 ε0 ) is more fasterthan plasma velocity, the second term on r.h.s. of Eq. (A.7) is neglected.
Darwin approximation means to neglect the transverse displacement current.
Ampere’s law with Darwin approximation is given as shown in Eq. (A.8).
()∂ϕ∇ × B = µ0 j − ε0 ∇
∂t

(A.8)

Current density is separeted to the density of solenoidal field jr (i.e.,∇ × jr = 0) and that of potential field jd . Equation (A.9) is obtained by
the curl of Eq. (A.8) and the relationship between curl and gradient (i.e.,∇ × ∇ϕ = 0).

APPENDIX A. AMPERE’S LAW WITH DARWIN APPROXIMATION125(∇ × (∇ × B) = ∇ × µ0 j − ∇ ×

1 ∂ϕ
∇
c2 ∂t

)

= µ0 ∇ × (jr + jd )∇ × B = µ 0 jd

(A.9)

By substituting Eq. (A.5) into Eq. (A.8), the following equation is obtained.
(
)}
{
∂A
∂
∇ × B = µ0 j + ε0
E+
∂t
∂t

(A.10)

Electric field is separeted to the part of solenoidal field Er and that ofpotential field Ed (i.e., ∇ · Ed = 0). By using Eqs. (A.9) and the divergenceof (A.10), the following equation is obtained.
}
{
∂A
∂
)
∇ · (∇ × B) = µ0 ∇ · j + ε0 (∇ · E + ∇ ·
∂t
∂t
(
)
1 ∂Er
0 = ∇ · µ0 jr + 2
c ∂t(A.11)

Equation (A.12) is Gauss’s law, where ρ is charge density. Here, it isassumed to be quasi-neutrality (ρ is constant).

∇·E =

ρ
ε0

(A.12)

By diﬀerenting partially both sides of Eq. (A.12), Eq. (A.13) is given.

∇·

∂E
∂Er
=∇·
=0
∂t
∂t

(A.13)

APPENDIX A. AMPERE’S LAW WITH DARWIN APPROXIMATION126
By substituting Eq. (A.11) into Eq. (A.13), ∇ · jr = 0 is obtained, thatis, jr = 0. Thus, j = jd . The Ampere’s law with Darwin approximation isexpressed as follows, where ji is the ion current density and je is the electron
current density.

∇ × B = µ0 (ji + je )

Appendix B
Components of matrix
equations
B.1

n+ 21
x component of vi
n+ 12

We derive x component of vi
n+ 12

for Eq. (3.18) in the Chapter 3. This isn+ 12much of the same for viy

and vizdv n
mni i = Z n e
dt

(
)
η
n
n
n
n
E + vi × B − ∇ × B
µ0

n+ 21

dvin
v
= i
dt

.

(B.1)

n− 12

− vi
∆t

n+ 12

By substituting Eq. (B.2), vi

(B.2)
n+ 1

n+ 1

n+ 1

(vix 2 , viy 2 , viz 2 ), B n (Bxn , Byn , Bzn ) andE n (Exn , Eyn , Ezn ) into Eq. (B.1), the following equation is obtained.

127

APPENDIX B. COMPONENTS OF MATRIX EQUATIONS

128

 n+ 1
n− 12
2
vix − vixmni 
 n+ 12n− 1 
viy − viy 2 
∆t 

n+ 12
n− 12
viz − viz n ∂B n 
 
∂Bz
n n
viy
Bz − vizn Byn
Exn
− ∂zy
 ηZ n e  ∂y

 


 ∂Bxn ∂Bzn 
 
n n
−
= Z n e Eyn  + Z n e vizn Bxn − vixBz − ∂x 
µ0  ∂z


 

n
n
∂B
∂B
y
n n
n n
vix
By − viy
Bx
Ezn
− ∂yx
∂x
 n+ 1  n+ 1
n+ 1
vix 2
viy 2 Bzn − viz 2 Byn n+ 12  Z n e∆t  n+ 12 n
n+ 12 n viy  −v
Bx − vix Bz 
2mni  iz



1
1
n+ 2
n+ 2 nn+ 12 nviz
vix By − viy Bx  n− 1
 n− 1 n− 1Exnviy 2 Bzn − viz 2 Byn
vix 2

 



 n− 12  Z n e∆t  n  Z n e∆t  n− 12 nn− 12 n 
= viy  +E  +
v
Bx − vix Bz 
mni  y 
2mni  iz



1
1
n− 2
n− 2 nn− 12 nnviz
Ez
vix By − viy Bx
 n ∂B n 
∂Bz
− ∂zy
∂y

ηZ n e∆t  ∂Bxn ∂Bzn 
(B.3)
−


−
∂x 
µ0 mni  ∂z
n
∂Byn
x
− ∂B
∂x
∂y
Space derivative is defined as follows.

∂Bzn= Dy Bzn∂y(B.4)

The second term of l.h.s and third term of r.h.s. in Eq. (B.3) are expressedby the following definition.
 
 
α
Bn
  Z n e∆t  x   nβ  =n By 2mi  
 
γ
Bzn

(B.5)

APPENDIX B. COMPONENTS OF MATRIX EQUATIONS

129

Equation (B.6) is given by the definitions Eqs. (B.4) and (B.5).
 n+ 1   n+ 1

n+ 1
vix 2
viy 2 γ − viz 2 β  n+ 12   n+ 12
n+ 12 −viy  viz α − vix γ 

 

n+ 12
n+ 12
n+ 12
viz
vix β − viy α n− 1    n− 1

n− 1
viy 2 γ − viz 2 β
vix 2
Exn
Dy Bzn − Dz Byn



  
 ηZ n e∆t  n− 12  Z n e∆t  n   n− 12
n− 12 nn
= viy  +Ey  + viz α − vix γ  −nn Dz Bx − Dx Bz mi   µ 0 mi 




n− 12
n− 12
n− 12
n
n
n
Ez
Dx By − Dy Bx
viz
vix β − viy α


−γ

1

  n+ 1 v 2
  ix 1   n+ −α viy 2 


n+ 12
1
viz
β



 γ
1

−β α
 n− 1



n− 12
n− 12
Z n e∆t nnn
2v+ mn Ex + γviy − βviz
Dy Bz − Dz Byi
 ix 1 ηZ n e∆t 

1
1
 n−

n
n− 2 n− 2nnn= viy 2 + Z me∆t
−


E
−
γv
+
αv
D B − Dx Bz 
n
y
iz
ix 
i
µ0 mni  z x


n− 21
n− 12
n− 12
Z n e∆t nnn
viz + mn Ez + βvix − αviy
Dx By − Dy Bxi(B.6)
Here, A and S are defined as follows, respectively.
1 −γ β




A= γ
1 −α


−β α
1


 n− 1

n− 12
n− 12
n
n
n
−
D
B
vix 2 + Ze∆t
+
γv
−
βv
D
B
E
z
y
y
z
x
iy
iz
mi

 ηZe∆t  n− 12 Ze∆t nn− 12
n− 12 S = viy + m Ey + αviz − γvix  −Dz Bxn − Dx Bzn iµ0 mi 



n− 21
n− 21
n− 12
Ze∆t nDx Byn − Dy Bxn
viz + mi Ez + βvix − αviy
n+ 12

vi

is expressed as follows by using A and S.

APPENDIX B. COMPONENTS OF MATRIX EQUATIONS

n+ 12

Avi

n+ 12

vi

130

=S
= A−1 S

1 + α αβ + γ αγ − β


1


2
A−1 =


αβ
−
γ
1
+
β
βγ
+
α
2
2
2
1+α +β +γ 2αγ + β βγ − α 1 + γ
2n+ 1Finally, vix 2 is obtained as follows.

1
n+ 1
vix 2 =
2
1 + α + β2 + γ2
[
{
[
]}
1
Z n e∆t
1 n n− 12
η
n− 12
2
n
n n− 2nn(1 + α ) vix +Ex + (Bz viy − By viz ) − (Dy Bz − Dz By )
mni
2
µ0
{
[
]}
1
1 n n− 12
Z n e∆t
η
n− 21
n
n
n
n n− 2+ (αβ + γ) viy +Ey + (Bx viz − Bz vix ) − (Dz Bx − Dx Bz )
mni
2
µ0
{
[
]}]
1
1
Z n e∆t
1
η
n− 12
n−
n−
+ (αγ − β) viz +Ezn + (Byn vix 2 − Bxn viy 2 ) − (Dx Byn − Dy Bxn )
mni
2
µ0

B.2

x component of E n

We derive x component of E n for Eq. (3.17) in the Chapter 3. This ismuch of the same for Eyn and Ezn .

1
E =− n n
Z eni
n

Exn =

{(

∇×B
jin −
µ0

n

}

)
× B + kB Te Z
n

n

∇nni

+

η
∇ × Bn
µ0

{( n
)
( n
) }
∂By
1
∂Bx
∂Bzn
∂Bxn
n
−
−
Bz −
Byn
µ0 Z n enni
∂z
∂x
∂x
∂y
(
)
) kB Te ∂nni
η ∂Bzn ∂Byn
1 ( n nn n+−
− n n jiy Bz − jiz By −Z enienni ∂xµ0∂y∂z

Appendix C
Neumann boundary condition
At the calculation boundary, the outflow particles are crossing simply
and the electromagnetic field is kept constant. The field at the boundary isinterpolated by only the known physical quantities of the inside cells. For an
orthogonal cell, this is simple, since the volume which has to be filled is a sliceof this cell, which is again simply an orthogonal volume. For a nonorthogonalcell, this volume is more complicated and it is diﬃcult to place the particles
inside this volume with the correct distribution function. Figure C.1 shows
the calculation boundary. The red and blue cells indicate orthogonal andnonorthogonal cells, respectively.
Here, the boundary condition adopted for the field quantities ϕ is thatthe spatial diﬀerences of the normal components ∂ϕ/∂n are set to be zero;i.e., a Neumann boundary condition is adopted. We describe two examples
for determining the physical quantity ϕ(P ) at the point P, as shown in Fig.
C.2.

131

APPENDIX C. NEUMANN BOUNDARY CONDITION

132

y

x

Figure C.1: Orthogonal and nonorthogonal cells.

n
T
R

P

U

Q B

N

L

∆x
T
∆y

N n
R

θ

P
L

Q
θ

U

B

Figure C.2: Two examples of Neumann boundary condition.

APPENDIX C. NEUMANN BOUNDARY CONDITION

133

1. From RQ = ∆x tan θ and QU = ∆y − ∆x tan θ, ϕ(Q) is expressed asfollows, where the over line indicates the distance between two points.

ϕ(Q) =

ϕ(R)(∆y − ∆x tan θ) + ϕ(U )∆x tan θ
∆y

(C.1)

From P Q = ∆x/ cos θ and Eq. (C.1), the normal component is expressed.

∂ϕ(N )
ϕ(P ) − ϕ(Q)=∂nPQ
ϕ(P ) cos θ ϕ(R)(∆y cos θ − ∆x sin θ) + ϕ(U )∆x sin θ
=
−
∆x
∆x∆y
(C.2)
From Eq. (C.2), ϕ(P ) is expressed when ∆x is equal to ∆y.

ϕ(P ) ==∆x ∂ϕ(N ) ϕ(R)(∆y cos θ − ∆x sin θ) + ϕ(U )∆x sin θ+
cos θ ∂n∆y cos θ∆x ∂ϕ(N )+ ϕ(R)(1 − tan θ) + ϕ(U ) tan θ
cos θ ∂n(C.3)

By substituing ∂ϕ(N )/∂n = 0 into Eq. (C.3), ϕ(P ) is obtained asfollows.

ϕ(P ) = ϕ(U ) tan θ + ϕ(R)(1 − tan θ)
2. From U Q = ∆x − ∆y/ tan θ and QB = ∆y/ tan θ, ϕ(Q) is expressed
as follows.
(
)
∆y
∆y
ϕ(U ) tan
+ ϕ(B) ∆x − tanθθϕ(Q) =∆x(C.4)

APPENDIX C. NEUMANN BOUNDARY CONDITION

134

From P Q = ∆y/ sin θ and Eq. (C.4), the normal component is expressed.

∂ϕ(N )
ϕ(P ) − ϕ(Q)=∂nPQ
ϕ(P ) sin θ ϕ(U )∆y cos θ + ϕ(B)(∆x sin θ − ∆y cos θ)
=
−
∆y
∆x∆y
(C.5)
From Eq. (C.5), ϕ(P ) is expressed when ∆x is equal to ∆y.

∆y ∂ϕ(N ) ϕ(U )∆y cos θ + ϕ(B)(∆x sin θ − ∆y cos θ)+
sin θ ∂n∆x sin θ()∆y ∂ϕ(N ) ϕ(U )1=++ ϕ(B) 1 −(C.6)
sin θ ∂n
tan θ
tan θ

ϕ(P ) =

By substituting ∂ϕ(N )/∂n = 0 into Eq. (C.6), ϕ(P ) is obtained asfollows.
(
)
ϕ(U )
1
ϕ(P ) =+ ϕ(B) 1 −
tan θ
tan θ

Appendix D
Normalization
All physical quantities used in 3D hybrid PIC code are normalized asfollows to reduce the calculation error and prevent the calculation from divergence.

x̃ =

x
c∆t

∆x̃ =

∆x
c∆t

Length;Mesh width;
Velocity;

ṽ =

v
c

Magnetic field;

B̃ = B ·Ze∆tmiElectric field;

Ẽ = E ·Ze∆t
cmi

Current density;j̃ = j ·c2 ∆t3
Ze

Particle density;ñ = n(c∆t)3Resistivity;

η̃ =

135

1
η
µ0 c2 ∆t

Appendix E
CFL condition
Time step ∆t and mesh width ∆x are limited by Courant-FriedrichsLewy (CFL) condition for solving partial diﬀerential equations numerically.
This CFL condition is easily derived from the wave dispersion relation. Here,
A(x, t) is expressed, where Â is complex magnitude, k is the wave numberand ω is the frequency.

A(x, t) = Â exp(ikx − iωt)Here, the spatial central diﬀerence at x0 is considered.

A(x0 + ∆x/2, t) − A(x0 − ∆x/2, t)
∆A
=
∆x
∆x
=

exp(ik∆x/2) − exp(−ik∆x/2)
A(x0 , t)∆x=isin(k∆x/2)
A(x0 , t)
∆x/2

(E.1)

By comparing Eq. (E.1) with ∂A/∂x, the wave number is expressed asfollows.

136

APPENDIX E. CFL CONDITION

K=

137

sin(k∆x/2)
∆x/2

(E.2)

This is much of the same for ω by using Ω.

Ω=sin(ω∆t/2)
∆t/2

(E.3)

The wave dispersion relation is ω 2 = c2 k 2 . c is the speed of light. Fromreplacing k and ω by Eqs. (E.2) and (E.3), Ω2 = c2 K 2 . Since the maximum
wave number is kmax = π/∆x,(2

sin (ω∆t/2) =

c∆t
∆x

)2
(E.4)

is obtained. When c∆t/∆x > 1, ω becomes complex and the calculation
becomes unstable and diverges. Since hybrid PIC code uses the implicit
method not the explicit method, c is replaced by the Alfven speed VA . Then,the CFL condition is ∆t < ∆x/VA .