MODELING AND SIMULATION OF STRONGLY COUPLED PLASMAS   By  Rahnuma Rifat Chowdhury               A THESIS  Submitted to Michigan State University in partial fulfillment of the requirements for the degree of  Electrical Engineering- Master of Science   2016   ABSTRACT MODELING AND SIMULATION OF STRONGLY COUPLED PLASMAS By Rahnuma Rifat Chowdhury The objective of this work is to develop new modeling and simulation tools for studying strongly coupled plasmas (SCP). Strongly coupled plasmas are different from traditional plasmas as potential energy is larger than the kinetic energy. The standard plasma model does not account for some major effects in SCP:  1) the change in the permittivity  2) the impact on relaxation of the charged particles undergoing Coulomb collisions in a system with weakly shielded long range interactions 3) the impact of statistical fluctuations in strongly coupled plasmas that leads to non-Markovian effects.  Proper modeling of such systems through consideration of Lévy flight processes gives rise to fractional derivatives in time that result in an incorporation of time history in the model. A Lévy flight is a random walk in which the steps are defined in terms of the step-lengths, which have a certain probability distribution, with the directions of the steps being isotropic and random. Lévy processes in the plasma give rise to fluctuations in medium through which the electromagnetic waves are propagating. Averaging over the Lévy processes will allow us to relate to other important parameters in the plasma.  iii               This thesis is dedicated to my parents, who have never stopped believing in me. iv  ACKNOWLEDGMENTS      I would first like to thank my thesis advisor Professor John Verboncoeur of the Department of Electrical and Computer Engineering at Michigan State University. He consistently steered me in the right the direction whenever he thought I needed it.  I would also like to acknowledge Professor Andrew J. Christlieb of the Department of Mathematics at Michigan State University as the co-advisor of this thesis. The door to Prof. Christlieb office was always open whenever I ran into a trouble spot or had a question about my research or writing. In addition, I would like to thank Assistant Professor Mohsen Zayernouri of the Department of Mechanical Engineering at Michigan State University as the committee member of this thesis. I am gratefully indebted to him for his very valuable comments on this thesis. Moreover, I would like to express my gratitude to Michael Murillo, Gautham Dharuman, Guy Parsey and Janez Krek. Without their passionate participation and input, this work could not have been successfully conducted. Finally, I must express my very profound gratitude to my parents and to my fiancé for providing me with unfailing support and continuous encouragement throughout my years of study and through the process of researching and writing this thesis. This accomplishment would not have been possible without them. Thank you. v  TABLE OF CONTENTS  LIST OF TABLES ......................................................................................................................vi  LIST OF FIGURES ...................................................................................................................vii ............viii   INTRODUCTION ........................................................................................................................1    Work Plan....................................................................................................3   MOLECULAR DYNAMICS.........................................................................................................5   THE PARTICLE IN CELL (PIC) METHOD ...............................................................................8   The PIC Algorithm .............................................................................................................8   THE PARTICLE MESH (PM) METHOD...................................................................................10   ................................................................................................11   THE PARTICLE-PARTICLE-PARTICLE MESH (P3M) METHOD .......................................13   .............................................................................................13   .............................................................................................14   METHOD OF CELL LINKED LISTS.........................................................................................15   ................................17   Steps for Executing the Code 17   Simulation Results for Non-Collisional Cases..18   Simulation Results for 24  BIBLIOGRAPHY .......................................................................................................................28vi  LIST OF TABLES  Table 7.1     Simulation Results for Both Non-Collisional (Case 1-4) and Collisional (Case 5-9)                            vii  LIST OF FIGURES      Figure 6.1      The division of a region of the simulation space into cells [8]...........15 Figure 7.1      Energy after the equilN=1000..............18 Figure 7.2      Temperature after the equilibration ...........19 Figure 7.3      ...19 Figure 7.4      Individual particle differential velocities , N=100020    Figure 7.5         Figure 7.6      Distribution function of individual particle velocities in log-log scale for                        dt.21  Figure 7.7      Data fitting for mean square differential velocities ,   Figure 7.8      Data fitting for mean square differential velocities in log-log scale                         22  Figure 7.9      .....22  Figure 7.10    Temperature after the equilibration p01, N=1000....23  Figure 7.11    ..24  Figure 7.12    Energy after the equ, N=1000..24  Figure 7.13    Temperature after the equi, ..25  Figure 7.14    Energy after the equilibration p01, N=1000....25  Figure 7.15    Temperature after the equilibration p0126  Figure 7.16    Energy after the equilibration p001, ...26  Figure 7.17    Temperature after the equilibration p001..27      viii  KEY TO SYMBOLS  Ti Ion temperature a Wigner-Seitz radius Z Degree of ionization  Coulomb coupling parameter r Radial distance to the particle  Inverse screening parameter m Mass of the affected particle D Debye length ne Electron temperature ni Ion density Te Electron temperature e Electron charge N Number of particles Neq Number of steps in the equilibration phase Npd Number of steps in the production phase dt Time step rcut Cut-off radius L Length of the simulation box wpi Ion plasma period mi Ion mass  1  INTRODUCTION  Lévy diffusion is a diffusion process with a non-linear relationship to time, in contrast to a typical diffusion process, in which the mean r2, of a particle is a linear function of time [1]. It can be described by a power law,   r2 ~ Dt (1) Where, D is the diffusion coefficient and t is the elapsed time. In a typical diffusion -icle undergoes sub-diffusion. For ultra-cold plasmas, with temperatures from 1K to 1000K, it is possible to capture the two body auto-correlation function for a low temperature ion-electron plasma through the use of a highly resolved Particle-In-Cell (PIC) model. Because one of the dominant interactions in ultra-cold SCP is the Coulomb force, a resolved PIC calculation with an average of one particle per cell can simulate relaxation in ultra-cold plasmas. But with moderate increases in density or length scales, resolved PIC calculations are not practical. This has inspired us to consider the extension of classical plasma models to study ultra-cold SCP, which will then be used as a benchmark tool for the development of models to describe the physics of SCP where it is not possible to provide an ultra-resolved particle based calculation. The key objective is to develop new plasma models that treat long range correlation as a subscale feature [2]. The unique aspect of SCP is that the potential energy exceeds the kinetic energy. Strong coupling is defined in terms of the dimensionless parameter, often referred to as the Coulomb coupling parameter, .                                                        2/kTia      (2) 2  Eq. (2) describes the ratio of the potential energy to kinetic energy [3]. In Eq. (2): Z = the charge number of the gas species  e = the unit charge k   Ti= the temperature of the ion species in Kelvin  a = [3/(4n)]1/3 = the Wigner-Seitz radius (mean inter-particle distance)  n = the density of the species << 1, the charged species in the plasma has no long range correlation and binary collisions characterize Coulomb scattering for that species. 1, the plasma species in question begins to exhibit long range correlaincreases in these systems, the plasma exhibits a collective behavior, giving the system properties resembling liquids and solids [3, 4, 5]. Plasmas with strong coupling have been created and studied using a range of experimental methods. But there are issues with studying SCP in the context of each of these experimental setups. In many of these systems, the ability to accurately determine the initial condition as well as take accurate non-invasive measurements is difficult. Furthermore, experimental methods do not permit simultaneous measurement of detailed phase space quantities, which can be critical to understanding wave-plasma interaction, collision dynamics, and even the details of anisotropic transport. In particle and atomic physics, a Yukawa potential (also called a screened Coulomb potential) is a potential of the form  UYukawa(r)= - 2 [exp(-] (3) the Coulomb coupling parameter which defines how strongly coupled the 3  system is.  As the system temperature decreases, the system becomes more and more strongly coupled. another scaling constant known as the inverse screening parameter. depends on Debye length, D as follows:   D (4) Debye length, D is the measure of a charge carrier's net electrostatic effect in solution, and how far those electrostatic effects persist. A Debye sphere is a volume whose radius is the Debye length. With each Debye length, charges are increasingly electrically screened. With every Debyelength, the electric potential due to a given charge will decrease by 1/e. Debye length is defined as follows:  D = (Te ee2)1/2     (5)             where, ne = Zni ne = electron temperature ni = ion density Te = electron temperature Hence, from ion density and ion temperature , and from electron temperature and ion  Work Plan  We are trying to look at particular Ultra Cold (UC) experiment with our existing set of simulation toolsDepending on that, we are going to validate the existing models or move towards developing a Fractional Calculus Model if there is any e 1. The 4  hypothesis is that jumps in v due to collisions plus long range interactions might produce Levy -simulations for both collisional and non-collisional cases and do least square fitting to the tail of the differential velocity distribution.  Levy Flight is a random walk in which step lengths have a probability distribution that is heavy tailed:  f(x)  x (6) If distribution of mean-2 has a fat tail for some effective Coupling, then Levy-flight like behavior is possible. Therefore, Levy flights might be possible in a non-equilibrium plasma system [6].   5  MOLECULAR DYNAMICS Popular molecular simulation techniques such as molecular dynamics or Monte Carlo are used to study the physical and chemical processes occurring in systems containing large numbers of particles. These methods require evaluation of either the total potential energy of a system of N particles (VTot) or the gradients of the potential energy. The total potential energy consists of terms that describe the various interactions among the particles in the system. These terms are usually functions of internal coordinates, such as internuclear distances between two particles, bond angles among three particles, or torsional angles among four particles. For condensed phase modeling, the total potential energy is often described as a sum of two-body interactions over all particle pairs.  The interaction terms are typically simple functions of the internuclear distance rij between particles i and j.   (7) The evaluation of Eq. (7) and the gradients are usually the most computationally demanding steps in a simulation, even if the functional forms for V(rij) are extremely simple. Brute force evaluation of Eq. (7distances. In a molecular dynamics simulation, each integration step often requires the evaluation of Eq. (7) and its gradients more than once depending on the integration scheme that is chosen [7]. It is clear that methods to reduce the computational burdens associated with numerous evaluations of Eq. (7) are required. The most obvious recent approaches are to modify the codes for scalable platforms. However, modifications of existing algorithms designed to reduce the computational burdens associated with evaluation of Eq. (7) can be made to increase the serial 6  performance and exploit scalable architectures to achieve enhanced performance. The algorithms that were developed can reduce potentially unnecessary computations of the internuclear distances for particle pairs used in the evaluation of Eq. (7). Common strategies to reduce the computational demands associated with Eq. (7) include the use of simple functions to describe the pair interaction potentials, and the assumption that the interaction between two particles is negligible beyond a certain cutoff distance rcut. The assumption of a cutoff distance in the interaction potential allows for a reduction in computational time, since the interaction between particles separated by distances exceeding rcut are not calculated. The easiest and most direct way to determine the set of internuclear distances that are within rcut is to evaluate all distances between all pairs, and eliminate those that exceed rcut. This step requires a potentially large number of unnecessary calculations, and might be the costliest computational step in such a simulation. A reduction of unnecessary calculations of internuclear distances can be accomplished through the use of the Verlet neighbor list [8]. This method requires the construction of a list of neighbors for each parrticle. A particlely defined to be all of the particles that are within a distance slightly greater than the range of the interaction potential. Information about the neighbors is stored in arrays. For the duration of the simulation or until the lists are updated, each particle is assumed to interact only with the particles on its neighbor list. The internuclear distances, interaction potentials, and forces are evaluated for each particle and its neighbors only. The list may be periodically updated to allow for the movement of particles into or out of the interaction range. Brute force construction or update of the list requires the when the system contains a relatively small number of particles [8, 9]. However, as the system 7  becomes larger, the memory requirements for maintaining the neighbor lists become prohibitive. Alternative methods for the efficient determination of the interacting neighbors for each particle include grid or cell approaches [8, 10-12]. These approaches partition the simulation space into grids or cells, to which the particles are assigned by virtue of their positions relative to the cells. Since each cell has an unchanging set of neighboring cells that contain the volume within the distance rcut of that cell, a particle associated with one of the cells has as its neighbors those particles assigned to the same or neighboring cells. The implementations of these methods usually assign the particles to the cells at each integration step. However, the same considerations used for the frequency of updating the Verlet neighbor-lists are applicable here. There is some overhead associated with these methods, and they are preferable only for systems that contain more than 1000 particles [8]. These methods substantially reduce the number of unnecessary internuclear distance calculations in evaluating Eq. (1), but do not completely eliminate unnecessary computations.   8  THE PARTICLE IN CELL (PIC) METHOD  The particle-in-cell (PIC) method refers to a technique used to solve a certain class of partial differential equations [13-15]. In this method, individual particles in a Lagrangian frame are tracked in continuous phase space, whereas moments of the distribution such as densities and currents are computed simultaneously on stationary mesh points. The method typically includes the following procedures:  Integration of the equations of motion.  Interpolation of charge and current source terms to the field mesh.  Computation of the fields on mesh points.  Interpolation of the fields from the mesh to the particle locations. Models which include interactions of particles only through the average fields are called PM (particle-mesh). Those which include direct binary interactions are PP (particle-particle). Models with both types of interactions are called PP-PM or P3M. The PIC Algorithm i. The plasma is described by a number of computational particles having position xp, velocity vp and each representing a fixed number Np of physical particles. ii. The equations of motion for unmagnetized particles are advanced by one time step using a leap frog integrator,   (8)   (9)             using the particle electric field from the previous time step. iii. The charge densities are computed in each cell using: 9    (10)  iv. The Poisson equation is solved:   (11)             and the electric field Ei in each cell is computed:   (12) v. From the field known in the cells, the field acting on the particles is computed as:   (13)             which is used in the next cycle. vi. The cycle restarts [16].     10  THE PARTICLE MESH (PM) METHOD  Particle Mesh (PM) is a computational method for determining the forces in a system of particles [17]. These particles could be atoms, stars, or fluid components and so the method is applicable to many fields, including molecular dynamics and astrophysics. The basic principle is that a system of particles is converted into a grid (or "mesh") of density values. The potential is then solved for this density grid, and forces are applied to each particle based on what cell it is in, and where in the cell it lies. Various methods for converting a system of particles into a grid of densities exist. Once the density distribution is found, the potential energy of each point in the mesh can be determined from the differential form of Gauss's law, which after identifying the electric field E as the negative gradient of the electric gives rise to a Poisson equation that is easily solved after applying the Fourier transform. Thus it is faster to do a PM calculation than to simply add up all the interactions on a particle due to all other particles for two reasons: firstly, there are usually fewer local grid points than particles, so the number of interactions to calculate is smaller, and secondly the grid technique permits the use of Fourier transform techniques to evaluate the potential, and these can be very fast. PM is considered an obsolete method as it does not model close interaction between particles well. It has been supplanted by the Particle-Particle Particle-Mesh method, which uses a straight particle-particle sum between nearby particles in addition to the PM calculation.   11  The PM Algorithm Huge model systems with 1/r potentials (celestial masses or microscopic charged 4-108 ions, electrons, or stars. Let us consider a square cell of Mshould contain an average of 10-100 particles. The equation of motion for the particle k is:    (14) r is the electrostatic or gravitational potential. It is determined by charge (or  (r, t) defined by the positions of all superparticles. r at time tn at the centers of the subcells, the given configuration of superions is first replaced by lattice-i,j. The easiest discretization method is the nearest grid point (NGP) rule:    (15) The next step is the calculation of the values of the potential produced by the charge lattice at all cell centers. The field strength at the position rk of some superparticle k in cell (i, j) is then:   (16)   (17) Next we integrate the equation of motion:    (18)   (19) 12  which completes the time step.  Combining the PM method and the molecular dynamics technique, we may take into account the short-range forces up to a certain interparticle distance, while the long-range contributions are included by the particle-mesh procedure. This combination of particle-particle and particle-mesh methods has come to be called the PPPM or P3M technique [17].   13  THE PARTICLE-PARTICLE-PARTICLE MESH (P3M) METHOD  ParticleParticle-ParticleMesh (P3M) is a Fourier-based Ewald summation method to calculate potentials in N-body simulations [18-21]. The potential could be the electrostatic potential among N point charges i.e. molecular dynamics, the gravitational potential among N gas particles in e.g. smoothed particle hydrodynamics, or any other useful function. It is based on the particle mesh method, where particles are interpolated onto a grid, and the potential is solved for this grid (e.g. by solving the discrete Poisson equation). This interpolation introduces errors in the force calculation, particularly for particles that are close together. Essentially, the particles are forced to have a lower spatial resolution during the force calculation. The P3M algorithm attempts to remedy this by calculating the potential through a direct sum for particles that are close, and through the particle mesh method for particles that are separated by some distance. The P3M Algorithm i. Mapping the particle xi to the mesh ii. Calculating   (20) iii. Subtracting   (21) from  iv. Computing Emesh and interpolating Emesh to particles v. Calculating the Ewald sum over particles   (22) 14    (23) Ewald Sum Ewald summation is a method for computing long-range interactions (e.g., Coulombic interactions) in periodic systems. It was first developed as the method for calculating electrostatic energies of ionic crystals, and is now commonly used for calculating long-range interactions in computational chemistry [22]. Ewald summation is a special case of the Poisson summation formula, replacing the summation of interaction energies in real space with an equivalent summation in Fourier space. In this method, the long-range interaction is divided into two parts: a short-range contribution, and a long-range contribution which does not have a singularity. The short-range contribution is calculated in real space, whereas the long-range contribution is calculated using a Fourier transform. The advantage of this method is the rapid convergence of the energy compared with that of a direct summation. This means that the method has high accuracy and reasonable speed when computing long-range interactions, and it is thus the standard method for calculating long-range interactions in periodic systems. The method requires charge neutrality of the molecular system in order to calculate accurately the total Coulombic interaction.          15  METHOD OF CELL LINKED LISTS  In the conventional method of cell-linked lists, the simulation space is partitioned into cells, the edges of which are no smaller than cutoff distance of the interaction potential. The particles are then assigned to the various cells, by virtue of their position in the simulation space. A linked-list of the particle indices is created during the sorting procedure. Also, at the beginning of a simulation, an array that contains a list of cell neighbors for each cell is created. The list remains fixed unless the simulation space changes during the simulation. This is dynamically adjusted to balance the particle numbers. A cell icell has as its neighbors any cell that contains at least one point that is within the distance rcut of any point within icell. Since the conventional method requires that the edges of each cell be no smaller than rcut, each cell has eight nearest neighbors (we are assuming periodic boundary conditions in all dimensions). These requirements ensure that all particles that are within the interaction range of any particle within icell are assigned to the eight nearest-neighbor cells of icell or icell itself. All particles occupying cells other than these are outside the interaction range of any particle located within icell [8].  Figure 6.1 The division of a region of the simulation space into cells [8] In Figure 6.1, both the x and y y hereafter) equal rcut. Evaluation of Eq. (7) occurs through looping over the cells using the linked-list of particles rather than accessing the particle indices sequentially as written in Eq. (7). This method 16  dramatically reduces the number of unnecessary internuclear distance calculations that would modifications can be made to further reduce the number of unnecessary distance calculations. In the conventional method, the distances between all particle pairs located within the rectangular arey are calculated. y = rcut, the area within which all distances are calculated is 9rcut2. The area within the cutoff radius for a single particle is only cut2. Thus, the traditional cell-linked list method calculates distances between all particle pairs within an area that is almost three times larger (or more, since cut, i = x or y) than that actually required for a particle [8].   17  RESULTS  Steps for Executing the Code   were set to the desired values. This defines the physical parameters.  Number of particles, N, in the input file yukawa.in was changed.  Cut-off radius, rcut was set following the mandatory relation to be satisfied: rcut < L/2. Example: rcut = L/5, L/10 etc. Interactions for r >rcut  are ignored.  According to the number of particles, the cut-off radius should be changed. The relationship between the number of particles, N and the length of the simulation box which is linked to rcut is as follows:  ni = N/L3 (24)  niaw3 = N/ ()3 (25)  = L/ aw   1/3 (26)  Yukawa_MD.py file was run to generate the output files i.e. irp_eq.out and irp_pd.out (saves particle positions, velocities and accelerations during the equilibration and production phase respectively over the simulation time), temp_eq.out and temp_pd.out (saves the temperature of the system during the equilibration and production phase respectively), energy_eq.out and energy_pd.out (saves the total energy, kinetic energy and potential energy of the system during the equilibration and production phase respectively vs. time).  The energy_plot.py was run to generate the plot of the total energy, kinetic energy and potential energy of the system vs. time.  The temperature_plot.py was run to generate the plot of the temperature of the system 18  during the equilibration and production phase.  The velocity_distribution.py was run to generate the distribution function of the particle velocities during the production phase.  The diff_velocity_distribution.py was run to generate the distribution function of the particle differential velocities and also to do the least square fitting that calculates the co-efficient a and b values.  All the distances are in the units of aw.  Times are in units of inverse ion plasma period wpi-1 where  wpi ie2/mi)1/2 (27)  All the energies and temperatures are in units of e2/aw.  Simulation Results for Non-Collisional Cases  The system was run for  1000 particles.    Figure 7.1 Energy after the equilibration p, N=1000 19   Figure 7.2 Temperature after the equilibration phase for , N=1000  We can see from Figure 7.1 and Figure 7.2 that the energy and the temperature of the system do not stabilize which means that the system has not been equilibrated. Also the kinetic energy exceeds the potential energy, which does not imply SCP.    Figure 7.3 Individual particle velocities , N=1000  Figure 7.3 represents all the particle velocities during the production phase. We can see jumps in particle velocities due to collisions and long range interactions.  20   Figure 7.4 Individual particle differential velocities , N=1000 Figure 7.4 represents all the particle differential velocities during the production phase. We can also see jumps in particle velocities due to collisions and long range interactions.   Figure 7.5 Distribution function of individual particle velocities , N=1000  21    Figure 7.6 Distribution function of individual particle velocities in log-log scale , N=1000  The distribution function was calculated by taking the average of all the particle velocities over all the time steps in the production phase. Jumps in velocities are more prominent in Figure 7.5 and Figure 7.6 which plots the distribution function of individual particle velocities during the production phase in linear scale and log-log scale respectively.  Figure 7.7 Data fitting for mean square differential velocities , N=1000 22    Figure 7.8 Data fitting for mean square differential velocities in log-log scale , N=1000  The mean square velocity change was calculated by taking the average of the differential velocities of all the particles over all the time steps during the production phase. Figure 7.7 represents only the tail of the differential velocity distribution and Figure 7.8 is a log-log representation of Figure 7.7. We can see from these figures that the distribution has a fat tail. After least square fitting to the differential velocity distribution, the value of the co-efficient b comes out as 1.46 which implies Lévy type behavior.  Figure 7.9 Energy after the equilibration p01, N=1000 23  Next, the system was run for 1000 particles with  but with reduced time step, dt=0.01. From Figure 7.9 and Figure 7.10, we can see that the energy and the temperature has been stabilized and the system has been equilibrated properly at much lower energy levels.   Figure 7.10 Temperature after the equilibration p01, N=1000     Figure 7.11 Individual particle differential velocities 01, N=1000  Figure 7.11 represents all the differential particle velocities over all the time during the 24  production phase. No large jumps in Figure 7.11 means the system has been stabilized properly. The value of the co-efficient b after least square fitting comes out as 5.49 which implies no Lévy type behavior. Simulation Results for Collisional Cases  The system was run for 1000 particles, dt=0.1 with collisions added into the system. At each time step, we took 5% of the particles and rotated the direction of the velocity by a random amount while keeping the magnitude the same. We did this rotation of velocity V by,  Vx= Csin(u)cos(t) (28)  Vy= Csin(u)sin(t) (29)  Vz= Ccos(u) (30)  Where, C= [Vx_old2+Vy_old2+Vz_old2]1/2 (31) We can see from Figure 7.12 and Figure 7.13 that the energy and the temperature have not been stabilized and the system has not been equilibrated properly.   Figure 7.12 Energy after the equilibration p, N=1000  25    Figure 7.13 Temperature after the equilibration p, N=1000  After tail fitting to the differential velocity distribution, the value of the co-efficient b comes out as 0.33 which implies Lévy type behavior. Next, the system was run for 01 with collisions included. We can see from Figure 7.14 and Figure 7.15 that the energy and temperature of the system have been stabilized and the system has been equilibrated properly.   Figure 7.14 Energy after the equilibration p01, N=1000 26    Figure 7.15 Temperature after the equilibration phase 01, N=1000  The value of the co-efficient b after least square fitting comes out as 0.65 which implies Lévy type behavior. Next, the system was run for 1 with collisions included. We can see from Figure 7.16 and Figure 7.17 that the energy and temperature of the system have been stabilized and the system has been equilibrated properly.   Figure 7.16 Energy after the equilibration p001, N=1000 27   Figure 7.17 Temperature after the equilibration p001, N=1000 The value of the co-efficient b after least square fitting comes out as 0.97 which implies Lévy type behavior. All the simulation results have been summarized in table 7.1. Table 7.1 Simulation Results for Both Non-Collisional (Case 1-4) and Collisional (Case 5-9) Cases  Case N   dt Neq Npd rcut b Lévy 1 1000 1 1 0.1 4000 2000 4 1.46 Yes 2 1000 1 1 0.01 4000 2000 4 2.69 No 3 1000 1 1 0.01 40000 20000 4 3.49 No 4 1000 10 1 0.1 4000 2000 4 3.96 No 5 1000 1 1 0.1 4000 2000 4 0.33 Yes 6 1000 1 1 0.01 40000 20000 4 0.65 Yes 7 1000 1 1 0.001 400000 200000 4 0.97 Yes 8 1000 10 1 0.1 4000 2000 4 1.39 Yes 9 1000 10 1 0.01 40000 20000 4 1.81 Yes    28                BIBLIOGRAPHY   29  BIBLIOGRAPHY  1. Liu, B., Goree, J. and Feng, Y. (2008). Non-Gaussian statistics and superdiffusion in a driven-dissipative dusty plasma.  2. Christlieb, A.J., Cheng, Y., Causely, M. and Verboncoeur, J. Modeling and Simulations of Strongly Coupled Plasmas.  3. Brush, G., Sahlin, H.L. and Teller, E. (1966). Monte carlo study of a one-component plasma.  The Journal of Chemical Physics 45: 2102.  4. Ichimaru, S. (1982). Strongly coupled plasmas: high-density classical plasmas and degenerate electron liquids. Reviews of Modern Physics 54 (4): 1017.  5. Morfill, G.E, Thomas, H.M, Konopka, U. and Zuzic, M. (1999) The plasma condensation: Liquid and crystalline plasmas. Physics of Plasmas 6: 17691780.  6. Christlieb, A.J., Choi, Y., Dharuman, G., Jain, M., Verboncoeur, J., Murrilo, M. and Roberts, J. (2016). Understanding Ultra Cold Plasmas Through Simulations [Powerpoint Slides].    7. Allen, M.P.  and Tildesley, D.J. (1990).  Computer Simulation of Liquids. NY: Oxford University Press.  8. Mattson, W. and Rice, B.M. (1999). Near-neighbor calculations using a modified cell-linked list method. Computer Physics Communications 119: 135-148.  9. Morales, J.J., Rull, L.F. and Toxvaerd, S. (1989). Efficiency test of the traditonal MD and the link-cell methods.  Computer Physics Communication 56(2): 129-134.  10. grid. Journal of Computational Physics 66: 1-20.  11. Lambrakos, S.G. and Boris, J.P. (1987). Geometric properties of the monotonic lagrangian grid algorithm for near neighbor calculations. Journal of Computational Physics 73: 183-202.  12. Brugé, F. (1993). Systolic calculation of pair interactions using the cell linked-lists method on multi-processor systems. Journal of Computational Physics 104: 263-266.  13. Birdsall, C. and Langdon, A. (1985). Plasma Physics via Computer Simulation. NY: McGraw-Hill.  14. Hockney, R. W. and Eastwood, J. W. (1988). Computer Simulation Using Particles. PA: Taylor and Francis.   30  15. Verboncoeur, J. (2005). Particle simulation of plasmas: Review and advances. Plasma Physics Controlled Fusion 47: A231A260.  16. Lapenta, G. Mathematical Derivation of the PIC method. pp. 21. in Particle In Cell Method, A brief description of the PIC Method. Centrum voor Plasma Astrofysica, Katholieke Universiteit Leuven.  17. -http://homepage.univie.ac.at/franz.vesely/simsp/dx/node48.html.  18. Toukmaji, A.J. and Board Jr., A.J. (1996). Ewald summation techniques in perspective: a survey. Computer Physics Communications 95(2-3): 73-92.   19. Eastwood, J.W., Hockney, R.W. and Lawrence, D.N. (1984). P3M3DP-the three-dimensional periodic particle-particle/particle-mesh program. Computer Physics Communications 35: C618-C619.  20. Deserno, M. and Holm, C. (1998). How to mesh up Ewald sums. II. An accurate error estimate for the particleparticleparticle-mesh algorithm. The Journal of Chemical Physics 109: 7694.  21. -body simulations (gravitatiohttp://www.scholarpedia.org/article/N-body_simulations#P3M_and_PM-Tree_codes.