Molecular Dynamics Simulations of Complex (Dusty) Plasmas

Complex or dusty plasmas are multi-component plasmas, which, in addition to the usual plasma components i.e. ions and electrons, contain micron-sized particles, also called grains or dust. These particles acquire high electric charges, and interact collectively over long distances. Like colloids, complex plasmas form solidand liquid-like structures with long range correlations and exhibit phase transitions. Unlike colloids, they exhibit a range of dynamic phenomena such as particle-mediated linear and nonlinear waves, shocks, wakes, instabilities, etc. Complex plasma crystallisation has been theoretically predicted (Ikezi, 1986) and subsequently discovered in the early 1990s (Chu & I, 1994; Hayashi & Tachibana, 1994; Melzer et al., 1994; Thomas et al., 1994) giving rise to the interest in the field among the whole physics community.

their experimental verification techniques will be detailed in Sections 3 and 4 respectively. Finally, numerical results and their comparison with experiments will be discussed in Section 5 and concluding remarks will be presented in Section 6.

Overview of complex (dusty) plasmas
The term "dusty plasmas" originated in astrophysics, where polydisperse (various size and shape) dust grains are exposed to various charged particles, ionised gases, and ionising radiation. Laboratory studies often involve high quality monodisperse (same size) microspheres added to gas discharges to study various complex and collective phenomena that do not occur in natural dusty plasmas. Thus the term "complex plasmas"  is often used instead, by analogy with "complex fluids", where similar complex phenomena are observed in multicomponent fluids.

Charging and forces
Mixed with a plasma, grains or microparticles collect ions and electrons and typically charge negatively due to higher mobility of electrons (Bronold et al., 2009;Goree, 1994;Melzer et al., 1994). However in the presence of ionising radiation, such as UV light, or thermionic emission, the particle charge may become positive. Typical charges are of the order of ∼ 10 4 electrons for ∼ 10 µm diameter particles. The charge value is proportional to the electron temperature and to the particle radius. The time it takes a particle to reach an equilibrium charge, the charging time, is inversely proportional to the particle size and plasma density (Goree, 1994). Its typical value is ∼ 100 ns for ∼ 10 µm grains in typical laboratory conditions. The particle charge is not constant and it fluctuates around an equilibrium. This may cause instabilities and if the charge value changes its sign, it may even result in particle coagulation.
Highly charged particles are affected by electric fields in the discharge and interact with each other electrostatically. Their interaction potential is usually assumed to be of a Yukawa (Debye-Hückel or screened Coulomb) type, if the background plasma is isotropic (Kennedy & Allen, 2003). This approximation has been also shown to be valid for particles levitating in a plasma sheath at the same height (Konopka et al., 2000) as in monolayer complex plasmas. Flowing plasma makes particle-particle interaction anisotropic with regions of negative and positive potentials (ion wake) (Melandsø& Goree, 1995;Vladimirov et al., 2003) and it also affects the charge of downstream particles.
Apart from the electrostatic force, grains are affected by other forces. Gravitational force becomes dominant for particles with a diameter 1 µm in the bulk of the discharge. Large particles are pushed by the gravity down into the plasma sheath, where the electric field is strong enough to levitate them. This effect makes it necessary to use microgravity conditions in order to produce large three dimensional (3D) structures. Neutral drag force results from collisions with the gas molecules. It is equivalent to friction and damps particle motion. Streaming ions affect grains via an ion drag. This force is responsible for a void formation (Goedheer et al., 2009;. Thermophoretic force arises due to a temperature gradient. It can be used for particle levitation (Rothermel et al., 2002). Intense light sources create a light pressure force, which is utilised for grain manipulation (Liu et al., 2003).
Complex plasmas can be characterised by two parameters. The coupling parameter Γ = U/T is the average ratio of the electrostatic potential energy to the kinetic energy of particles. The screening parameter κ = a/λ D is the ratio of the interparticle distance to the screening (Debye) length. These parameters determine if the complex plasma is in the crystalline or a liquid state (Hamaguchi et al., 1997;Ikezi, 1986;Vaulina et al., 2002). Crystalline plasmas, which have long range correlations, are characterised by large values of the coupling parameter Γ 170 (Ikezi, 1986). Liquid phase state has smaller values 1 Γ 170 and short range correlations. Solids and liquids are strongly coupled states. Gaseous state is weakly coupled (Γ < 1), with uncorrelated particle positions.

Natural occurrence and significance
Dust is abundant in space, where it is found in planetary rings, comet tails, interstellar clouds, and planetary nebulae. The sources of ionisation are also present, such as charged particles from cosmic rays and stellar winds, gas ionised by the stellar radiation, stellar radiation itself, and various radioactive elements in the dust, which emit charged particles and ionising gamma rays. Dusty plasma effects are believed to be involved in formation of dark spokes in Saturn rings (Hartquist et al., 2003). Charge fluctuations are known to enhance coagulation of particles (Konopka et al., 2005). This effect may influence the models of planet formation. Spacecraft and satellites are often charged by the solar wind (Whipple, 1981). This can cause their malfunction due to electrical breakdown. It also increases the drag force due to enhanced collisions with ions. Lunar and Martian dust can be charged by the solar radiation and levitate above the planet surface (Sternovsky et al., 2002). It poses threat to machinery and spacesuits because of its abrasive properties. Charged ionospheric aerosols affect radio wave propagation (Cho et al., 1996) and often disrupt communications. Ultrafine charged particles in the Earth atmosphere influence cloud formation by providing centres of condensation and thus affect the Earth radiative budget (Boulon et al., 2010) with implications for climate models. Aerosol charging modifies the atmospheric chemistry as well as the formation and transport of pollutants (Aikin & Pesnell, 1998).

Industrial applications
Particles with designed properties are grown in a plasma environment for various technological applications (Boufendi et al., 2011).
These include production of fine powders for ceramics and catalysts, phase separated materials, coatings for solar cells, and nanocoatings for optics. Undesirable dust growth has been observed in ultra-clean etching reactors. As the size of the semiconductor device features has approached 22 nanometers, a single dust particle of a similar size can destroy a whole device, significantly reducing the yield of the manufacturing process. This requires strict measures to prevent dust formation. Fusion devices were found to produce metal dust in significant quantities by evaporating or sputtering their walls. This dust is flammable, radioactive and poses safety hazard. It can contaminate and quench the plasma reducing energy yield.

Complex plasmas as model systems
One of the most interesting laboratory uses of complex plasmas is to study properties of solids and liquids at the microscopic or kinetic level. Complex plasmas possess a unique 247 Molecular Dynamics Simulations of Complex (Dusty) Plasmas www.intechopen.com combination of properties and share many of them with other model systems such as colloids and granular media. Grains in complex plasmas can be easily observed. They have sizes of a few microns and separation distances of almost a millimetre. Thus complex plasmas remain optically thin over many interparticle distances. Illuminated with a laser, particles can be imaged with a video camera producing images with high contrast. Typical timescales of particle motion are of the order of 10 ms, meaning that a moderately high speed camera is adequate. The damping rate due to the background gas can be below 1 s −1 .T h i sm a k e si t possible to observe grain-mediated wave motion. The particle-particle interaction potential is continuous and long range, qualitatively similar to that in real solids and liquids. Complex plasmas exhibit a multitude of dynamic phenomena such as waves Zhdanov et al., 2003), solitons (Durniak et al., 2009;Samsonov et al., 2002), shock waves (Luo et al., 1999;Samsonov & Morfill, 2008), melting and crystallisation (Knapek et al., 2007), Mach cones , diffusion (Nunomura et al., 2006), heat transport (Nosenko et al., 2008;Nunomura et al., 2005a), and shear flows (Hartmann et al., 2011). Here we will focus on MD simulations of these phenomena in strongly coupled complex plasmas.

