Kinetic simulation of spherically symmetric collisionless plasma in the inner part of a cometary coma

Uppsala University

Kinetic simulation of spherically symmetric collisionless plasma in the

inner part of a cometary coma

Student: Pavel Dogurevich Supervisor: Anders Eriksson

July 11, 2019

Contents

1 Introduction 1

2 Cometary plasma environment 1

3 Model 2

3.1 Overview . . . . 2

3.2 Neutral gas . . . . 3

3.3 Plasma . . . . 4

3.4 Absorbing nucleus . . . . 5

3.5 Electric field . . . . 5

3.6 Electrostatic potential . . . . 6

3.7 Equations of motion . . . . 6

4 Simulation 7 4.1 Overview . . . . 7

4.2 Motion . . . . 7

4.3 Ionization . . . . 8

4.4 Plasma absorption . . . . 9

4.5 Metaparticles . . . . 9

4.6 Normalization . . . . 10

4.7 Summary . . . . 11

4.8 Implementation . . . . 13

5 Results 13 5.1 Simulation parameters . . . . 13

5.2 Charge of nucleus . . . . 15

5.3 Plasma number density . . . . 16

5.4 Electrostatic potential . . . . 18

5.5 Electron motion . . . . 20

6 Discussion 21

2 COMETARY PLASMA ENVIRONMENT 1

Figure 1: Plasma environment of a highly active comet. Credit: NASA/JPL

1 Introduction

This report describes the work done for a short research project concerned with the simulation of cometary plasma. The next section gives a brief overview of the cometary plasma environment. It is followed by the descrip- tion of the model and its numerical application as the simulation. The final sections of the report present the results obtained from the simulations and discuss the applicability and possible improvements to the simulation and the model.

2 Cometary plasma environment

Comets are small bodies made of dust and ice orbiting the Sun. When a

comet come closer to the Sun, the ices in the cometary body (the nucleus)

are evaporated by the solar radiation. This effect creates a large, expanding

envelope of gas and dust around the nucleus: the coma. The constituent

particles of the coma become photoionized by the solar radiation and via

other processes, leading to a complex plasma environment surrounding the

comet that is yet to be fully understood. An example of a cometary plasma

environment is depicted in Figure 1. For an in-depth discussion of the plasma

3 MODEL 2

environment consult Odelstad (2018).

The main focus of the project is the inner part of the cometary iono- sphere, or the inner coma, which roughly corresponds to the diamagnetic cavity. Neither solar wind particles nor external magnetic fields can pene- trate into this region. The environment in the inner coma is dominated by unmagnetized dusty plasma produced from the neutral gas and dust par- ticles lifted from the nucleus. Main ionization processes occurring in the region are photoinization by solar radiation, photoelectron impact ionization and charge transfer between ions and neutral molecules. The inner coma is relatively simple in comparison to the other regions of the plasma envi- ronment, which could allow a simple model to produce physically relevant results. This project attempts to do exactly that.

3 Model

We chose to apply one of the simplest, idealised analytical models to describe the conditions in the inner coma. This section outlines the different parts of the model. Details concerning the numerical application and time evolution are contained in Section 4.

3.1 Overview

This model describes the production of the neutral gas at the nucleus, the ionization of the neutral molecules, and the interaction of the ions and elec- trons via the electric field. Figure 2 depicts the main parts of the model graphically.

We applied several major approximations to simplify the processes in the inner coma.

• The plasma is considered collisionless and unmagnetized. The only interaction between the ions and electrons is the electrostatic force.

The plasma does not affect the neutral gas.

• Only a single ionization process is included in the model. The ion- ization rate is constant and limited to a spherical volume surrounding the nucleus. Outside of the sphere the ionization rate is 0. The elec- trons and ions are produced at constant energies. That is, the energy distribution function is the delta function.

• The neutral gas density profile is spherically symmetric and time inde-

pendent.

3 MODEL 3

Figure 2: A sketch of the model. Pictured are the nucleus of radius R _n and charge q _n , the ionization region of width ∆R _i and the simulation particles with positions r 1 , r 2 , r 3 .

• The electric field generated by plasma particles is spherically symmet- ric. The charge of a particle is considered to be distributed over a sphere centered at the origin (see Figure 2). Although plasma parti- cles technically move in 3 dimensions, this approximation allows us to reduce the equations of motion to the radial coordinate and the radial velocity. It also greatly simplifies the equation for the electric field.

3.2 Neutral gas

The nucleus is modeled as a perfect sphere of radius R _n located at the origin.

