Plasma Simulation in the Limit of Strong Radiation Reaction Master’s thesis in Master Engineering Mathematics and Computational Science Nils Tornberg DEPARTMENT OF PHYSICS CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2022 www.chalmers.se www.chalmers.se Master’s thesis 2022 Plasma Simulation in the Limit of Strong Radiation Reaction Nils Tornberg Department of Physics Chalmers University of Technology Gothenburg, Sweden 2022 Plasma Simulation in the Limit of Strong Radiation Reaction Nils Tornberg © Nils Tornberg, 2022. Supervisor and Examiner: Arkady Gonoskov, Department of Physics Master’s Thesis 2022 Department of Physics Chalmers University of Technology SE-412 96 Gothenburg Telephone +46 31 772 1000 Typeset in LATEX, template by Kyriaki Antoniadou-Plytaria Gothenburg, Sweden 2022 iv Plasma Simulation in the Limit of Strong Radiation Reaction Nils Tornberg Department of Physics Chalmers University of Technology Abstract Plasma is encountered in many situations and thus investigating plasma properties is of interest. Simulating plasma dynamics can be computationally demanding, especially when high-intensity electromagnetic fields are present, such as in modern high-intensity lasers. A common approach to such simulations is to use a particle-in- cell code. In high-field situations the force of radiation reaction becomes important. Recently it has been shown that in the limit of strong fields this radiation reaction causes particles to travel along a radiation free direction. In this work a particle-in- cell code based on the principle of particles following this radiation free direction is presented and evaluated. The code is publicly available at https://github.com/ TornbergNils/rfd_solver. The code is found to produce sensible predictions in line with what is found by more conventional particle-in-cell codes. Keywords: Plasma, Particle-in-Cell, PIC, Radiation reaction v https://github.com/TornbergNils/rfd_solver https://github.com/TornbergNils/rfd_solver Acknowledgements This work could not have been produced without the help of my advisor Arkady Gonoskov. His guidance has ben invaluable in the course of this work. I also want to thank my friends and family for their support. Nils Tornberg, Gothenburg, June 2022 vii Contents List of Figures xi 1 Introduction 1 2 The Radiation Free Direction and Plasma 3 2.1 Radiation-free direction . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.2 Quantum limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.3 Strong fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.4 Radiation reaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.5 Plasma properties and waves . . . . . . . . . . . . . . . . . . . . . . . 7 2.6 The plasma z-pinch . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3 Structure of a particle-in-cell code 11 3.1 Particle-in-cell approach to plasma simulation . . . . . . . . . . . . . 11 3.2 Advancing the field . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.3 Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.4 Moving the particles using the Boris pusher . . . . . . . . . . . . . . 14 3.5 RFD pusher and treatment of singularities . . . . . . . . . . . . . . . 15 4 Implementation using C++ and Python 19 4.1 Numerical solver in C++ . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.2 Visualization and plotting in python . . . . . . . . . . . . . . . . . . 19 5 Validation 21 5.1 Langmuir oscillations . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 5.2 The RFD pusher . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 6 Comparison of the RFD model with Lorentz Force case 25 6.1 Plasma oscillations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 6.2 Plane wave propagating through plasma . . . . . . . . . . . . . . . . 26 6.3 Plane wave with gaussian perturbation . . . . . . . . . . . . . . . . . 28 6.4 Plane wave with shape factor . . . . . . . . . . . . . . . . . . . . . . 30 6.5 Radiation of a plasma slab . . . . . . . . . . . . . . . . . . . . . . . . 30 6.6 Z-pinch effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 7 Discussion 35 7.1 Behaviour of the model . . . . . . . . . . . . . . . . . . . . . . . . . . 35 ix Contents 7.2 Approximations and numerical sources of error . . . . . . . . . . . . . 36 7.3 Further work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 8 Conclusion 39 Bibliography 41 x List of Figures 2.1 The plane spanned by the vectors E and B is vizualised here. The value of B is fixed to unity while E is determined by the coordinates. The arrows represent the in-plane component of the RFD while the color represents the orthogonal component. Note the discontinuity along the middle as E ·B changes signs. . . . . . . . . . . . . . . . . 4 3.1 An illustration of the simulation cycle in a particle-in-cell code. A single cycle moves the simulation state forwards by dt, and the effect of a different particle pusher is investigated in the results section. . . 12 3.2 An illustration of how the components of the electromagnetic field are arranged in time and space. The components are simulated in this staggered way to facilitate the leap-frogging and thus second-order accuracy of the algorithm. . . . . . . . . . . . . . . . . . . . . . . . . 13 5.1 The normalized total particle momentum as a function of time. A sine is fitted to extract the oscillation frequency. The frequency of the fitted wave ωf = 1.26 · 1014 has a relative error of around 1% compared to the analytical plasma frequency of 1.275 · 1014. . . . . . 21 5.2 The electric field sampled at a point as a function of time. A sine is fitted against the field to extract the oscillation frequency. The fitted frequency of ωf = 1.54 · 1014 has a relative error of 6% compared to the value given by eq 2.24. . . . . . . . . . . . . . . . . . . . . . . . 23 24figure.caption.11 6.1 Momentum in Lorentz force case to the left and RFD case to the right. In the Lorentz case the Langmuir oscillations occur at the correct frequency, while in the RFD simulation the particles move essentially randomly after an extremely short time. . . . . . . . . . . 26 6.2 Plane wave case. Initial condition for the Ey field, which is equal to the Bz field (not shown). This results in a plane wave propagating in the positive x-direction. Also seen is the uniform initial positron distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 xi List of Figures 6.3 Model comparison for the plane wave case. In the top figures the initially uniformly distributed positrons are gathered up by the prop- agating wave. This yields lines of higher positron-density plasma. In the RFD case there are more of these lines. In the bottom fig- ure the trajectories of the particles are compared. The RFD model seems to predict sharp turns that are not present in the Lorentz force simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 6.4 Illustration of the electromagnetic IC for the Gaussian perturbation case. The gaussian is visible in the top figure, showing the Ex field while the Ey field is shown at the bottom. The overlap between the plane wave and gaussian components can be seen in the middle. . . . 28 6.5 Charge and positron densities for the Gaussian case. The upper plots show positron density while the lower ones show the charge density. As the simulation proceeds the particles being pushed against the Gaussian perturbation separates differently charged particles, pro- ducing clusters with higher charge density. There is charge separation in both cases but more strongly in the RFD than in the Lorentz case for the negative charges. The negative charge of the clusters can be more than 10 times greater in the RFD case. . . . . . . . . . . . . . 29 6.6 Charge densities for the Pulse case. The upper and lower plots show densities at earlier and later times respectively. Charge separation can be seen, which is more pronounced in the RFD case . . . . . . . . 30 6.7 Initial condition for the plasma slab case. Two linearly polarized plane waves of the same amplitude propagate into the slab, compressing the plasma between them. . . . . . . . . . . . . . . . . . . . . . . . . . . 31 6.8 Model comparison for the slab case. In both figures we see the particle x position as a function of simulation time. . . . . . . . . . . . . . . 32 6.9 Initial condition for the magnetic field and positron density when the Benett condition is reached. In the simulation using the Lorentz force the particles will compress until reaching an equillibrium state at the density specified by the Benett condition. . . . . . . . . . . . . . . . 33 6.10 The RFD simulation of the z-pinch effect. As Ez increases relative to Bθ the rate of the plasma compression decreases. In stark contrast to the Lorentz case the compression will not stop at any specific density, but continue seemingly indefinitely. . . . . . . . . . . . . . . . . . . 33 xii 1 Introduction Plasma is often an object of interest in physics. Plasma is present in a wide range of circumstances from astrophysical environments to terrestrial laboratories. Often referred to as the fourth state of matter, a plasma is a medium where charged particles interact and give rise to collective behaviour. An example is a gas heated to the point were a considerable number of the atoms are ionized. [1] Another, more exotic example is where the plasma is formed by electrons and positrons, which are charges with equal mass. This type of plasma can be gener- ated by Breit-Wheeler pair production in environments with strong electromagnetic fields, such as those generated by modern high-intensity lasers [8] or near neutron stars [9]. Such strong-field environments also cause moving charges to experience intense radiation losses. This radiation-reaction can be modelled either as a clas- sical force, or semi-classically by probabilistic recoils. The recoil occurs when a photon is emitted, and the occurrence rate is computed within strong-field quantum electrodynamics (SFQED). Simulating plasma dynamics accurately can be very demanding. An effective ap- proach to simulation is to use the kinetic model of plasma. To evaluate the time evolution of the plasma then requires a six-dimensional representation of the parti- cles in phase space, taking into account the particle positions and momenta. Solvers based on this principle, such as particle-in-cell codes, have been very effective his- torically for various studies. This representation is computationally demanding, however, and the demands increase in exotic environments where SFQED effects become more frequent and dramatically restrict the time step of simulation. The Breit-Wheeler pair production is one such effect. Under some circumstances the energy of a photon with an energy above twice the rest energy of the electron can be converted into matter as an electron-positron pair. As an example of the difficulties caused when trying to simulate this, pair production can result in the growth of the number of particles by several orders of magnitude. A common way of taking into account the effect of pair production is to utilize an event generator. The event generator uses an analytical expression to pseudo-randomly generate new pairs appropriately. This accounts for the effect, but takes additional computation time. It also opens up for some surprising worst-case scenarios. A cascade of pairs 1 1. Introduction being produced could end up exhausting the available memory. [10] For many particle-in-cell codes the most computationally intensive part of the sim- ulation will still be the particle mover, however. [5] For the particle dynamics radiation-reaction becomes relevant in strong-field environments. A particle accel- erating in a electromagnetic field will emit radiation. At lower energies, the energy of the photon emitted will be negligible. At higher energies, the energy of these photons can approach the energy of the electrons themselves. This effect must be taken into account, as the emission changes the momentum and trajectory of the electron. It has been shown that in the limit of extremely strong electromagnetic fields the particles tend to travel in the so-called radiation-free direction (RFD). [11] As the particle loses momentum in every other direction but can be accelerated along the RFD, the particles will orient themselves to this direction. In the limit this alignment will be instantaneous. The model is to assume the particles travel entirely or mostly along the RFD. Replacing the mover of a particle-in-cell code with a particle mover built on this model then has the potential to produce results taking into account the effect of radiation reaction. In this thesis a comparison of particle-in-cell simulations based on both this RFD mover and the conventional Boris mover is presented. Results taking into account some quantum effects without using tools such as event generators are presented. Several quantum effects, most notably pair-production, have been left outside the scope of this investigation. The two models predict overall similar behaviour, with some differences in particular cases. 2 2 The Radiation Free Direction and Plasma This section clarifies what is meant by strong fields in the context of this work. The expression used for determining the radiation free direction is the most crucial part of this section. 2.1 Radiation-free direction The approach to plasma simulation described in this work build upon the idea of the radiation free direction (RFD) as a way of describing particle motion in the limit of strong fields. The particles propagating along the RFD yields a complete model of plasma dynamics. This approach was first suggested by A. Gonoskov and M. Marklund in [11]. As the field in the particle’s rest frame becomes stronger radiation reaction starts to play a bigger role in the particle dynamics. The particle traveling in this direction results in the Lorentz force being co-directed with motion resulting in significant suppression of radiation losses. Motion in any other direction will result in a loss of momentum for the particle from radiation reaction. Meanwhile, the field can accelerate the particle in the radiation free direction. This results in particle motion aligning towards the radiation free direction. In the limit of strong fields the strength of the radiation reaction and the accelerating field makes this alignment instant, which is the basis for using the model to describe particle dynamics. The motion is also assumed to be ultra- relativistic. Under these assumptions the particles move with nearly the speed of light independently of the actual momentum gained along the RFD. [11] This direction is described by the expression n±RFD = (√ u− u2(E × B)± [(1− u)|B|E + u(E ·B) B |B| ] ) √ [E2B2 − u(E × B)2] (2.1) 3 2. The Radiation Free Direction and Plasma Figure 2.1: The plane spanned by the vectors E and B is vizualised here. The value of B is fixed to unity while E is determined by the coordinates. The ar- rows represent the in-plane component of the RFD while the color represents the orthogonal component. Note the discontinuity along the middle as E · B changes signs. where the plus and minus signs in the superscript indicates the sign of the particle’s charge. u and w are determined by the expressions u = 2 B2 E2 +B2 1− √ 1− w w , (2.2) w = 4(E × B)2 (E2 +B2)2 (2.3) and E and correspond to the electric and magnetic fields respectively. As the function described in equation 2.1 is a vector field in six dimensions it can be difficult to visualise. It is instructive to imagine the plane spanned by the E and B vectors. The E and B vectors can be normalized so that B has unit length. Varying the size relative to B and direction of E in this plane reveals most of the interesting behaiviour of the RFD. The RFD in this plane is vizualised in figure 2.1. A notable feature is the line discontinuity along the center, where a component of the RFD switches direction. 2.2 Quantum limit It is of some interest to more closely investigate which scales require different models to describe them accurately. To better understand the domain in which the RFD is applicable it is instructive to consider classical synchrotron radiation and when the classical expression starts to fail. A synchrotron is an apparatus in which particles circle in a magnetic field orthogonal 4 2. The Radiation Free Direction and Plasma to their plane of rotation. As the particle rotates, it accelerates, and a particle accelerating in an electromagnetic field emits radiation. As derived in several texts on electromagnetic field theory, e.g [3], the frequency of the radiation typically emitted by a synchrotron is given by ωc = 3qeH 2mec ( γ mec2 )2 (2.4) Where qe and me are the charge and mass of the electron respectively, H is the magnetic field orthogonal to the plane of electron motion, c the speed of light and γ the Lorentz factor. [3] The energy of the electron is Ee = meγc 2 (2.5) and thus the ratio between the typical energy of the synchrotron photon and the electron energy is ~ωc Ee = 3 2 ~qeHγ m4c6 (2.6) where ~ denotes the reduced Planck constant. This ratio can be used as a benchmark for when the classical synchrotron theory needs to be adapted to quantum effects. The ratio increases in proportion to γ and H. After a certain point, the particle is predicted to emit photons with energies greater than its own kinetic energy. When predicted photon energies reaches this point the classical description becomes invalid and a semi-classical description is needed. In the semi-classical description the emitted photon energy distribution is instead computed within SFQED. SFQED instead indicates emission of photons with energies greater than the particle energy is unlikely. Another point is that the electron must experience a recoil from this emission for momentum to be conserved. At low energies the recoil is negligible, but at higher energies the description of the particle dynamics must be reconsidered, as must the photon energies. This recoil is what results in the force of radiation reaction, which gives rise to the particle alignment with RFD. In conclusion: as long as ~ωe Ee << 1, the classical theory should be fine. This is not a restriction for the RFD, as it is still possible for alignment to occur. The main criteria for the validity of the RFD model is that the timescale for alignment to the RFD is small. For estimates on the requirements for rapid alignment in both classical and semi-classical situations, see [11]. 2.3 Strong fields To meaningfully talk about strong fields one needs to compare the intensity with some reference field or quantity. A convenient quantity is the electron mass, yielding 5 2. The Radiation Free Direction and Plasma a characteristic energy, wavelength and charge. When you compare these quantities you get the Schwinger-Sauter field or critical field of quantum electrodynamics Ee λeqe = mc2 ~ mc qe = m2c3 qe~ = Es ≈ 1.323 · 1018 V m . (2.7) As the electrodynamics is a relativistic theory, the comparison needs to be done in an Lorentz-invariant fashion. Two invariants which allow for this comparison of particle and field properties are then χe = γ √ (E + v × B)2 − (E · v/c)2 Es (2.8) χγ = ~ω mec2γ √ (E + (c2k/ω) × B)2 − (E · (ck/ω))2 Es (2.9) Where E and B are the electromagnetic fields, v the velocity of the electron, and ω and k are the frequency and wavevector of the field. χe compares the field strength in the electron rest frame with the typical energy of its emission spectrum. Note the similarity of this comparison of two quantities with the one in the previ- ous section of synchrotron radiation energy with electron energy. In a similar way χe describes the importance of quantum effects from strong electromagnetic fields. Meanwhile, χγ describes the regime where quantum pair production becomes rel- evant. [4] Again, when these ratios are small the classical models for simulating plasma dynamics are useful, and as they approach unity the importance of quantum effects increase. 2.4 Radiation reaction Radiation reaction is the force caused by the emission of photons by an accelerating charged particle. It is this force that gives rise to the particle propagation along the RFD. In the relativistic, yet still classical case the force is approximately given by the expression f clRR = 2 3 q2 em 2 ec ~2 χ2 ev (2.10) which is the dominating term in the Landau-Lifschitz form of the force [4]. This force can be considered as the classical limit of the quantum electrodynamical (QED) effects which give rise to the alignment with the RFD. As the external fields increase in strength and the particles in momentum radiation reaction needs to be considered in a QED framework. Quantum electrodynamics are beyond the scope of this thesis but the result for the average force in the case where χe is large is 6 2. The Radiation Free Direction and Plasma f qRR ≈ −0.37q 2 em 2 ec ~2 χ2/3 e v. (2.11) These expressions can then be used to show that particle momentum aligns with the RFD at a timescale that depends on the field intensities through χe. When this timescale is short enough the RFD model starts to become useful for determining particle motion. As a heuristic the criterion δr = (E/Es)(ω/ωc)−2 > 3 · 108 can be used, as discussed in section 3 of [11], for when the RFD starts to be relevant. As the importance of the RR is not discrete, but a function of the field, we can also consider the empirical results as presented in [11]. These considerations are helpful when evaluating what situations the RFD is accurate for. Note however that the RFD-based model itself is agnostic to the field strengths. 2.5 Plasma properties and waves A plasma is a semi-ionized gas consisting of neutral and charged particles. At larger length scales the plasma is mostly neutral. Any charged region will result in particles of the opposite charge quickly rearranging themselves to shield off the non-neutral region and restore neutrality. The length and time scales for this behaviour are related to the Debye length and the plasma frequency. A plasma can be characterized by its density, electron temperature and has charac- teristic time and space scales given by the plasma frequncy and Debye length. The electron temperature is given by the equipartition theorem and measures the average thermal velocity of the electrons. The Debye length is defined by the equation λ2 D = kBTe n0q2 e (2.12) where kB is the Boltzmann constant and n0 the electron number density. The Debye length is typically the largest length scale at which a small charge perturbation is noticed. Finally, when the plasma is perturbed from its equilibrium state the plasma will oscillate around its equilibrium at the plasma frequency. The plasma frequency ωp is determined by the equation ω2 p = 4πneq2 e me (2.13) There are many possible waves in plasma. One of the arguably simplest waves is the longitudinal electrostatic wave. To examine this wave we follow the approach outlined in [1] and [2]. To start out we assume that the plasma can be modelled as a fluid. This is called the hydrodynamic model of plasma. Assuming the hydrody- namic description of plasma is valid, the one-dimensional plasma-dynamics can be 7 2. The Radiation Free Direction and Plasma described by ∂ne ∂t +∇ · (neve) = 0 (2.14) me [ ∂ve ∂t + (ve · ∇)ve ] = −qeE− ∇pe ne (2.15) pen −γe e = const. (2.16) Further, if the ions present in the plasma are assumed to be stationary relative to the lighter particles, the gradient of the electron pressure can be written as ∇pe = 3pe∇ne. The total solution can also be considered a perturbation against the background solution of ne = n0, ve = v0, E = 0 as ne = n0, ve = v0 +v 1, E = E1 where the subscript 1 indicates the small perturbation. Inserting this into the description above and discarding terms that are second order in the perturbation yields ∂n1 ∂t + n0 ∂v1 ∂dx = 0 (2.17) me ∂v1 ∂t = −qeE− 3kbTE n0 ∂n1 ∂x (2.18) ∂E ∂x = qen1. (2.19) Inserting the monochromatic wave solution X ∼ exp{i(kx+ ωt)} into the above yields the algebraic system −iωn1 + ikn0v1 = 0 (2.20) −imeωv1 = −qeE1 + 3kbTE n0 ikn1 (2.21) ikE1 = −en1 (2.22) and solving this system yields the equation (ω2 − n0e 2 me − 3kbTe me k2)n1 = 0 (2.23) To get any nontrivial solutions from this requires n1 6= 0. This gives us the dispersion relation ω2 = n0e 2 me + 3kbTe me k2 = ωp + 3v2 thk 2 2 (2.24) for longitudinal electrostatic plasma waves. This holds given the assumptions made and additionally as long as k2λd << 1, where λd is the debye length. 8 2. The Radiation Free Direction and Plasma 2.6 The plasma z-pinch By Maxwell’s equations any current will induce a magnetic field. Consider then a cylindrical coordinate system. If a current is assumed to flow in the z-direction the magnetic field will be generated along in the direction of θ. This yields a J × B force which compresses the conductor/plasma, if the current is traveling in a conductor/plasma. In the case of a plasma, there is an equilibrium which is reached when the plasma pressure balances the compressing force. The condition for this to occur is called the Bennet condition [12]. The condition is, as adapted for a electron-positron plasma with stationary ions 4NkBTe = πI2 (2.25) where N is the particle density, assuming equal densities for both particle species, kB is Boltzmann’s constant and I the current traveling along the z-axis. 9 2. The Radiation Free Direction and Plasma 10 3 Structure of a particle-in-cell code The code is a particle-in-cell (PIC) code. For the simulations intended to evaluate the use of the RFD model for high-field or extreme plasma dynamics, the common Boris particle mover is replaced by the RFD particle mover as described in sec- tion 3.5. The mover calculates the RFD function as described in section 2.1 and propagates the particles in the direction at speed c. A point is that while both the implemented particle movers are about as fast per iteration, the RFD mover describes different physics. The RFD mover accounts for the effect of strong radiation reaction and for situations where this effect is relevant the RFD mover should be faster than the alternative ways of simulating these physics. That is, a fair performance comparison should be between the RFD mover and something closer to a SFQED PIC code. 3.1 Particle-in-cell approach to plasma simulation Plasma often consists of a large number of particles. This means some measures are needed to make relevant problems tractable by computer simulation. The particle- in-cell method that has been very successful for this. The particle-in-cell method is based on the idea of simulating a smaller number of macroparticles, each representing 100 − 108 or even more real particles. In a real plasma there are in most cases too many particles to simulate. By using macroparticles the simulator only needs to track perhaps 105 − 106 macroparticles. Care needs to be taken to include enough macroparticles for statistical validity. [5] [6] In a particle-in-cell code, the computation proceeds in four steps for each timestep, as seen in figure 3.1. In order to propagate the simulation state by a single timestep, the code needs to calculate the electromagnetic field on a grid, interpolate the field to particle positions, move the particles, interpolate the current from the particles and repeat. The charge and masses of the particles in the code can be weighted by a factor to simulate higher densities, as the weight cancels out when calculating the dynamics. Again, a sufficient number of particles are needed to faithfully represent the plasma in phase-space. 11 3. Structure of a particle-in-cell code Figure 3.1: An illustration of the simulation cycle in a particle-in-cell code. A single cycle moves the simulation state forwards by dt, and the effect of a different particle pusher is investigated in the results section. 3.2 Advancing the field A relatively simple approach to advancing the field is to use Yee’s algorithm, which works well for solving Maxwells equations. The two equations necessary for the algorithm are ∂H ∂t = −c∇× E (3.1) ∂E ∂t = −c∇×H − 4πJ. (3.2) It is an useful property of the total PIC model that it satisfies the two remaining equations without them being explicitly considered, as long as charge and current are conserved. We wish to simulate a two-dimensional region, so we assume we have a infinite extent in the z-direction with no change in any quantity. This yields all derivatives of z to be 0. Writing out the equation for the Ez time derivative yields ∂Ez ∂t = ∂Hy ∂x − ∂Hx ∂y − 4πJz (3.3) The clever trick used in the Yee algorithm is to stagger the locations at which the E and H fields are stored in time and space. This allows for second order accuracy in time and space when propagating the field. The values of Ez are stored at integer multiples of the space step dx. The values of Ex, Hy, and then Hx,Ey are stored at spatial locations offset by dx/2 or dy/2 respectively, while finally the Hz component is offset by both dx and dy. In addition, the fields are staggered in time. The fields are stored at times staggered by half a time step dt to, again in order to reach second order accuracy. To see the 12 3. Structure of a particle-in-cell code Figure 3.2: An illustration of how the components of the electromagnetic field are arranged in time and space. The components are simulated in this staggered way to facilitate the leap-frogging and thus second-order accuracy of the algorithm. layout of the Yee grid, see figure 3.2. Discretizing the equation yields the update scheme En+1/2 z = En−1/2 z + c dt ( Hynm,k+1 +Hynm,k dx + Hxnm+1,k +Hxnm,k dy ) − dt 4πJz (3.4) for the Ez field, and similar ones for the other fields. In the above the superscribed index denotes the time step, while the sub-indices indicate the location along x- and y-axes respectively. For a detailed and high-quality exposition of this algorithm see e.g. [7]. For this scheme to be stable requires that the Courant condition be fulfilled. Assuming an equal space step dx in the x and y direction the condition is c dt dx < √ 2. (3.5) This is a quite restrictive condition. Another aspect of the algorithm is that it results in numerical dispersion, which can undesirable. Electromagnetic waves with different frequencies will propagate at different rates. This effect can be mitigated by increasing the number of grid cells per wavelength, or by other means, e.g by spectral methods. In this work the effect should not greatly influence the results. 3.3 Interpolation It is necessary to know the field quantities at the particle positions and the current at the grid positions. This can be done by particle in cell interpolation. This works by assigning a contribution from the four closest grid points to the particle position, weighted linearly using how close the particle is to the grid point. In a similar manner a portion of the current produced by the moving particle is assigned to the four closest points. 13 3. Structure of a particle-in-cell code Explicitly this is done by finding the lower left corner of the grid cell the particle occupies. Any grid quantity A (such as E, B or J) is then calculated by using the displacement x and y from the lower left corner to calculate the four weights w00 = (dx− x)(dy − y) (dxdy) (3.6) w10 = x(dy − y) (dxdy) (3.7) w01 = (dx− x)y (dxdy) (3.8) w11 = xy (dxdy) . (3.9) (3.10) Using these equations any particle quantity can be calculated as Aij = wijAp (3.11) where Ap is the corresponding particle quantity. In a similar way any particle quantity can be calculated as Ap = w00A00 + w10A10 + w01A01 + w11A11. (3.12) This interpolation scheme was adapted from [5]. This scheme is a first order method and there are higher order methods for interpolation that can be used. It is impor- tant to note that the interpolation from the grid points to the particles and vice versa must be done in same manner to avoid numerical problems such as particle self-acceleration. For example, using a linear interpolation scheme in one direction and a quadratic one in the other will result in numerical artifacts. 3.4 Moving the particles using the Boris pusher The particles are propagated in time by using the Boris scheme. The Boris scheme is a leapfrog scheme that works by separating the momentum update over a single time step due to the electrical field and treating the effect of the magnetic field as a rotation. A property of the scheme is that it conserves the energy of the particles under the effect of the magnetic field. For a more detailied description of the algorithm and some details on its accuracy see e.g [5]. By introducing the quantity u = γv Newtons equation with the Lorentz force can be written as un+1/2 − un−1/2 dt = qe me [ En + 1 c un+1/2 + un−1/2 2γn × Bn ] . (3.13) 14 3. Structure of a particle-in-cell code Introducing the helper quantities u− and u+ defined by un−1/2 = u− − qeEndt 2me (3.14) un+1/2 = u+ − qeEndt 2me (3.15) and inserting them into 3.13 gives the expression u+ − u− dt = qe 2γnmec (u+ − u−) × Bn. (3.16) To advance u in time is then done by calculating first u− according to the above. The quantities t and s are then calculated as t = qeBdt 2γnmec (3.17) s = 2t 1 + t2 (3.18) and used to calculate the rotation as u′ = u− + u− × t (3.19) u+ = u− + uť × s. (3.20) Adding the final part of the impulse to u+ results in the desired un+1/2 which can be used to move the particles according to xn+1 = xn + un+1/2dt γn+1/2 (3.21) 3.5 RFD pusher and treatment of singularities In contrast to the Boris pusher, the RFD pusher works by simply propagating the particles, which are assumed to be highly relativistic, along the RFD at c as xn+1 = xn + n±RFDc dt. (3.22) When running this pusher there is no need for the velocity to be stored, saving memory. However, as this model is discontinuous, there are some singularities in eq. 2.1. These need to be resolved for use in the numerical scheme. While the analytical expression is well-behaved, inputting the expression as it is into a code will result in division by zero errors. The cases or singularities which result in such errors are dealt with as described in the following paragraphs. 15 3. Structure of a particle-in-cell code The first singularity of importance is in the expression for u when w is close to zero and is dealt with by Taylor expanding around w = 0 as 1− √ 1− w w = 1 2 + w 8 +O(w2). (3.23) Discarding terms of second order and higher yields a good approximation, which resolves the issue. Another problem occurs along the discontinuity where E is orthogonal to B and u approaches unity. As the RFD expression describes a vector field in a six-dimensional space this is difficult to visualise. A way of viewing the function is to fix the B component and to set it to unit length. Taking the plane spanned by E and B and plotting the vector field for many values of E yields figure 2.1. As the more complicated behaviours of the RFD depend on the relative orientation of E and B this view is appropriate. The mentioned discontinuity is then found along the line at the center of 2.1 were E·B changes sign and |E| < |B|. The RFD can be broken into three terms. Consider the first term √ u− u2 (E × B)√ E2B2 − u(E ×B)2 . (3.24) Using relation |E × B| = |E||B| sin θ gives an expression in terms of theta, where theta is the angle between the fields in their common plane. Note that this expression now only gives a magnitude of a vector, with a direction defined by the cross product. This expression is √ u− u2 |E||B| sin θ√ E2B2 − u(|E||B| sin θ)2 . (3.25) If E then approaches the direction orthogonal to B in the plane then θ → π 2 , u(θ)→ 1, sin(θ)→ 1 and the expression becomes undefined. A way to resolve this is to first consider the expression for u u = 2 B2 E2 +B2 1− √ 1− w w = B2(E2 +B2) 2(E × B)2 1− √√√√(1− 4(E × B)2 (E2 +B2)2  (3.26) where the last expression comes from inserting w and simplifying. If θ ≈ π/2 is the angle between E and B then we can define α = π/2 − θ. Taylor expanding then yields E × B = |E||B| cos(α) = |E||B|(1 − α2/2). This approximation then eventually results in 16 3. Structure of a particle-in-cell code B2(E2 +B2) 2E2B2(1− α2/2) ( E2 +B2 − √ E4 +B4 + 2E2B2 − 4E2B2(1− α2/2)2 ) E2 +B2 = 1 2E2(1− α2/2) ( E2 +B2 − √ (B2 − E2) + 2E2B2α2 ) . In the second step the term containing the 4th power of α has been discarded. In the next steps the fact that (B2−E2) > 0 is used to factor out the term from under the root sign. The approximation √ 1 + x ≈ 1 + x/2 can then be used as seen in the following calculations. After this the approximation 1 1−α2 ≈ 1 + α2 is used and several simplifications are made. This results in 1 2E2(1− α2/2) E2 +B2 − (B2 − E2) √√√√1 + 2E2B2α2 (B2 − E2)2  ≈ 1 2E2(1− α2/2) ( E2 +B2 − (B2 − E2)(1 + E2B2α2 (B2 − E2)2 ) ≈ 1 + α2 2E2 (E2 +B2 −B2 + E2 − E2B2α2 (B2 − E2)) = 1 + α2 2E2 (2E2 − E2B2α2 (B2 − E2)) = 1 + α2 − α2 B2 B2 − E2 +O(α4). Define η = E2 B2−E2 > 0. Working with the above expressions and discarding again the higher order term the final expression for u is u = 1− α2η. (3.27) Using this approximation we can write the denominator D of the RFD expression as D ≈ EB|α| √ 1 + η. (3.28) The first term of the RFD can then be written as √ u− u2(E × B) EB|α| √ 1 + η = d̂u √ 1− uEB(1− α2) EB|α| √ 1 + η = d̂u √ η 1 + η (1− α) (3.29) where d̂ points in the direction of the cross product. The second term can approaches zero as α does (1− u)BE EB|α| √ 1 + η = α2η |α| √ 1 + η → 0. (3.30) 17 3. Structure of a particle-in-cell code The final term results in the discontinuity discussed earlier. To show this we use the approximation E ·B = αEB. This yields u(E ·B) EB|α| √ 1 + η = α |α| 1√ 1 + η (3.31) which implies the discontinuity as the fraction α/|α| behaves as the sign function in the expression. Thus all the difficulties in the region where E ·B changes sign have been dealt with. There are some further problematic areas, which can be resolved in the same manner. For the case E = B,E ·B ≈ 0, the RFD proceeds approximately along the E × B direction. For the cases E = 0 and B = 0 the directions used are sign(α)B and E respectively. These directions can be shown to be appropriate using the same methods used earlier in this chapter. 18 4 Implementation using C++ and Python The steps of the simulation procedure described in the previous section must be implemented in a code to produce results. Additionally the data produced must be visualized. For this thesis the numerical solver was implemented in the C++ pro- gramming language, while the visualisation was done using the matplotlib package for the Python programming language. To compile code, handle data, figures and some control logic, bash and GNU make were also used. The code is publicly available at https://github.com/ TornbergNils/rfd_solver. 4.1 Numerical solver in C++ The code written in C++ has four main parts. A solver class containing most of the implemented algorithms, a RFD class which computes the RFD using a given field, an experiment class containing the initial conditions for a simulation and finally the main function. The main function receives arguments and uses these to decide what experiment to run, what model to use and in what directories to save the simulation data. The data is saved into a binary format which can then be loaded into Python. In order to facilitate learning and produce quality code an attempt has been made to use good programming practices such as polymorphism, compartmentalization, clear code and object oriented programming. There is still desirable for improvements in this direction in several areas of the code. Software can always be made better. 4.2 Visualization and plotting in python The python plotting program consists of a file defining several functions which vi- sualise different aspects of the data. As there are three spatial components of each of the fields and the current, and also the positions and velocities for both particle species, there are very many different ways of visualising the data. The main file 19 https://github.com/TornbergNils/rfd_solver https://github.com/TornbergNils/rfd_solver 4. Implementation using C++ and Python can then be used to run any number of selected functions from this file allowing the user to view a desired facet of the data. The focus of the plotting programs are to visualise as many aspects as possible of the results. There are extremely many possible diagnostics and ways of viewing the result, as the problem is relatively high-dimensional. This has as a trade-off increasing the complexity of the produced results and the plotting program, but seems to be worth it or even necessary when evaluating results or debugging. 20 5 Validation 5.1 Langmuir oscillations A critical step in the development of any numerical code is validation. It needs to be shown that the code has no errors in the implementation. For PIC codes a relatively simple test is to excite Langmuir oscillations and to estimate the frequency of the oscillations. A correct oscillation frequency suggest that all parts of the code are correct, that is, the field integration, current and field interpolation and particle mover. For the Langmuir validation the Boris particle mover was used instead of the RFD mover. This ensures the correctness of all parts except the RFD mover. The RFD mover must then be verified separately. Figure 5.1: The normalized total particle momentum as a function of time. A sine is fitted to extract the oscillation frequency. The frequency of the fitted wave ωf = 1.26 ·1014 has a relative error of around 1% compared to the analytical plasma frequency of 1.275 · 1014. The Langmuir oscillations have a frequency ωp determined by the relation ( in CGS 21 5. Validation units ) ω2 p = 4πneq2 e me (5.1) where ne is the electron number density, qe is the electron charge andme the electron rest mass. To excite these oscillations electrons are placed against a positive ion background. They are then set in motion along the x-axis. The electrons are initialized with a thermal velocity vth according to a Gaussian distribution v0 ∼ N(0, vth) and a small sinusoidal perturbation as ve = v0 + v1 = v0 + 0.1 vth sin 2πx λ . (5.2) This results in a current density Jx = nev1qe (5.3) and a field, which in the hydrodynamic approximation has the amplitude Ex = 4πqenev1 ωp . (5.4) Parameter Value ne 5.11e18 vth 1.41e9 ωp 1.275e14 λ 2e− 4 Table 5.1 The values used for the parameters can be seen in table 1. Running the simulation then yields an accurate oscillation of the total particle momentum. Normalized particle momentum versus time is seen in figure 1. The degree of correspondence for the current density and electric field are harder to determine due to noise. The amplitude of the field does however compare favorably to theory. A sine wave fitted to the electric field at over time is seen in figure 3 Extracting the frequencies from these curve fits yield a deviation from theory of less than 1% for the total momentum and less than 6% for the field. There are phys- ical reasons to expect some deviation from theory. The non-zero amplitude of the oscillations and temperature of the plasma should yield slight non-linearity. There is also Landau damping to consider. Thus this comparison is sufficient evidence for the correctness of the code. 22 5. Validation Figure 5.2: The electric field sampled at a point as a function of time. A sine is fitted against the field to extract the oscillation frequency. The fitted frequency of ωf = 1.54 · 1014 has a relative error of 6% compared to the value given by eq 2.24. 5.2 The RFD pusher For the code using the RFD pusher there is no obvious default test case. However, as the code for the two solvers is essentially identical except for the mover, it should be sufficient to benchmark the developed numerical code using independently obtained numerical results. For such a comparison see figure 5.3. As the rest of the solver is valid by the previous section, that is the interpolation steps and FDTD field propagator, the results from the combined solver should then be faithful to the model based on the RFD. 23 5. Validation Figure 5.3: Comparison of numerical RFD in this work versus a similar figure in [11]. As it is impossible to propagate the particles in two directions at once the in- plane component of the numerical RFD is arbitrarily defined to point in the positive direction at the B · E = 0 discontinuity, resulting in the discrepancy between the centers. Parts of the figure are adapted from [11] by A. Gonoskov and M. Marklund under CC BY 4.0 24 6 Comparison of the RFD model with Lorentz Force case The main result of this thesis is the evaluation of the RFD pusher for simulating extreme plasma dynamics. In order to evaluate model performance and behaviour, comparisons are made for several simple test cases with a extensively used and well-tested simulation approach, which is particle-in-cell simulations using the Boris pusher. Simulations using this pusher use the Lorentz force without radiation reac- tion. [3] [6] In this section the two models are compared and evaluated. An important point to consider is that the parameters used in the simulation are scalable. It is the relationship and interplay between the parameters that define the behaviour. For example, if the number and energy of particles, the field strengths and the size of the simulation domain are scaled in a cohesive manner the same effects will appear at many different scales. As the model performance is evaluated here not much attention has been afforded to simulating real scenarios. Many of the scales that yield the same behaviours probably do not occur with any sensible frequency in the real world. If there is an interest in modelling a specific real scenario, more attention must be paid to the specific parameters and not just their relationships. 6.1 Plasma oscillations The first comparison is done for the Langmuir oscillations as described in the vali- dation section. In the Lorentz case the oscillations are excited and are then active until they are dampened out. In the RFD case they will however die out extremely quickly, not completing a single period. As an oscillation depends on momentum to continue, this is actually a reasonable result. There is no notion of momentum in the RFD model. When the particles in the Lorentz case have their maximum momentum, the field is zero. They thus have no potential energy, but only kinetic. They continue along their paths and the energy is converted back to potential energy, restoring the field. In the RFD case, the coherent motion of the particles instead stops when the field is reduced to 25 6. Comparison of the RFD model with Lorentz Force case Figure 6.1: Momentum in Lorentz force case to the left and RFD case to the right. In the Lorentz case the Langmuir oscillations occur at the correct frequency, while in the RFD simulation the particles move essentially randomly after an extremely short time. Figure 6.2: Plane wave case. Initial condition for the Ey field, which is equal to the Bz field (not shown). This results in a plane wave propagating in the positive x-direction. Also seen is the uniform initial positron distribution. zero. Thus the Langmuir oscillations will not occur as they would in the Lorentz case. A conclusion is that the RFD model does not show unstable behaviour such as exponential growth under these conditions. 6.2 Plane wave propagating through plasma The effect of a linearly polarized plane wave in uniform plasma was then inves- tigated. Positrons and electrons were placed uniformly with a temperature of 1 TeV in the simulation region. At this temperature almost every single particle is highly relativistic. Two periods of a plane wave with amplitude 1014 statV/cm were also initialized in the region. The plane wave was set to propagate in the positive x-direction. 26 6. Comparison of the RFD model with Lorentz Force case When simulating the system for the Lorentz and RFD cases the particles tend to bunch up along the maxima and minima of the waves. This yields regions or lines with increased particle density. However as both positrons and electrons bunch up the charge density is still mostly neutral. Figure 6.3: Model comparison for the plane wave case. In the top figures the initially uniformly distributed positrons are gathered up by the propagating wave. This yields lines of higher positron-density plasma. In the RFD case there are more of these lines. In the bottom figure the trajectories of the particles are compared. The RFD model seems to predict sharp turns that are not present in the Lorentz force simulation. The Lorentz simulation predicts the particles mostly traveling in curved lines, while the RFD predicts similar behaviour with some more abrupt turns along the way. A comparison of some sampled particle trajectories can be seen in 6.3. Another difference between the two simulated positron densities it that there seems to be more high-density regions in the RFD case. It is important to note that the behaviour of the Lorentz simulation depends strongly on the temperature of the electrons and the amplitude of the wave. Low wave amplitude relative to the particle temperature yields only thermal behaviour from the particles. Meanwhile there is no notion of temperature in the RFD model, as the particle velocities are always highly relativistic. If the electric fields produced by particle movements are on the same order or greater than the wave amplitude the dynamics are strongly influenced. In this case, many field configurations are dampened out or overwhelmed by noise. If the amplitude of the wave is greater than than the amplitude of the fields generated by particle motion value then changing the amplitude essentially has no effect. 27 6. Comparison of the RFD model with Lorentz Force case Figure 6.4: Illustration of the electromagnetic IC for the Gaussian perturbation case. The gaussian is visible in the top figure, showing the Ex field while the Ey field is shown at the bottom. The overlap between the plane wave and gaussian components can be seen in the middle. Another numerical difference between the two models is that the RFD must be run at a much smaller timesteps to yield robust results. Above a certain timestep size the results of the simulation depend on the timestep size, suggesting numerical instability. Thus the RFD simulation is done with a timestep dtRFD = 1 25dtBoris. A likely reason for this is that the particles are very sensitive to changes in the field in the RFD model. A particle being stepped past a change in the field without changing direction is then a discretization error and should be avoided. 6.3 Plane wave with gaussian perturbation An electrostatic Gaussian perturbation with a peak amplitude of 0.2 times the am- plitude of the plane wave is added to the case in the previous section. The pertur- bation is centered in the simulation region. The initial conditions for the Ex and Ey fields are shown in figure 6.4. This results in the particles being repelled from or attracted to the center depending on the sign of their charge. This sign-dependence results in regions of higher charge as the neutrality of the positron-electron plasma is disturbed. The typical lines of higher particle density seen when a plane wave is present are deformed by the perturbation. The positron density is changed in different ways in the Lorentz and RFD cases. As the particles in the Lorentz case keep their momen- tum, they are in a sense less affected by the perturbation. They do not instantly 28 6. Comparison of the RFD model with Lorentz Force case Figure 6.5: Charge and positron densities for the Gaussian case. The upper plots show positron density while the lower ones show the charge density. As the simulation proceeds the particles being pushed against the Gaussian perturbation separates differently charged particles, producing clusters with higher charge density. There is charge separation in both cases but more strongly in the RFD than in the Lorentz case for the negative charges. The negative charge of the clusters can be more than 10 times greater in the RFD case. turn when the field changes direction as they pass through the perturbation, but continue until their momentum is exhausted. This seems to result in the charge distribution in the Lorentz case being affected across more of the simulation domain, as can be seen in figures 6.5. As the RFD particles turn around instantly as soon as the perturbing field is no longer domi- nant, the effect of the perturbation turns out to be more local than in the Lorentz case. However, this leads to a very strong concentration of RFD electrons at the perturbation while a node of the wave is at the center. During this time the per- turbation seems to dominate the dynamics. Then, as the wave propagates it will carry this concentration of electrons along with it. This seems a likely explana- tion for the behaviour, but further work would be necessary to exclude alternative explanations. The trend seems to be for the RFD the lead to stronger effects locally while the Lorentz model predicts the overall particle distributions to be affected less intensively but more over all. The differences in how the positrons are distributed over time can be seen in figure 6.5 for the Lorentz and the RFD case. This gives a different view of the results, as 29 6. Comparison of the RFD model with Lorentz Force case Figure 6.6: Charge densities for the Pulse case. The upper and lower plots show densities at earlier and later times respectively. Charge separation can be seen, which is more pronounced in the RFD case the strongly negative regions indicate the electrons are clumped together. For the positrons there is similar behaviour. Again the particle peak density is predicted to be higher by the RFD model. 6.4 Plane wave with shape factor The initial conditions of the plane wave case as described previously are used again with a modification. The wave is modulated with a shape function φ(x, y) = e −( x−x0 Lx )2−( y−y0 Ly )2 , where x0, y0, Lx and Ly are all xmax/4 = ymax/4. This results in a "pulse" which propagates in the x-direction. As the pulse is no longer completely symmetrical like the plane wave, charge separation occurs. The electrons and positrons are pushed to the sides due to their differently signed charge. Again, the RFD predicts a stronger effect with more extreme charge densities and higher particle densities. 6.5 Radiation of a plasma slab Particles were placed along the middle of the simulation region as a plasma slab. To the sides of the middle two plane waves were initialized to propagate into the slab. The particles were given a temperature of 1 TeV. This configuration can be seen in figure 6.7 In both the RFD and the Lorentz case the waves propagate into the slab, creating a standing wave in the central region. Initially the waves push the particles along. However, as the wave propagates in the Lorentz case it carries particles along with 30 6. Comparison of the RFD model with Lorentz Force case Figure 6.7: Initial condition for the plasma slab case. Two linearly polarized plane waves of the same amplitude propagate into the slab, compressing the plasma between them. it and out of the central region. The RFD simulation instead predicts the particles being pushed together and oscillating in-place in the standing wave region. The particles oscillate both in the x- and y-directions, the x-component of the oscillation can be seen in figure 6.8. As the wave then propagates away the particle displacement from the center is lower in the RFD case. The Lorentz case predict most of the particles leaving the center of the simulation area while the RFD simulation predicts they stay. Again this can be seen in figure 6.8. 6.6 Z-pinch effect In this scenario a circular region of uniform density electron-positron plasma is initialized in the simulation area. This circular region represents a plasma column. For the Lorentz force simulation the particles are then given a positive or negative velocity along the z-direction depending on the sign of their charge. This gives rise to a current, which is used to calculate the initial magnetic filed in the Bθ direction. These initial conditions give rise to the z-pinch effect, compressing the plasma column described by the circle in 2D geometry. A cutout of the By initial conditions can be seen in 6.9 along with the equilibrium state of the z-pinch as described by the Benett condition. This is to be contrasted with the result of the RFD simulation. In the RFD simulation there needs to be a field in the Ez direction to initiate particle movement. It turns out the plasma behaviour is entirely dependent on the ratio of this Ez field and the strength of the Bθ field. 31 6. Comparison of the RFD model with Lorentz Force case Figure 6.8: Model comparison for the slab case. In both figures we see the particle x position as a function of simulation time. The particles are then instantly set into motion along the z-axis by the Ez field, as described by the RFD model. An effect similar to the z-pinch then occurs. By running the RFD simulation for several different value of the Ez field a qualitative relationship between field strength and model behaviour can be established. For low values of the Ez field, the particles will very quickly compress to a point in the center of the disk. Increasing the value of Ez decreases the rate of this compression. In the limit the particles will move only along the z-axis and not compress at all. This seems to be due to the normalization required for the propagation of the par- ticles by the RFD. The Ez field contribution seems to dominate the magnetic field generating the particle movement. The current, and thus magnetic field generated by the particles is not directly dependent on Ez. This results in the particles moving mostly in the z-direction and ignoring the contribution from the magnetic field. 32 6. Comparison of the RFD model with Lorentz Force case Figure 6.9: Initial condition for the magnetic field and positron density when the Benett condition is reached. In the simulation using the Lorentz force the particles will compress until reaching an equillibrium state at the density specified by the Benett condition. Figure 6.10: The RFD simulation of the z-pinch effect. As Ez increases relative to Bθ the rate of the plasma compression decreases. In stark contrast to the Lorentz case the compression will not stop at any specific density, but continue seemingly indefinitely. 33 6. Comparison of the RFD model with Lorentz Force case 34 7 Discussion 7.1 Behaviour of the model A fundamental assumption of the RFD model is that it relies on extremely high fields for any correspondence of predictions to real world results. However, the model as implemented does not take into account the absolute strength of the fields in any way. The model only considers the electric and magnetic field strengths and orientations relative to each other. In the case where one field component dominates another, the dynamics will be entirely determined by the dominating component. The normalization involved when computing the RFD ensures this. A second property to highlight is that the model has no "memory", as the momentum of the particles is not taken into consideration. This seems to result in many of the divergences from the predictions made using only the Lorentz force. The Langmuir waves being absent from simulations is almost certainly an effect of this. The tendency of the RFD model to predict more extreme charge and particle den- sities could also be a result of this factor. Consider for example the case where a gaussian perturbation is introduced to the plane wave. In the Lorentz simulation the particles are always affected by the perturbation in the sense that it has an effect on their momentum. Meanwhile, the perturbing field has a lesser effect in the RFD model if it is dominated by the wave. When the Gaussian dominates, the parti- cles instead completely discard their previous momentum and change direction. For the electrons, this results in them being strongly concentrated at the center of the Gaussian, which can be thought of representing a collection of positive, stationary ions. As the plane wave propagates and becomes the dominating contribution again, this collection of particles is shoved forward in the direction of the wave. One would expect the particles, which have equal charge, to strongly repel each other. This tendency does not seem to be as strong in the RFD model, however. The strongly negatively charged collection of particles remains. In the pulse case we see a similar tendency of the RFD model to predict higher charge and particle densities than the simulations using the Lorentz force. 35 7. Discussion The other main deviation is for the z-pinch simulation. The Lorentz case has the particle pinch balanced in a sense by the plasma pressure. This pressure is absent during the simulation of this situation using the RFD model. Thus, overall the model behaves similarly to the Lorentz case, but predict higher densities. 7.2 Approximations and numerical sources of er- ror For any simulation approximations must be used. As the results presented in this work are qualitative, it seems unlikely that the mostly small approximations made will fundamentally alter the nature of the results. However this could still happen. For example, while very small timesteps are used, the RFD is calculated by assuming the E and B fields are stored at the same time. Due to the nature of the FDTD scheme this is a simplification that can be bypassed by modifying the numerical scheme. Another problem is the numerical dispersion introduced by the FDTD scheme, alias- ing and spurious heating. As the RFD model is sensitive to noise these are all things that could effect the results of the model. To diminish the likelihood of such effects corrupting the results the simulations have been run at different scales and using different time and spacesteps. As no changes to the character of the results have been detected this indicates numerical stability of the results. 7.3 Further work There are some interesting possibilities for further work. A path is to conjecture a set of well thought out initial conditions, corresponding to some realistic astrophysical scenario. In this scenario the role of radiation reaction would need to be prominent while the effect of other effects (such as pair production) not accounted for negligible. Simulating such a scenario using the RFD model could then provide realistic data for particle behaviour and dynamics. As the movement of the particles under the RFD is due to the emission of radiation, the particle trajectory data could conceivably be used to calculate a characteristic emission spectrum. It would then be possible to use this spectrum to search for signs of the corresponding simulated behaviour in astrophysical or experimental data. Another possibility is to incorporate a version of the RFD mover into a more sophis- ticated solver. This could be used to expand the range of scenarios for which the RFD could be used, and to remove sources of numerical noise, as the RFD mover deployed here seems very sensitive to noise. It is also possible to implement a hybrid model using a weighted average of the Lorentz and RFD velocities. The weighting of the two model velocities could be done as a function of the field strength which 36 7. Discussion would enable the modeling of problems with field strengths of different orders. There is also the approach of incorporating an RFD mover into a more mature framework. There are more sophisticated ways to approach PIC simulations than the one used in this work. Either developing such methods and reducing numerical errors or adding functionality such as a event generator for pair production could prove useful and and allow for new results. 37 7. Discussion 38 8 Conclusion The numerical RFD solver predicts regions with higher densities as compared to the Lorentz force based solver. The removal of momentum from consideration cre- ates some troubles in low-field simulation regions. This is to be expected from a model based on the assumption of extremely strong fields. Most of the troubles of the RFD approach seem to be of this type. To create an entirely useful solver for extreme plasma dynamics, future work should focus on accounting for the assump- tions not solved by the RFD model. These should be low-field situations and pair production. While there seems to be many questions still remaining regarding the behaviour of the model, no extreme numerical problems have been detected in the course of this work. This does not mean that no such problems exist, however. The numerical performance and predictions of the model have been in line with what can be expected when looking at the analytical expression RFD. In conclusion, the RFD model seems to have some promise as a model for plasma dynamics. 39 8. Conclusion 40 Bibliography [1] D. Anderson, M. Lisak, P. Johannisson, M. Marklund, Basic Plasma Physics: Theory and Applications, Gothenburg, Sweden, 2003. [2] F. Chen, Introduction to Plasma Physics and Controlled Fusion, 3rd ed, Los Angeles CA, USA: Springer International, 2015 [3] L. Landau and E. Lifschitz, The classical theory of fields, 4th Rev. English ed., University of Minnesota, USA: Butterworth-Heinemann, 1980. [4] Gonoskov, Arkady Blackburn, Tom Marklund, Mattias Bulanov, S.. (2021). "Charged particle motion and radiation in strong electromagnetic fields", Preprint. [5] C.K. Birdsall, A.B. Langdon, Plasma physics via computer simulation, New York, USA: Taylor and Francis, 2004. [6] J.M. Dawson, "Particle simulation of plasmas", Rev. Mod. Phys., vol 55, nr. 2, ss 403-447, Apr. 1983, doi:10.1103/RevModPhys.55.403. [7] A Taflove, S C. Hagness, Computational Electrodynamics: The Finite- Difference Time-Domain Method, 3rd ed., Norwood, USA: Artech House, 2005. [8] A. R. Bell, J G. Kirk, "Possibility of Prolific Pair Production with High- Power Lasers", Phys. Rev. Lett. vol. no. pp. 101, 200403, Nov. 2008, doi:10.1103/PhysRevLett.101.200403. [9] A. K. Harding and D. Lai, "Physics of strongly magnetized neutron stars", Rep. Prog. Phys. vol 69 no. 9 pp. 2631 Aug. 2006, doi:10.1088/0034-4885/69/9/RO3. [10] A. Gonoskov et. al., "Extended particle-in-cell schemes for physics in ultrastrong laser fields: Review and developments", Phys. Rev. vol. 92, no. 2, pp. 023305, Aug. 2015, doi:10.1103/PhysRevE.92.023305 [11] A. Gonoskov and M. Marklund, Radiation-dominated particle and plasma dy- namics, Physics of Plasmas vol. 25, 093109, Sep. 2018, doi:10.1063/1.5047799 41 Bibliography [12] M. G. Haines, A review of the dense Z-pinch, Plasma Phys. Control. Fusion vol. 53 093001, June. 2011, doi:10.1088/0741-3335/53/9/093001. 42 DEPARTMENT OF PHYSICS CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden www.chalmers.se www.chalmers.se List of Figures Introduction The Radiation Free Direction and Plasma Radiation-free direction Quantum limit Strong fields Radiation reaction Plasma properties and waves The plasma z-pinch Structure of a particle-in-cell code Particle-in-cell approach to plasma simulation Advancing the field Interpolation Moving the particles using the Boris pusher RFD pusher and treatment of singularities Implementation using C++ and Python Numerical solver in C++ Visualization and plotting in python Validation Langmuir oscillations The RFD pusher Comparison of the RFD model with Lorentz Force case Plasma oscillations Plane wave propagating through plasma Plane wave with gaussian perturbation Plane wave with shape factor Radiation of a plasma slab Z-pinch effect Discussion Behaviour of the model Approximations and numerical sources of error Further work Conclusion Bibliography