Numerical models for complex plasmas
Computer simulations are used to describe and predict the behaviour of systems whose complexity makes analytical treatments impossible or very difficult. Numerical models advance our understanding of what basic processes are responsible for different observable phenomena. They can be rerun using exactly the same initial conditions  with altered physical processes (e.g. forces) or parameters (e.g. damping rate). Examples of numerical simulations of complex plasmas include dynamics of bilayers , self diffusion in two dimensional (2D) liquids , phase transition between solid and liquid (Farouki & Hamaguchi, 1992), phonons in a linear chain (Liu & Goree, 2005a), diffusion in 2D liquids (Ott & Bonitz, 2009a), defect dynamics , shear flows (Sanbonmatsu & Murillo, 2001), and nonlinear wave propagation (Durniak et al., 2009). There are two basic types of numerical simulations: stochastic and deterministic. The Monte Carlo method (Sheridan, 2009b) belongs to the first type, while the molecular dynamics (Allen & Tildesley, 1987) and the fluid model (Goedheer et al., 2009) to the second. MD simulations often incorporate stochastic elements in order to simulate Brownian motion and the effects of finite temperature. All these techniques are used to simulate complex plasmas, however here we will only consider the MD method.
In order to simulate thermodynamic processes correctly, canonical or microcanonical ensembles are used. This requires conservation of certain thermodynamic quantities (volume, particle number, temperature, energy, etc.) by the algorithm. Numerical errors tend to accumulate and increase the total energy of the system. Energy conserving integrators, such as Verlet and Beeman, are often utilised in order to keep the total energy of the system constant.
Since a relatively small number of particles is involved in the simulation, their kinetic energy and thus temperature tend to have significant fluctuations. Constant temperature is maintained by deterministic techniques such as Nosé-Hoover thermostat (Liu et al., 2006) and velocity rescaling (Ohta & Hamaguchi, 2000;Ott & Bonitz, 2009a), or by stochastic methods such as Andersen thermostat (Nelissen et al., 2007) or Langevin dynamics (Schveȋ gerteta l ., 2000). Langevin dynamics is an extension of the MD method, in which a random (Langevin) force and a damping term are added to the equations in order to simulate the effects of liquid or gaseous background on micron-sized particles. This accounts for random kicks by fast moving molecules as well as for the friction force caused by the liquid or gas drag.
Computing pair interactions for all possible pairs of particles in an ensemble is very costly, when the number of particles is large. Since the interaction force decreases with the distance, it is possible to simplify the calculations. Short range interactions allow introduction of a cut-off distance (Vaulina & Dranzhevski, 2006) beyond which the forces are neglected. This method needs to maintain a list of neighbours to keep track of interacting grains in order to increase efficiency. Long range forces can not be easily truncated, however they can be averaged to reduce the computational cost. Several methods are used such as the Ewald summation (Ott et al., 2011), the fast multipole, the tree code, and particle-mesh-based techniques (Frenkel & Smit, 2002;Rapaport, 1995), and an example of the latter the particle-particle particle-mesh method . The number of particles in the simulation can be reduced using periodic boundary conditions: a particle exiting a simulation box from one side will re-enter on the opposite side. Particles interact not only with other particles in the simulation box but also with particles in image boxes Hamaguchi, 1999). Free boundaries are also often used in complex plasma simulations.