The neutral molecules are produced at the surface of the nucleus at the total production (outgassing) rate of Q (in particles per unit time). The newly produced molecules are distributed uniformly over the surface of the nucleus.

They travel radially outwards with constant velocity u. Both Q and u are constants, and the model does not include processes that affect the neutral gas. Therefore, the density of the neutral gas is given by

n(r) = Q

4πr ² u . (3.1)

3 MODEL 4

3.3 Plasma

Several main processes contribute to the production of ions in the inner coma in reality. This model includes a single ionization process with constant and uniform ionization rate, which can be used to reproduce processes such as photoionization. Photoionization occurs when a neutral molecule absorbs a high energy photon, producing a singly-ionized ion and an electron (we neglect processes that produce more than one electron). In the model, all the extra energy of the photon is transferred to the electron, which is ejected at a random direction. The ion continues to travel radially outwards with the initial velocity equal to the velocity of the neutral molecule. The electron and ion are created with the initial kinetic energy of E _0,e and E _0,i , respectively.

The number of ions created per unit time depends on the ionization rate of the neutral gas ν. We assume that the ionization rate is constant and uniform in the spherical shell extending from ∆R _n to R _n + ∆R _i . The ionization rate outside the ionization volume is 0. Then the ion production rate (particles per unit time and unit volume) for r ∈ [R _n , R _n + ∆R _i ] is given as

p(r) = ν n(r) = ν Q

4πr ² u . (3.2)

The number of ions produced in time ∆t in the ionization shell is

∆P (∆t) =

Z t+∆t t

Z

V

p(r) dV dt = ν Q

u ∆t∆R _i . (3.3)

Note that the number of ions produced in the shell per time unit does not depend on the distance from the nucleus.

The maximum ion number density can be estimated in a simple way using the ion production rate, if the electrostatic interaction in the plasma is neglected (”neutral plasma”). Then the total flux of ions through a sphere of radius r at time t < ^r−R _v

i

is given as Γ(r, t) = 4π

Z r r−v

t

p(r) r ² dr = ν Q

u v _i t, (3.4)

where v _i is the initial velocity of the ions, which would continue travel radially outwards with a constant velocity in neutral plasma. Not all ions produced closer than r to the nucleus in time t have reached r, hence the flux is increasing with time. But when t is large enough for all such ions to reach r (in other words, t ≥ ^r−R _v

i

), the flux is given by Γ(r) = 4π

Z r R

p(r) r ² dr = ν Q

u [r − R _n ] . (3.5)

3 MODEL 5

In this case the flux is constant in time and depends only on the distance from the nucleus. The equation for the number density is then

n ⁰ (r) = Γ(r)

4πr ² v _i = ν Q 4πruv _i



1 − R _n r



, (3.6)

which has the maximum at 2R _n . So the maximum number density is given as

n ⁰ _max = ν Q

16πR _n uv _i . (3.7)

This number is useful for choosing time and length scales of the simulation, as will be described in Section 4.6.

3.4 Absorbing nucleus

The model allows the nucleus to carry charge. Initially, the charge of the nucleus q n is 0. Any plasma particle that reaches the surface of the nucleus is absorbed into the nucleus. The particle is considered to be removed from the model, while its charge is added to the charge of the nucleus q _n . Since the charge is fully transferred from the removed particles to the nucleus, the total charge in the model is conserved (the total charge is zero).

3.5 Electric field

Ions and electrons in the plasma give rise to an electric field. Due to the spherical symmetry of the model, Gauss’ law can be used to obtain the expression for the electric field at the position of each particle (Boella et al.

2018).

Consider a sequence of ions and electrons arranged by their distance to the nucleus in the ascending order. Then, the component of electric field in the radial direction at particle i is

E _i (r) = 1 4πε ₀ r ²

i−1

X

j=1

q _j + q _n

!

. (3.8)

Denoting k = _4πε ¹

0

and Q i = P i−1

j=1 q j + q n , the equation becomes E _i (r) = k Q i

r ² . (3.9)

The electric field does not depend on the distances between the particles,

only on the distance of the particle to the origin.

3 MODEL 6

In contrast to the equation for the electric filed in the shell method from Boella et al. (2018), we decided against including the self-interaction term + ¹ ₂ q _i . Although it is justifiable within the ideal model of shell-like particles, real particles do not interact with themselves. We believe that the omission of the term brings the model closer to reality.

3.6 Electrostatic potential

Using the equation for the electric field, the equation for potential can be obtained. Consider N particles ordered by increasing r. Taking potential at infinity as 0, potential at the most distant particle N is

φ _N = Z ∞

r

E _N dr = kQ _N Z ∞

r

1