Interaction and confinement
Since plasmas contain a mixture of species moving on very different time scales, it is impossible to simulate a meaningful number of them as individual particles even using a supercomputer. Thus complex plasmas are typically simulated as microparticles interacting via an effective potential. The ion-electron plasma is not explicitly included and enters only as a screening parameter in the interaction potential. The neutral gaseous background is simulated as a friction force. This approach is valid in most cases, as comparisons with experiments show. However there were some notable attempts to include all plasma components (Ikkurthi et al., 2009;Joyce et al., 2001). It is also frequently assumed that the particles have a constant charge, which is often a very good approximation for large and highly charged grains. The most frequently used particle-particle interaction potential in complex plasma simulations is Yukawa, however Coulomb (Lai & I, 1999) as well as a repulsive-attractive 6 Will-be-set-by-IN-TECH potentials have been used (Chen , Yu & Luo, 2005). Other forces can be added e.g. due the influence of magnetic fields , ion drag (Ikkurthi et al., 2009), or external excitations (Durniak et al., 2009). External excitations applied to an equilibrated structure allow to investigate non-equilibrium dynamics , and various dynamic phenomena (Jefferson et al., 2010).
Since Yukawa potential is purely repulsive, an additional confinement is needed, if free boundary conditions are used. The most common confining potential is parabolic (Durniak et al., 2009;Ma & Bhattacharjee, 2002;Nelissen et al., 2007), but other shapes have been used such as 4 th order (Lai & I, 1999), 10 th order , soft-wall (Ott et al., 2008), and hard-wall (Klumov et al., 2009). Special shapes of confinement are used in order to obtain particular particle arrangements, e.g. a 2D annular potential to obtain rings (Sheridan, 2009a) and anisotropic potentials for elliptical clusters (Cândido et al., 1998), linear chains (Liu & Goree, 2005a), or flat disks .
The Langevin dynamics method includes a stochastic force and a dynamic damping force into the equations of motion. The stochastic force does not depend on the particle momentum. It has a zero mean value and a Gaussian probability distribution with a correlation function: where ν is the damping coefficient, k B is the Boltzmann constant, T the temperature of the system and the indices i and j are linked to particles. The damping force is often chosen to be equal to the gas drag force, which depends on the gas pressure, the kind of gas and the momentum exchange between the gas molecules and the grain surface (Epstein, 1924).

Our simulation code
The code used to simulate the complex plasmas presented in this chapter is based on an objects-oriented multi-threaded programming. It is assumed that the microparticles comprising a complex plasma interact with each other via a Yukawa potential, their motion is damped by a neutral drag, and that they are confined in an external potential well with free boundaries. A Langevin force is used to study effects of finite temperature. The ions and electrons are not explicitly included in the model. We take into account the interaction of every microparticle with every other one (particle-particle code), thus there is no cut-off for the potential. The equations of grain motion are solved using a fifth order Runge Kutta method with the Cash Karp adaptive step size control (Press et al., 1992). This makes the code precise, simple and stable at the expense of computational efficiency. It is used to simulate a wide range of dynamic phenomena in a system of several thousand particles.
The model is based on the Newtonian equations of motion written for each microparticle: (2) where m is the particle mass, r = x + y + z is the particle coordinate (z = 0for2D)withthe subscripts s and j denoting different particles, f fr s is the friction force due to collisions with neutrals, f int s is the grain-grain interaction force, f ext s is the external excitation force, L s (t) is the Langevin force (1), ν is the damping rate, Ω h and Ω v are respectively the horizontal and vertical confinement parameters of the well (Ω v = ∞ in 2D case), U 0 is the Yukawa interaction potential, λ D is the Debye screening length, Q is the particle charge and r sj = |r s − r j | is the intergrain distance. The overdots denote time derivatives. Q and λ D are kept constant during the simulation. Dimensionless units are used: the lengths are normalised in terms of the screening length λ D and the time in terms of t = 4πǫ 0 mλ 3 D /Q 2 . The units are converted into the dimensional units after the simulation is completed, to facilitate comparison with the experiments. The code records the position, velocity, and potential energy of each particle at specified time steps.
After seeding the grains randomly, the code is run with the external and Langevin forces switched off until the equilibrium is reached and a monolayer crystal lattice or a solid 3D cluster is formed. The structural properties of the resulting crystals are characterised before they are utilised as inputs for simulations of dynamic phenomena. These simulations are performed by applying various excitation forces. A random (Langevin) force is used to simulate the thermal Brownian motions of particles and to obtain phonon spectra. Pulsed excitations are applied to investigate nonlinear waves and structural properties of complex plasmas.
We use the following parameters m = 5 × 10 −13 kg, λ D = 1 mm, Q = −16000e (where e is the electron charge) in all our simulation runs. Other parameters are listed in Table 1.

Data and structural analysis
The results of complex plasma simulations are analysed in order to determine their basic microscopic and macroscopic parameters, structural and dynamical properties using standard methods, which are also used for analysing the experimentally obtained particle tracks.
The local orientation of 2D crystalline cells is characterised using the local bond orientational order parameter for each lattice cell: ,whereN is the number of nearest neighbours, θ j is the angle between the x-axis and the bond connecting the central particle with its neighbour j. The average bond orientation angle θ 6 is used to highlight crystal grains separated by strings of defects as well as lattice deformations. The value of |ψ 6 | gives the local order parameter, which is equal to one for an ideal crystal.
Delauney triangulation of the lattice is used in order to find the nearest neighbours of each particle and determine their numbers. An ideal hexagonal lattice would have particles with 6 nearest neighbours. Lattice defects are defined as lattice cells that have other numbers of neighbours, such as 5 or 7. Dislocations are pairs of 5-and 7-fold cells, they are also called penta-hepta defects and are characterised by their Burgers vectors, perpendicular to the axis formed by the two defective cells. Defects and dislocations play an important role in melting and plastic deformations (Durniak & Samsonov, 2011;Knapek et al., 2007).
phonon 2D  , ν is the damping rate, κ is the screening parameter, is the initial crystal diameter, F ex0 is the amplitude of the excitation force F ex , x 0 is the force offset with respect to the centre of the lattice, w is the width and τ is the duration of the force. Continuous randomly changing excitation force is used to simulate phonon spectra (case 1). Pulsed excitation force is applied on one (cases 3-5,10) or both (cases 2,6) sides of the lattice. Its spatial profile is profile (cases 2-5,10) is a parabola (inverted and truncated at negative values) It is modified in case 6 to keep a constant value after the maximum is reached (F ex = const for t ≥ τ). For cases 1-8 and 10, the excitation force is independent of y. The Mach cones are simulated using an anisotropic Gaussian excitation force The number density of the lattice is computed using the Voronoi analysis (Voronoi, 1908). A Voronoi cell is defined as a set of points for which a given particle is the nearest. The local number density is proportional to the inverse area in 2D (volume in 3D) of a Voronoi cell. The compression factor of a lattice is calculated as the ratio of the number density n of a strained lattice to its unperturbed number density n 0 .
The kinetic temperature of a complex plasma is determined from the velocities of individual particles. The lattice is split into bins and the average bin velocity v is calculated. It is then subtracted from the speeds of all particles. The average kinetic energy E in the bin is determined using the mean square random velocity wherek B is the Boltzmann constant and d is the number of the degrees of freedom.
A correlation analysis is performed in order to assess the structure of complex plasmas. The pair correlation function g(r) gives the probability of finding a specific distance between two particles in the system relative to the probability of finding that distance in a completely random particle distribution of the same density (Crocker & Grier, 1996;Quinn et al., 1996). It measures the translational order of the lattice. For a perfect crystal at zero temperature g(r) is a series of δ-functions. At non-zero temperature peaks have finite widths and decaying amplitude. The position of the first peak of the pair correlation function is determined by the interparticle spacing and is often used to measure the average interparticle distance.

Experimental verification
Like any other simulations, numerical models of complex plasmas have to be verified by comparison with the experiments. The code should reproduce the same phenomena as observed in the experiment with similar quantitative characteristics.

Laboratory complex plasmas
The most frequently used method to obtain high quality complex plasmas in the laboratory is to add premanufactured monodisperse plastic microspheres to the plasma (Thomas et al., 1994). Since the crystalline state is most often the subject of interest, relatively large (∼ 10 µm in diameter) particles are used. The reason for that is the lower boundary of the coupling parameter (Γ∝Q 2 /T 170) required for crystallisation (Sect. 2.1). The value of the kinetic temperature can not be lower than the temperature of the neutral gas, which is close to the room temperature. Thus the coupling parameter can be only increased by increasing the potential energy of the grains or their charge Q, which is a function of the grain size. The dependence of the charge on the particle size makes it also necessary to use monodisperse particles in order to make their charges, potential energies, and levitation heights identical.
Most ground-based experiments utilise a capacitively coupled radio-frequency (rf) gas discharge Durniak & Samsonov, 2011;Käding et al., 2008). The particles are suspended in the plasma sheath, where a strong electric field counteracts the gravity. They typically form between one (Nosenko & Zhdanov, 2009) and a few layers (Quinn & Goree, 2001) creating 2D and quasi-2D structures. In order to obtain 3D structures, gravity has to be compensated (Sect. 2.1). This can be done in a strong electric field such as in striations of a 253 Molecular Dynamics Simulations of Complex (Dusty) Plasmas www.intechopen.com direct current glow discharge (Fortov et al., 2005b), using a thermophoretic force (Arp et al., 2004), or under microgravity conditions on parabolic flights or on board International Space Station (Konopka et al., 2005;Seurig et al., 2007). A Q-machine (Luo et al., 1999) has also been used to create weakly coupled complex plasmas using polydisperse particles.

The experimental setup
Our experiments are performed in a capacitively coupled rf discharge vacuum chamber as shown in Fig. 1. An argon flow (a few sccm) maintains a constant working gas pressure (1-2 Pa) in the chamber. An rf power is applied to the lower disc electrode, which is 20 cm in diameter. The chamber itself is the other grounded electrode. Due to different area of the electrodes and different mobility of ions and electrons, the powered electrode has a negative self-bias voltage, which helps to suspend the particles in the plasma sheath against the gravity. The particles used are monodisperse plastic microspheres 8.9 or 9.19 µmi nd i a m e t e r . T h e y are injected into the plasma through a particle dispenser, levitated in the plasma sheath, and confined radially by a rim on the outer edge of the electrode, forming a monolayer hexagonal lattice of approximately 6 cm in diameter. The particles are illuminated by a horizontal thin (0.2-0.3 mm) sheet of laser light and imaged by a top-view digital camera. Two parallel horizontal tungsten wires, both 0.1 mm in diameter are placed below the particles. Negative pulses applied to one or both wires excite compressional disturbances and deformations.
The experimental results are analysed by identifying the particle positions in all video frames using the intensity weighted moment method (Feng et al., 2007;Ivanov & Melzer, 2007). The particle velocities are calculated by tracking their positions from one frame to the next.