r ² dr = k Q _N

r _N . (3.10)

Similarly, the potential at particle i is φ _i =

Z ∞ r

E _i dr

=

N −1

X

j=i

Z r

r

E _j dr + Z ∞

r

E _N dr

| {z }

φ

=

N −1

X

j=i

kQ _j  1 r _j − 1

r _j+1



+ k Q _N

r _N . (3.11)

Hence, the difference in potential between two adjacent particles is φ _i − φ _i+1 =

N −1

X

j=i

kQ _j  1 r _j − 1

r _j+1



−

N −1

X

j=i+1

kQ _j  1 r _j − 1

r _j+1



= kQ _i  1 r i

− 1 r i+1



. (3.12)

Using the formula above the potential can also be defined in a recurrent form

(φ _N = k ^Q _r

N

,

φ _i = φ _i+1 + kQ _i 

1 r

− _r ¹

i+1



, i < N. (3.13)

3.7 Equations of motion

The only force acting on the plasma in the model is due to the electric field.

Because of the symmetry, the motion of the plasma particles is essentially

4 SIMULATION 7

one-dimensional in the radial direction. The electric field acts along the radial direction, hence the angular momentum of particles is conserved. Only electrons have a non-zero angular momentum, as they are produced with the velocity vector oriented in a random direction. The ions, like the neutral molecules, have angular momentum of zero.

We applied the shell method from Boella et al. (2018) to describe the motion of the plasma. As usual, the particles are ordered by the increasing r. Denoting the angular momentum as p _ϕ and mass of the particle as m, the equations of motion for a single particle are given below (the subscript i is omitted for convenience).

dr

dt = v r , (3.14)

dv _r dt = qE

m + p ² _ϕ

m ² r ³ = k q m r ²

i−1

X

j=1

q _j + q _n

!

+ p ² _ϕ

m ² r ³ (3.15) The presence of the 1/r ³ term could cause issues when r approaches 0. But in our case the nucleus prevents particles from reaching the origin.

4 Simulation

4.1 Overview

A computer simulation code was implemented using the model described in Section 3 to observe the time-development of plasma under chosen condi- tions. The code implements one step (or iteration) of the simulation, which advances the state of the simulated system in time by ∆t. A single simulation runs for the required number of iterations τ _max . One iteration mainly consists of production of new plasma particles, evolution of the position and velocity of plasma particles in time, and absorption of plasma by the nucleus. The subsections below describe the implementation of the simulation in detail.

4.2 Motion

The equations of motion require the particles to be ordered by their posi- tion, which changes in time, leading to particles overtaking each other and changing the order. This makes it impossible to find an analytical solution.

Instead, the simulation uses the leapfrog numerical integration method to

iteratively evolve the position and velocity of particles. This method is easy

to implement and it suits the problem very well. It is applicable to second-

order differential equations. It conserves the dynamical energy of the system.

4 SIMULATION 8

And it is stable for oscillatory motion when the time step is smaller than the period of oscillations.

The idea behind the method is that the position is shifted by half time step relative to the velocity. Mathematically, the iteration step τ looks as

r _{τ +1/2} = r _{τ −1/2} + v _τ ∆t, (4.1)

v _{τ +1} = v _τ + a(r _{τ +1/2} ) ∆t. (4.2) At each iteration of the simulation the values for velocity correspond to the current simulation time, while the values for position are a half time step behind. Note that the method requires that both velocity and position do not depend on themselves, which is true for our model.

Leapfrog method imposes an important restriction on the value of the time step ∆t. The time step must be constant and less or equal to 2/ω, where ω is the angular frequency of the oscillations. The electron plasma oscillations are the oscillations with the highest frequency in the plasma.

Therefore, we use the electron plasma frequency ω _pe to limit ∆t.

A complication of our model is the need for the particles to be ordered by r for the acceleration term to evaluate correctly. Thus a re-ordering step is included in the implementation just after the position update, before the updated velocity is evaluated.

4.3 Ionization

The simulation is initially devoid of plasma particles. At the beginning of each simulation iteration new ions and electrons are created via the ionization process.

During each ionization step ∆N ion-electron pairs are created at positions chosen randomly in [R _n , R _n + ∆R _i ]. The ion and electron from the same pair are created at the same position.

Conforming to the model, the ions are created with zero angular momen- tum and with initial radial velocity v _r = p2E _0,i /m _i . The ions can never gain angular momentum as it is constant in time.

The electrons are created with their velocity oriented in a random direc- tion. That is, denoting the angle between the radial vector and the velocity as α, cos α is randomly chosen from [−1, 1]. Given the total initial veloc- ity v = p2E _0,e /m _e , the initial radial component of the electron velocity is v _r = v cos α, and the angular momentum is p _ϕ = r ² ₀ v ² sin ² α (where r ₀ is the position at creation).

The final part of the ionization process is to advance the velocities of the

newly created particles by a half of the time step, as required by the leapfrog

integration method.

4 SIMULATION 9

4.4 Plasma absorption

The final part of the simulation step is the absorption of plasma into the nucleus. The nucleus is represented as a sphere of radius R _n fixed at the origin, acting as the lower spatial boundary for the simulation particles.

Initially, the nucleus is electrically neutral. During the simulation, any particle that reaches the nucleus is absorbed by it: the charge of the particle is transferred to the nucleus, and the particle is removed from the simulation.

Mathematically, the evolution of the charge of the nucleus q _n can be expressed as

q _{n,τ +1} = q _n,τ + X

r

≤R

q _i . (4.3)

4.5 Metaparticles

Substituting typical values from the observations (for example, from comet 67P) into the formula of the ion production rate ∆P gives very large num- bers: more than 10 ²⁰ ions per second. Simulation of this amount of particles is not feasible on a regular PC (and, arguably, on any computer existing nowadays). To make the computation possible, the simulation code runs with metaparticles.

Each metaparticle represents N _s physical particles, which is referred to as the weight of metaparticles. The weight is constant in time and identical for all metaparticles. Consequently, the transition to metaparticles keeps the equations of motion unmodified, as N _s is canceled in q/m fraction, but an N _s factor is introduced into the equation for the electric field. The metaparticles behave identical to individual physical particles dynamically, but produce an electric field which is N _s times stronger.

The weight is simply the ratio of the number of ions ∆P to the number of metaparticles ∆N produced per time step in the ionization volume,

N _s = ∆P

∆N = ν Q u

∆t ∆R _i

∆N . (4.4)

The weight partially depends on the physically-based quantities ν, Q and u, and partially on the simulation parameters ∆t, ∆R _i and ∆N , which can potentially take arbitrary values. Ideally, large values would be chosen for the time step (to make the simulation time pass faster) and for the ioniza- tion radius (to better reflect the physical conditions around the nucleus), while keeping ∆N value low. However, since the electric field is proportional to the weight, the potential difference between two adjacent metaparticles could become so large that most particles would not be able to pass trough.

Metaparticles would effectively become bound to each other, which is not

4 SIMULATION 10

representative of the conditions we aim to simulate. Hence, an upper limit must be imposed on the weight, limiting possible values of the simulation parameters.

The potential drop between two adjacent metaparticles has to be much less than their initial energy, to allow the particles to move freely and the simulation to evolve smoothly. Using Equation 3.12, the potential difference between adjacent metaparticles separated by distance ∆r is calculated as

∆φ = N _s k |Q|  1

r − 1

r + ∆r



. (4.5)

First, consider which values of r and ∆r would result in the maximum po- tential drop. The value of r should be the smallest possible, which is R n . For ∆r, the maximum possible value should be chosen. We estimated this value by considering the case where particles are stationary and distributed uniformly over the ionization radius. Every simulation step ∆N ions are produced, so the separation between them is _τ ^∆R

0

∆N , where τ ₀ is the number of time steps passed since the start of the simulation. The resulting equation for the maximum expected potential drop is

∆φ _max = k e ν Q u

∆t ∆R _i

∆N

1 R n

− 1

R _n + _{τ ∆N} ^∆R

!

, (4.6)

and the limiting condition can be expressed as E _0,min  ∆φ _max .

The limiting condition gives insight into the possible values and coupling of the simulation parameters. From the equation, it follows that to maintain the constant value of ∆φ _max the changes in the time step or ionization radius must be accompanied by the proportional change in ∆N . If the values of ∆t and ∆R _i are fixed, the minimum ∆N can obtained by solving the equation using a numerical method, as the equation is transcendental in ∆N . In this case, the value of τ ₀ should be chosen to be much less than the total number of simulation steps.

4.6 Normalization

Quantities involved in the simulation are normalized by the scaling param- eters instead of being used as is. Normalization reduces the floating point precision errors during the calculation, as the quantities involved become closer to unity. It also simplifies the interpretation and comparison of the resulting data.

The normalization process is simple: for each variable in the model a

scaling parameter is chosen, and the variable is divided by the scaling pa-

rameter. If the variable is x and the scaling parameter is x _s , then the scaled

4 SIMULATION 11

Table 1: Scaling parameters

Parameter Value

Distance r s λ D =

q E

ε

n

e

Time t _{s ω} ¹

pe

= q

m

ε

n

e

Velocity v _s

q E

m

Mass m _s m _e

Charge q _s e

Energy E _s E _ref

Electrostatic potential φ _s ^E

_e Electrostatic coefficient k s φ

r

q

=

q E

ε

n

e

variable is x ⁰ = x/x _s . This process is applied to input parameters and to the model equations.

The four main variables in the model are the radial coordinate, time, mass and charge. The choice of scaling parameters for mass and charge is straightforward: the electron mass m _e and elementary charge e. For space and time, we selected the characteristic plasma quantities: Debye length λ _D and the inverse of electron plasma frequency w _pe . These parameters require reference values of number density n _ref and energy E _ref for evaluation. The maximum neutral plasma number density n ⁰ _max (Equation 3.7) is a good value for n _ref , since the actual number density in a simulation is expected to be lower than n ⁰ _max . This ensures that the time and space scales are less or equal to the actual plasma period and Debye length in the simulation. Finally, it is practical to choose the initial electron energy as the reference energy, as the normalized electron energy then becomes unity.

Other scaling parameters such as velocity and electrostatic potential scales are derived from the main scaling parameters (i.e. v _s = r _s /t _s ). Table 1 pro- vides the complete list of the scaling parameters used in the simulation.

The application of normalization does not change the form of the equa- tions of motion of the plasma.

4.7 Summary

Table 2 lists all input parameters that control the initial state of the simu-

lation. Most of the parameters are based on physical quantities that can be

4 SIMULATION 12

Table 2: Input parameters Parameter Description

R n Radius of the nucleus.

Q Neutral gas production rate.

ν Ionization frequency of the neutral gas.

u Neutral gas velocity.

m i Mass of an ion.

E _0,e Initial electron energy.

E 0,i Initial ion energy.

E _ref Reference energy.

n ref Reference number density

∆R _i Ionization radius.

∆N Metaparticle production rate.

∆t Time step.

τ max Number of iteration steps.

measured in the inner coma of a comet. The bottom four (∆R _i , ∆N, ∆t, τ _max ) are the consequences of the numerical application of the model. These pa- rameters are not bound to observable physical quantities, and only concerned by computational performance and intent of the simulation.

The input parameters initialize a particular run of the simulation, which then normalizes the inputs (Section 4.6) and runs for τ _max number of iter- ations. The outline of a single iteration is listed below. Note that all the equations and variables involved are normalized.

1. Create ∆N ion-electron pairs of metaparticles at random positions with energies E _0,i and E _0,e . Advance v _r of the particles by ∆t/2.

2. Update positions of each particle

r _{τ +1/2} = r _{τ −1/2} + v _r,τ ∆t.

3. Order particles by increasing r.

5 RESULTS 13

4. Update velocities for each particle

v _{r,τ +1} = v _r,τ +

"

k N _s q m r _{τ +1/2} ²

i−1

X

j=1

q _j + q _n,τ

!

+ p ² _ϕ m ² r _{τ +1/2} ³

#

∆t.

5. Absorb particles with r ≤ R _n into the nucleus q _{n,τ +1} = q _n,τ + X

r

≤R

q _i .

4.8 Implementation

The simulation code was implemented with Python 3 using Numpy, SciPy and Matplotlib packages. The source code is available on GitHub:

https://github.com/pdogurevich/spherical-coma

5 Results

This section describes the different simulation runs and the results we ob- tained from them.

5.1 Simulation parameters

The input parameters for the simulations were based on data from the ob- servations of the comet 67P during Rosetta mission. Specifically, the values around January 2015 were chosen from Heritier et al. (2018) as the references.

Although most of the parameter values are representative of the physical conditions at comet 67P, the mass of the ions is not. In reality, the plasma in the inner coma is dominated by singly-ionized water ions. We chose a much lighter ion while keeping the initial kinetic energy unchanged. The ions became much faster, making their transport time across the ionization region smaller. This allows for quicker saturation of the number density close to the nucleus, bringing the simulation to a quasi-stable state faster and decreasing the required simulation run time. It is important to mention that the change in ion mass does not affect the neutral gas velocity or number density.

The change in ion mass is accompanied by the change in ion number

density. According to the neutral plasma number density (Equation 3.6),

the faster the ions are moving the lower the number density will be. It

follows that the number density observed in the the saturated state of the

simulation will be lower for lighter ions.

5 RESULTS 14

Table 3: Typical values used for the input parameters for the simulations and the scaling parameters derived from them.

Parameter Value Quantity Value

R _n 1 km r _s 1.5 m

Q 10 ²⁶ s ⁻¹ t _s 1.1 × 10 ⁻⁶ s

ν 10 ⁻⁷ and 1.9 × 10 ⁻⁶ s ⁻¹ v _s 1.3 × 10 ³ km s ⁻¹

u 1 km s ⁻¹ m _s m _e

m _i 100 m _e q _s e

E _0,e 10 eV E _s 10 eV

E _0,i 0.1 eV φ _s 10 V

E _ref E _0,e N _s 1.1 × 10 ¹² and 2.4 × 10 ¹³

n _ref 250 cm ⁻³ n ⁰ _max 11 and 202 cm ⁻³

∆R _i 10 km ∆φ _max 1.6 × 10 ⁻⁴ and 3.4 × 10 ⁻³ eV

∆N 10 per timestep t 0.06 s

∆t 0.1 t _s = 1.1 × 10 ⁻⁷ s τ _max 500’000

However, the input parameters can be modified to increase n ⁰ _max to the same value that would be expected for water ions. This change would make the results more representative of the conditions in the inner coma. But it would increase the maximum expected potential drop (Equation 4.6), which could lead to noisy and even incorrect results.

In the end we ran two types of simulations: with low and high ion num- ber density. For the high density simulation the ionization rate is increased proportionally to the increase in the ion velocity due to lighter mass. Ta- ble 3 shows the typical input and scaling parameters for both types of the simulations. Additionally, other characteristic simulation quantities are in- cluded for reference. Note that we used the same reference number density n _ref for both types of simulations, leading to same scaling parameters applied to both.

Consider the values of the input parameters that are not directly based on

observations. The time step ∆t must be smaller than t _s /2 for the integration

5 RESULTS 15

Figure 3: The evolution of the charge of the nucleus in time in low (left) and high (right) ion number density simulations.

method to be stable. The value of t _s /10 was chosen to ensure accuracy of the results. The ionization radius were chosen to be sufficiently small to avoid high metaparticle weights, but still making sense physically for the physical extent of the inner coma. The metaparticle production rate ∆N was fixed at a small number to decrease the computation complexity, but kept high enough to avoid high potential difference between particles. The number of iterations τ _max is selected such that the simulation finishes in reasonable time. For the hardware the simulations ran on, 500’000 iterations finished in about five days. The values for these parameters were chosen empirically, after running several simulations and tweaking the inputs.

5.2 Charge of nucleus

Figure 3 shows the evolution of the charge of the nucleus in time. The plot can be separated in two stages. At the beginning of the simulation, the nucleus quickly accumulates negative charge by absorbing electrons. Ions initially travel away from the nucleus, so the nucleus needs time to attract them via electric field. Electrons are faster than ions and may travel towards the nucleus after ionization, making them likely to be absorbed. Then the second stage occurs when the nucleus obtains the charge strong enough to attract ions. At the same time the negative charge repels electrons, decreasing the number of electrons absorbed in time. The charge begins to oscillate, while the mean value appears to slowly approach a constant: −2.4 × 10 ¹⁴ e and

−2 × 10 ⁻¹⁵ e for the low and high density simulations respectively.

During most of the simulation run the charge of the nucleus appears to oscillate with a period similar to the period of electron plasma oscillations.

In high density simulation, n _e ≈ 175 cm ⁻³ close to the nucleus, giving the

5 RESULTS 16

Figure 4: The oscillations of the charge of the nucleus over a short time period in low (left) and high (right) ion number density simulations..

period of oscillation T pe = q

m



n

e

= 8.4 × 10 ⁻⁶ s = 7.7 t s . In low density simulation n e ≈ 1.1 cm ⁻³ , giving T pe = 1.1 × 10 ⁻⁴ s = 97 t s . Both values are close to the periods of oscillations displayed in Figure 4. This effect is likely caused by the oscillations of the electrons close to the nucleus. As the electrons are displaced toward the nucleus during the oscillations, they hit the nucleus and get absorbed.

The difference in the slope of the plots between the low and high density simulations in Figure 3 can be explained by the higher metaparticle weight in the high density simulation, in addition to the difference in number densities.

The dynamics of metaparticles is identical in both simulations, but in the high density simulation metaparticles produce an electric field that is about 20 times stronger. Although as many electron metaparticles reach the nucleus in the same period of time, the electric field produced by the nucleus is stronger. It starts attracting ions earlier, and the oscillations begin earlier in the simulation.

5.3 Plasma number density

The ion and electron number density values follow each other closely, as evident in Figure 5. The only exception is the higher ion number density close to the nucleus, which is caused by the negatively charged nucleus attracting ions and repelling electrons. The extra ions shield the more distant particles from the effects of the nucleus charge. Thus the plasma in the simulation is quasi-neutral, as expected from real plasma.

The deviations of number densities from the neutrality are within expec-

tations. The deviations along the length L are expected to be proportional

5 RESULTS 17

(a)

Figure 5: Number densities of ions and electrons over the ionization region at the end of the simulations. The deviation from the neutrality is plotted in green with the vertical axis located on the right of the plot. Left: low ion number density simulation. Right: High ion number density simulation.

to (λ _D /L) ² for L  λ _D (for the derivation consult e.g. Wiesemann (2014)).

In our case,

n _i − n _e n _i + n _e    

≈ 2  λ _D L

 2

. (5.1)

In the low density simulation, at r = 3000 the formula gives 1.7 × 10 ⁻² , while the deviations are 2 × 10 ⁻² . In the high density simulation, the expectation at r = 1900 is 8 × 10 ⁻⁴ : three times lower than the actual deviations of 2.5 × 10 ⁻³ at the same position. In both cases, the values are reasonably close.

The number density for both types of simulations decays approximately as 1/r, as evident from Figure 6. The neutral plasma model (Equation 3.6) predicts the same when r  R _n , so this behavior is expected.

In comparison to the neutral plasma model, the number densities are

lower, and their peaks are closer to the nucleus. Only in the high density

simulation the number density reaches the values expected from the neutral

plasma model, around the peak. Otherwise, the densities are 2 - 6 times

lower. Although the difference is to be expected (plasma is not electrically

neutral, after all), the gap would be less if the simulations ran for longer

time to allow slow ions to saturate the ionization region and increase the

number density. Comparison of Equations 3.4 and 3.5 reveal that the number

density would not reach the maximum value at a position r until time t =

(r − R _n )/v _i had passed. Since we are interested in the ionization region,

the ion transport time required for the saturation of the the whole region is

5 RESULTS 18

Figure 6: Best fit to the number density profile and comparison to the neutral plasma model profile (Equation 3.6). Left: low ion number density simula- tion. Right: high ion number density simulation.

∆R _i /v _i . Unfortunately, within the short time available for this project, we could only run the simulations for a tenth of this time.

The transport time across the ionization region is 0.005 s for electrons and 0.5 s for ions. The total simulation run time is 0.06 s, meaning that the electrons have enough time to cross the region ten times. The ions, however, only travel one-tenth of the region. Taking in account that they travel only in radial direction, at the end of a run the only part of the simulation that contains a saturated number density of ions is about 1 km away from the nucleus.

5.4 Electrostatic potential

The electrostatic potential oscillates rapidly, as Figure 7 demonstrates. The oscillations are generally quicker closer to the nucleus, following the number density profile. This oscillations are the electron plasma oscillations, with the period given by the formula

T _pe = 2π/ω _pe = 2π r m _e  ₀

n e e ² . (5.2)

The actual periods of potential oscillations in the data are in agreement with

the formula. For example, in the low density simulation at r = 3000 the

oscillation period is approximately 90 t _s , while T _pe = 83 t _s . The period is

7.7 t s at r = 640 in the high density simulation, with respective T pe = 7.2 t s .

The amplitude of oscillations is much higher in the high density simulation,

probably due to a stronger electric field produced by metaparticles.

5 RESULTS 19

Figure 7: The mean and the standard deviation of the electrostatic potential in time. The mean is calculated at the end of the simulation over 5000 and 1000 time steps for low and high density simulations, respectively. For 4 selected positions the evolution of the potential in time over a short period at the end of the simulation is presented at the bottom. Left column: low ion number density simulation. Right column: high ion number density simulation.

The slope of the potential is approximately logarithmic (Figure 8). This agrees with the expected potential, assuming the number density is propor- tional to 1/r. Given the equation of motion for the electron gas

m e n e

dv

dt = −∇p e − en e E, (5.3)

applying the ideal gas law p _e = n _e kT _e and assuming equilibrium state ^dv _dt = 0 results in

E = − kT _e

en _e ∇n _e . (5.4)

Since n _e ∝ 1/r, it follows from Equation 5.4 that E ∝ 1/r as well. Then, as

E = −∇φ, the potential φ is proportional to ln r.

5 RESULTS 20

Figure 8: Best fit to the mean potential over the ionization region. Left: low ion number density simulation. Right: high ion number density simulation.

5.5 Electron motion

We injected several test electrons into the simulations to gain an under- standing of the electron motion on the level of individual particles. These electrons did not affect the simulation in any way: they did not interact with the other particles and the nucleus could not absorb them. The evolution of the electrons was observed in two stages. First, the motion was recorded while the simulation was running. Then, the electrons were evolved in a stopped simulation, in its frozen final state. The second stage is much less expensive computationally, so it allowed us to observe the electron motion for much longer time. Figure 9a shows the path of a test electron in the low density simulation plotted in phase and configuration spaces. The other test electrons behaved similarly.

Consider the trajectory of the test electron in Figure 9a. The path in phase space is quasiperiodic. But it is not forming a closed path while the simulation is running, because of the changing electric potential. Conversely, the path is closed for the motion of the electron in the frozen potential. In addition, the changes of potential on short time scales due to plasma oscilla- tions make the orbit of electron deformed, jagged. In contrast to the smooth Kepler orbits in Figure 9b, which arise from a smooth 1/r potential. The trajectory of the electron in configuration space is periodic as well. However, it does not close even in the frozen potential. According to Bertrand’s theo- rem (Santos, Soares, and Tort 2007), the only potentials that allow for closed orbits are proportional to r ² or 1/r. Therefore, the open orbits of the electron are due to the near-logarithmic electrostatic potential in the simulation.

The precession of the electron orbit in configuration space of Figure 9a

indicate that the spherical symmetry is an acceptable approximation, as over

6 DISCUSSION 21

(a) The trajectory of a test electron in phase and configuration spaces.

(b) Reference Kepler orbits in phase and configuration spaces.

Figure 9: The time evolution of a test electron in time in the low density simulation, along with reference Kepler orbits. The electron is created with v _r = 0. The solid lines indicate the data obtained during the final 200’000 time steps of the simulation. The dashed lines are the evolution of the elec- tron in the frozen potential after the simulation finished. Different colors were used for the electron path to make it easier to follow its evolution.

the time of many periods the path will cover a circular region in space.

6 Discussion

The results of the simulations are promising. Even with the significant simpli-

fication of the spherically symmetric electric field and the others, the model

reproduced the expected plasma behavior such as quasi-neutrality, plasma

oscillations, the logarithmic profile of the electrostatic potential. This gives

hope that the results of the simulations can be potentially applied to real

REFERENCES 22

plasma.

However, the simulations run time was too short. Not enough time has passed for the ions to saturate the ionization region (see Section 5.3), so the state of the simulation has not fully developed. Additionally, the run time was too short for low frequency phenomena to manifest, such as fluctuations of kinetic energy populations of electrons.

To increase the run time, improvements in the implementation and the model are required. The major flaw of the current implementation is its inability to benefit from parallelization. Running the simulation on several computational cores will make it possible to reach ion transfer time in a mat- ter of days. From the model perspective, decoupling the electron and ion time scales or using ergodic integral invariants might improve the computational time significantly.

The artificial ionization boundary is not realistic. Possible improvement would be to increase the ionization range to a high number (a few thousand km), so it would have less effect on the inner region of the simulation. An- other possibility is to remove the boundary all together, replacing it with a metaparticle weight that is dependent on the distance from the nucleus.

However, this idea requires further investigation to asses its feasibility.

Finally, the model can be made more complex to make it more realis- tic. The first obvious choice is to introduce Maxwellian energy distribution of newly produced electrons and ions instead of constant energy. Further, adding more plasma processes to the model, such as charge transfer between ions and neutral molecules, would potentially produce better results. The downside of the improved realism is, of course, higher computational com- plexity.

As the conclusion for the project, the model and the simulation are promising. But more effort is required to produce potentially physically relevant results.

References

Boella, E. et al. (2018). “Gridless particle technique for the Vlasov-Poisson system in problems with high degree of symmetry”. In: Computer Physics Communications 224, pp. 136–143. doi: 10.1016/j.cpc.2017.11.004.

Heritier, K. L. et al. (2018). “Plasma source and loss at comet 67P during

the Rosetta mission”. In: A&A 618, A77, A77. doi: 10 . 1051 / 0004 -

6361/201832881.

REFERENCES 23

Odelstad, E. (2018). “Plasma environment of an intermediately active comet.

Evolution and dynamics observed by ESA’s Rosetta spacecraft at 67P/Churyumov- Gerasimenko”. PhD thesis. Uppsala University.

Santos, F. C., V. Soares, and A. C. Tort (2007). “An English translation o

Bertrand’s theorem”. In: arXiv e-prints. arXiv: 0704.2396 [physics.class-ph].

Wiesemann, K. (2014). “A Short Introduction to Plasma Physics”. In: arXiv

e-prints. arXiv: 1404.0509 [physics.plasm-ph].