Structure of complex plasmas
Complex plasmas are formed into different shapes by the confinement potential (see Sect. 3.2). Parameters of the confinement should be selected to ensure that stable structures are formed. The structures of interest include linear 1D chains, monolayer 2D lattices, and 3D balls. Different parameter values result in formation of zigzag lines instead of chains (Sheridan, 2009a), or multiple layers instead of monolayers. The only structure stable against shear perturbations in crystalline monolayers is hexagonal . Several different structures exist in bilayer lattices (Donkó & Kalman, 2001) and transition between them is controlled by the layer separation. As the separation increases, the structure changes from hexagonal to staggered square and then to staggered rhombic. More crystal structures are observed in multilayer and 3D strongly coupled systems. Face centred cubic (fcc), body centred cubic (bcc), and hexagonal close packed (hcp), illustrated in Fig. 2, have been identified in (Hamaguchi et al., 1997;Klumov et al., 2009). Phase diagrams of the phase transitions between liquid, bcc and fcc phases have been simulated by (Hamaguchi, 1999).

Linear and nonlinear waves
Complex plasmas sustain waves, which are analogous to those in ordinary solids and liquids. Small amplitude linear waves as well as nonlinear waves, solitons and shocks are observed.
Most wave experiments in crystalline complex plasmas are performed in monolayers, since they are easy to obtain at low gas pressure. This results in very low damping and therefore underdamped wave motion can be studied in great details.

Small amplitude waves
Complex plasmas can be in the solid or liquid state (Sect. 2.1). Crystalline monolayer complex plasmas sustain acoustic compressional and shear wave modes, as shown experimentally and numerically in Zhdanov et al., 2003). The phonon spectra are isotropic for long wavelength phonons, but strongly anisotropic in the short wavelength case. The

255
Will-be-set-by-IN-TECH compressional mode has a higher frequency than the shear mode. As the complex plasma crystal melts, the shear mode disappears, while a dust thermal mode appears (Nunomura et al., 2005). Brownian simulation has been used to calculate the cut-off wavenumber for the shear mode in a liquid state and to compare the wave spectra in 2D complex plasma solids and liquids with those from different analytical theories (Hou, Misković, Piel & Murillo, 2009). It confirmed the existence of the dust thermal mode first observed by (Nunomura et al., 2005). Lattices formed of particles with two different sizes exhibit optical compressional and shear modes. Simulation of bilayers revealed the existence of the optical modes and verified the existence of the k = 0 energy (frequency) gap ).
The method to compute phonon spectra from either experimental or simulated data relies on the Fourier transform of particle velocities v(x, t) both in time t and in the x-direction Zhdanov et al., 2003): where k and ω are the wave number and the frequency, N and M are the numbers of data points in space and time respectively, L o is the length of the field of view and T o is the recording period. The longitudinal and transverse modes are resolved by using v x and v y components of the velocity respectively. The resulting phonon spectra are shown in Fig. 3. They are obtained using a MD simulation (Table 1, case 1) with a stochastic force. The theoretical dispersion relations Hou, Misković, Piel & Murillo, 2009) are calculated for a perfect hexagonal lattice by solving the eigenvalue problem for the dynamical matrix D µν : where the summation is performed over all particles j with coordinates r j and mass m; U 0 (r j ) is the Yukawa potential (Eq. 2), and coordinates µ, ν take the values {x, y}. The spectra depend on the propagation angle or the direction of k. The longitudinal L and transverse T modes are given by: In the limit of long wavelengths the dispersion relations (Eq. 3) become ω L,T = c L,T k,w h er e c L,T are the longitudinal and transverse wave speeds, which depend on the particle charge and the screening parameter. This can be used to determine the values of Q and κ in experiments.

Nonlinear pulses and solitons
Pulsed excitation applied to a lattice with a laser or a biased wire results in a localised propagating disturbance. This disturbance can be compressional (Nosenko et al., 2002;Samsonov et al., 2002) or shear . is its width) and a relation between the propagation speed and the amplitude (faster solitons have larger amplitudes). This has been observed in a 2D experiment and a linear chain simulation (Samsonov et al., 2002). Interactions between nonlinear waves is another subject of interest. Experimental and numerical investigation of collisions of counter-propagating solitons in complex plasma monolayers find that solitons with larger amplitude experience larger delays and that the amplitude at the collision point is different from the sum of the initial soliton amplitudes (Harvey et al., 2010). Figure 4 shows a head-on collision of two solitons simulated with the parameters of Table 1, case 2. The amplitudes of the pulses slightly decrease due to the neutral damping as they propagate. The amplitude of the overlapping solitons is lower than the sum of the initial amplitudes.
Waves can gain amplitude due to nonlinear effects even in the presence of damping. A soliton propagating in a lattice with decreasing number density gains amplitude as observed in simulation (Table 1, case 3) and experiment (Fig. 5). It is found that the measured amplitude gain is higher than that predicted by the KdV equation with damping included (Durniak et al., 2009).

Shock waves
Shock waves are propagating discontinuities arising from large amplitude perturbations, therefore they cannot be treated as small amplitude waves. Shocks can cause phase transitions and present a challenge for simulations since they are likely to cause numerical instabilities. Experimental observations of shock waves in complex plasmas were reported in (Fortov et al., 2005b;Samsonov et al., 2003;. The structure of a simulated shock (Table 1, case 4) is shown in Figure 6. The shock front has a thickness of a few interparticle distances and an oscillatory structure. It propagates from left to right at a velocity decreasing from 57 to 46 mm/s or a Mach number varying from 1.9 to 1.6. The lattice, initially in the solid state, is melted after the shock. There is a discontinuous jump at the shock front in compression factor, number density, kinetic temperature, and defect fraction Samsonov & Morfill, 2008).

Defect dynamics and plastic deformation
Lattice defects and dislocations determine mechanical properties of crystals and may be responsible for material fatigue and catastrophic failure. The effect of temperature on defects in 2D Coulomb clusters has been studied numerically by (Nelissen et al., 2007). It was found that the defect mobility strongly depends on the neighbouring defects, that the geometrical defects have different dynamics than the topological defects, and that a fast cooling rate favours formation of a non-equilibrium glass-like state with many defects. Dislocations have been observed to propagate in crystalline complex plasmas supersonically (Nosenko et al., 2007). They interact with nonlinear waves  and they are generated during shear slips in plastic deformations (Durniak & Samsonov, 2011). The direction of the wave excitation force F ex (and thus the wave propagation direction) is shown with the same colour as the corresponding wave trajectories. The dislocations move either in small increments due to elastic deformations of the lattice, or in large jumps due to lattice structural changes. One can see that the defect jumps occur in directions almost parallel to the Burgers vector regardless of the wave propagation direction. The crystal lattice (b) is initially characterised by a screening parameter κ = 0.725. The defects are marked by (7-fold defect) and (5-fold defect) and the Burgers vector by the pink arrow.
The interaction of a dislocation with a wave is simulated using parameters listed in Table 1, case 5. Figure 7a shows the trajectories of an isolated dislocation as the wave passes by from different directions (runs 1-4). The excitation wave has an amplitude between 4.8 and 6.1 mm/s in all runs and propagates at about 38 mm/s. The dislocation either stays at the same lattice site, or jumps to a neighbouring pair of particles. In the former case it displaces roughly in the direction of the wave, while in the latter case it moves almost parallel to its Burgers vector. This result agrees with the experiment .
Plastic deformation under a uniaxial compression is numerically modelled using the conditions of Table 1, case 6. As the strain increases and exceeds the elastic limit, shear slips occur causing stress relaxation. This happens because a uniaxial compression is a superposition of a uniform compression and shear. Complex plasmas are very compressible and do not change their structure under a uniform compression. However their shear strength is not very high, thus their structural failure results in a shear slip. Shear slips are initialised by generation of a pair of dislocations, which move in opposite directions (Durniak & Samsonov, 2011). (b) Time evolution of the particle trajectories shows that the local stress is relaxed by a shear slip (trajectories twist at time 2 s). The crystal defects are marked in blue (7-fold defect) and green (5-fold defect). Only the particles along the slip line are shown. Voronoi maps (c-e) with defects marked by △ (5-fold) and (7-fold) visualise the lattice structure at (c) t=1.5 s, (d) t=2.04 s, and (e) t=2.3 s. The colour scale shows the bond orientation angle θ 6 . A pair of dislocations is generated at t=2.04 s, they separate and move in opposite directions as the slip progresses. Velocity vector maps show the particle velocities at (f) t=1.5 s, (g) t=2.04 s, and (h) t=2.3 s. The colour scale indicates the angle between the velocity vector and the horizontal axis. The slip is shown as particle rows moving in opposite directions.

Coulomb (Yukawa) clusters
Coulomb (or Yukawa) clusters are systems of up to a few hundred charged particles confined in a 2D or 3D well and interacting via a Coulomb (or Yukawa) potential. They occur in systems of trapped electrons, ions, colloids, and complex plasmas. Since they comprise a small number of particles, they are easy to model numerically (Nelissen et al., 2007). Properties of clusters are different from those of the bulk material and they depend on the cluster size (Yurtsever et al., 2005), the interaction and confinement potentials. The emergence of bulk behaviour in a strongly coupled Yukawa cluster has been studied in (Sheridan, 2007). Clusters are often simulated using the Monte Carlo method, however we will focus here on the MD technique. Figure 9 shows the structure of 2D and 3D clusters simulated using the parameters of Table 1, cases 7 (2D) and 8 (3D). The particles in these clusters interact via a Yukawa potential and are confined in an isotropic parabolic potential. The structure of the clusters results from the interplay between the particle-particle interaction and the global confinement. The interaction potential favours hexagonal order whereas the confinement induces a circular symmetry. Thus large clusters have a hexagonal inner core and circular outer shells, while small ones tend to contain only shells (Lai & I, 1999). The same effect is observed in 3D clusters (Arp et al., 2004;Totsuji et al., 2002). Clusters assume particularly stable configurations at certain "magic" numbers of particles (Tsuruta & Ichimaru, 1993). The effect of the screening length on the structure of Coulomb balls has been studied in (Bonitz et al., 2006;Käding et al., 2008): particles moved from the outer shells to the inner ones as the screening length increased.
Metastable states of 3D Yukawa clusters can occur with a significantly higher probability than the ground state. The results strongly depend on the screening parameter and the damping coefficient. Slow cooling favours the ground state over the metastable ones . The effects of anisotropic confinement and interaction potentials have been studied by (Killer et al., 2011). The structure of spherical clusters was found to be unaffected by a weak ion focus unlike the structure of elongated clusters.
Since clusters have a finite size, they have a finite number of normal modes of oscillations.
Knowing the normal modes of a system allows to determine its response to external excitations, i.e. its dynamics. Clusters of N particles in a harmonic potential well have 2N normal modes. Oscillations of clusters are not purely compressional or shear, however they can be described as compression-like or shear-like (Melzer, 2003). It was shown that in asymmetric potentials both the rotational and breathing modes of elliptical clusters were robust . Melting and defect excitation in Coulomb clusters have been also investigated (Lai & I, 2001;Nelissen et al., 2007).

Mach cones
Mach cones are propagating V-shaped disturbances. Their existence in a dusty plasma of Saturn rings was first predicted by (Havnes et al., 1995). They were then observed in a monolayer complex plasma , where they were generated by fast particles moving parallel to the main layer. They had a multiple V-shaped structure formed by compressional waves. The properties of Mach cones were confirmed using laser excitation . The opening angle of Mach cone obeys the relation µ = sin −1 (1/M), M = v/c being the Mach number of an object moving at speed v through a medium with an acoustic speed c. Measurements of this angle can therefore be used as a diagnostic tool for complex plasmas to detect inhomogeneity . Since there are two acoustic wave modes in plasma crystals, compressional and shear (see Section 5.2.1), two types of Mach cones exist. They are distinguished by the particle motion in the cone front. The particles move  The red lines correspond to the particle bonds calculated by a Delauney triangulation, which reveals the static structure of the clusters. The numbers of particles N in the 2D clusters are shown in (a). The 3D cluster (b) has been been generated using 150 particles. Its shell structure is visualised in (c) by plotting the particle positions in cylindrical coordinates (r = x 2 + y 2 , z). In (b) the colour scale corresponds to the particle's vertical position z.
perpendicular to the cone front in compressional Mach cones and parallel to the front in shear Mach cones (Ma & Bhattacharjee, 2002;Nosenko et al., 2003).  Figure 10 shows a numerically generated Mach cone with parameters listed in Table 1 case 9. The excitation force has been chosen similar to that reported by (Ma & Bhattacharjee, 2002; 262 Molecular Dynamics -Studies of Synthetic and Biological Macromolecules www.intechopen.com Nosenko et al., 2003). We visualise the cones by plotting the velocity vector map (Fig. 10a) and also the vorticity and divergence maps of simulated particle positions, which highlight the shear (Fig. 10b) and compressional (Fig. 10c) cones respectively. The compressional cone has a multiple structure while the shear cone is single. This is due to the fact that the shear wave is slower and thus it can propagate a shorter distance than the compressional before being damped. The angles of the inner and outer compressional cones might be slightly different due to nonlinear effects .

Phase transitions
MD simulations of charged grains in plasmas have been used to study phase transitions before this became possible experimentally (Farouki & Hamaguchi, 1992). Solid and liquid phases have been predicted as well as a hysteresis at the transition between them. This hysteresis corresponds to a superheated solid and supercooled liquid. Solid superheating was later observed experimentally (Feng et al., 2008). Phase diagrams of Yukawa systems have been computed in (Robbins et al., 1988) predicting liquid as well as solid fcc and bcc structures.
Simulations show that strongly screened Yukawa systems have a triple point or a point on the phase diagram where liquid, bcc, and fcc phases coexist, whereas weakly screened Yukawa systems do not have a triple point (Hamaguchi et al., 1997). In experiments, phase transitions can be induced by stochastic laser heating, shear flows (Nosenko & Zhdanov, 2009), shock waves (Knapek et al., 2007) or by changing the discharge power (Rubin-Zuzic et al., 2006). In the latter case a propagating crystallisation front has been observed. It was shown that this process is fundamentally 3D (Klumov et al., 2006). Complex plasma recrystallisation has been simulated by . It was found that the sizes of crystal domains have a power-law time dependence. Figure 11 shows melting and recrystallisation of a complex plasma lattice excited by a shock wave. The simulation parameters are given in Table 1 case 10. As the shock propagates, the kinetic temperature increases and defects are generated ( fig. 11a)). The pair correlation function calculated at different times indicates that the order in the system decreases during melting and then increases reaching almost the initial level at 6 s as the lattice recrystallises.

Transport phenomena
Transport phenomena are irreversible statistical processes which are responsible for transfer of mass, momentum, or energy in matter. These processes use similar mathematical formalisms and are described by similar equations. The three most commonly considered transport phenomena are diffusion (mass transfer), heat conduction (energy transfer), and viscosity (momentum transfer). Their fundamental nature makes them very important for understanding basic properties of matter. Complex plasmas offer a possibility to study these processes at the level of individual particles and compare experiments directly to MD simulations.
Two simulation techniques are used to study transport phenomena: the equilibrium and nonequilibrium methods . The first method calculates the particle trajectories of a system in a state of statistical equilibrium (Vaulina et al., 2008). The second method applies a perturbation to an equilibrated system and measures the changes it causes Voronoi maps visualise the lattice structure at (c) t=0 s, (d) t=0.5 s, (e) t=2 s, and (f) t=6 s. The lattice defects are marked by (7-fold), (5-fold), and * (other). (Sanbonmatsu & Murillo, 2001). If the investigated system exhibits small deviations from the statistical equilibrium, which is the case for equilibrium simulations, the transport coefficients are found from the Green-Kubo relations . Nonequilibrium simulations calculate the transport coefficients directly, by computing the diffusion coefficient from the mean square displacement (Ohta & Hamaguchi, 2000), viscosity coefficient from the velocity profile (Sanbonmatsu & Murillo, 2001), or the thermal conductivity from the temperature gradient ). Both equilibrium and nonequilibrium methods might produce artifacts due to the finite size of the system or insufficient recording time, thus care should be taken. Another possible problem is the existence of the transport coefficients in particular systems. It was found that in a 2D Yukawa liquids the diffusion coefficient exists at high temperature and the viscosity coefficient at low temperature but not in the opposite limits ). The thermal conductivity did not appear to exist at high temperature and it could not be evaluated at low temperature due to computational limitations.
The diffusion is governed by the Fick's law: J = D∇C,whereJ is the diffusion flux, D is the diffusion coefficient, C is the molecular concentration. The diffusion coefficient determines the time dependence of the mean square displacement of all particles |r(t) − r(0)| 2 = 4Dt α . If α = 1 the diffusion is normal, otherwise (α = 1) it is anomalous: α > 1 corresponds to superdiffusion and α < 1 to subdiffusion. The temperature dependence of the diffusion coefficient for 3D Yukawa systems has been studied by (Ohta & Hamaguchi, 2000). It was found to be independent of the screening parameter. Different methods of computing the diffusion coefficient at short observation times have been compared by (Vaulina et al., 2008) and the required minimum observation time has been estimated. Anomalous diffusion has been reported in 2D Yukawa systems in some simulations  as well as experiments (Juan & I, 1998), which contradicts some other numerical (Vaulina & Dranzhevski, 2006) and experimental results (Nunomura et al., 2006). As shown by (Ott & Bonitz, 2009b) the diffusion exponent depends critically on the neutral gas friction, which makes anomalous diffusion a transient effect in simulations. This certainly does not rule out superdiffusion in experiments, however it is possible that limited fields of view and insufficiently long observation times might obscure anomalous effects and the role of nonequilibrium states in the experimental dissipative-driven systems. The effect of coherent transport on the diffusion exponent, e.g. waves and flows is also difficult to rule out.
Shear viscosity in strongly coupled 3D Yukawa system has been studied by (Sanbonmatsu & Murillo, 2001). It was found that the viscosity coefficient has a minimum and that it is nonlocal, with the scale length consistent with the correlation length. The decay of the velocity profile deviates from the one predicted by the Navier-Stokes equation. The minimum of the viscosity coefficient is also observed in 2D Yukawa liquids (Liu & Goree, 2005b). The wavenumber-dependent viscosity, which characterises the viscous effects at different length scales has been computed by (Feng et al., 2011). They have also verified the accuracy of the Green-Kubo relation for static viscosity in the presence of damping. Dynamic shear viscosity has been studied experimentally and numerically by (Hartmann et al., 2011) and shown to exhibit strong frequency dependence. A shear-thinning effect has been demonstrated under static shear.
Heat conductivity in a 2D strongly coupled system has been simulated using nonequilibrium methods by ). The results show that the heat conductivity coefficient depends on the damping rate and indicate that it might be temperature-dependent. The experimental results, however, show no temperature dependence within the experimental uncertainty (Nosenko et al., 2008;Nunomura et al., 2005a).
Since the Newton's law of viscosity, the Fourier's law of heat transfer, and the Fick's law of molecular diffusion are very similar, relations between the transport coefficients can be found. The Stokes-Einstein formula relates the diffusion and viscosity coefficients. It has been thoroughly tested in 3D liquids even at the molecular level. In 3D strongly coupled complex plasmas the Stokes-Einstein relation has been verified for a wide range of temperatures down to the solidification point . However this relation is violated in 2D complex plasmas near the disordering transition, remaining valid at higher temperatures (Liu et al., 2006).

Conclusion
Methods of molecular dynamics simulations of complex plasmas and their results have been reviewed in this chapter and illustrated by examples of such simulations. Complex plasmas belong to the soft matter class of materials; they are formed by mesoscopic and highly charged particles immersed in a plasma. They are underdamped due to their gaseous background and thus ideally suitable for studying dynamic effects at the particle level. Complex plasmas 265 Molecular Dynamics Simulations of Complex (Dusty) Plasmas www.intechopen.com 22 Will-be-set-by-IN-TECH resemble colloids and some effects such as phase transitions, diffusion and shear flows can be observed in both. Some phenomena such as waves, shocks, kinetic temperature, and phonon heat transfer exist only in complex plasmas. The numerical results have been compared with experiments in order to make sure that simulated effects are observable in physical systems. Unfortunately due to space constraints, some effects not exclusive to complex plasmas were left out (such as lane formation) as well as some of those not extensively simulated (such as coupling of vertical and horizontal modes in monolayers). The authors hope that their choice of dynamic effects well illustrates the beauty and complexity of complex plasmas. Molecular Dynamics is a two-volume compendium of the ever-growing applications of molecular dynamics simulations to solve a wider range of scientific and engineering challenges. The contents illustrate the rapid progress on molecular dynamics simulations in many fields of science and technology, such as nanotechnology, energy research, and biology, due to the advances of new dynamics theories and the extraordinary power of today's computers. This second book begins with an introduction of molecular dynamics simulations to macromolecules and then illustrates the computer experiments using molecular dynamics simulations in the studies of synthetic and biological macromolecules, plasmas, and nanomachines. Coverage of this book includes: Complex formation and dynamics of polymers Dynamics of lipid bilayers, peptides, DNA, RNA, and proteins Complex liquids and plasmas Dynamics of molecules on surfaces Nanofluidics and nanomachines