Dynamical correlations and collective excitations of Yukawa liquids
Abstract
In dusty (complex) plasmas, containing mesoscopic charged grains, the graingrain interaction in many cases can be well described through a Yukawa potential. In this Review we summarize the basics of the computational and theoretical approaches capable of describing manyparticle Yukawa systems in the liquid and solid phases and discuss the properties of the dynamical density and current correlation spectra of three and twodimensional strongly coupled Yukawa systems, generated by molecular dynamics simulations. We show details of the dispersion relations for the collective excitations in these systems, as obtained theoretically following the quasilocalized charge approximation, as well as from the fluctuation spectra created by simulations. The theoretical and simulation results are also compared with those obtained in complex plasma experiments.
1 Introduction
Strongly coupled plasmas – in which the average potential energy per particle dominates over the average kinetic energy – appear in a wide variety of physical systems: dusty plasmas, charged particles in cryogenic traps, condensed matter systems such as molten salts and liquid metals, electrons trapped on the free surface of liquid helium, astrophysical systems, such as the ion liquids in white dwarf interiors, neutron star crusts, supernova cores, and giant planetary interiors, as well as in degenerate electron or hole liquids in twodimensional or layered semiconductor nanostructures [1]. Many of these systems share some properties, which allows modeling them by considering explicitly only a single type of charged species and using a potential that accounts for the presence and effects of other types of species. This latter may be thought of as forming a chargeneutralizing background, which is either nonpolarizable or polarizable. In the first case the interaction of the main plasma constituents can be expressed by the Coulomb potential, while in the case of polarizable background the use of the Yukawa potential is appropriate to account for screening effects ( is the Debye length). Perhaps the most important realizations of systems lending themselves to the approximation of the interaction by the Yukawa potential are charged colloids [2, 3, 4, 5] and dusty (complex) plasmas [6] (for comprehensive review on dusty plasmas see e.g. [7, 8]).
In the case of 2D colloidal systems the microscopic particles move in thin liquid films or between two closely separated glass plates.
In dusty plasmas both threedimensional (3D) and twodimensional (2D) settings appear in nature and in laboratory environments. In laboratory experiments 2D systems appear as a layer of dust particles levitated in gaseous discharges. While most of the studies on this latter system have been carried out in the crystalline state (for early references to “plasma crystals” see [9, 10, 11, 12]), the liquid state is receiving more current attention [13, 14, 15, 16, 17, 18]. The important difference between colloidal and dusty plasma systems is in the damping rate in particle dynamics and concomitantly in the wave dynamics. In colloidal suspensions the background liquid exerts a large friction on the moving charged particles, while in dusty plasmas the background is gaseous and therefore the friction is lower and the damping of the waves is weak. For this reason our focus in this Review will be directed at dusty plasmas. The Review will cover studies of strongly coupled plasmas mostly in the liquid state, where both free motion and localization intervene. The principal observation is that from the point of view of collective behavior it is the localization – even though imperfect localization, or quasilocalization – of particles that plays the principal role. In contrast to the Vlasov plasma where the collective modes arise from the fluidlike continuum behavior, in the strongly coupled liquid they are more related to the normal modes of the interacting quasilocalized particles. This, of course, suggests a link with the harmonic phonon theory of crystal lattices. At the same time, one has to allow for the randomness of the distribution of the particles and for the finite lifetime of the localization in the constantly changing potential landscape. This latter process is expected to be primarily responsible for the damping of the collective modes, in contrast both to Vlasov plasmas, where Landau damping dominates and to weakly correlated plasmas where collisional damping is the principal damping mechanism. This physical picture suggests a microscopic equationofmotion model where the particles are trapped in local potential fluctuations. The particles occupy randomly located sites and undergo oscillations around them. At the same time, however, the site positions also change and a continuous rearrangement of the underlying quasiequilibrium configuration takes place. Inherent in this description is the assumption that the two time scales are well separated and that for the description of the fast oscillating motion, the time average (converted into ensemble average) of the drifting quasiequilibrium configuration is sufficient. Here the distinction between the ”direct” and ”indirect” thermal effects should be emphasized: the former are responsible for the actual motion and migration of the particles, the latter refer to the accessibility of the possible configurations of the random sites and to the temperature dependence of the probability of a particular configuration.
The interaction potential energy of particles in Yukawa liquids is given by (e.g. [19]):
(1) 
where is the charge of the particles, is the permittivity of free space, and the Debye length accounts for the screening of the interaction by other plasma species. The main (dimensionless) parameters, which fully characterize the systems are: (i) the coupling parameter (defined in the same way as for Coulomb systems):
(2) 
where is the WignerSeitz (WS) radius and is the temperature, and (ii) the screening parameter
(3) 
is the customary measure of the ratio of the average potential energy to the average kinetic energy per particle; the strong coupling regime, relevant here, corresponds to 1. In the limit the interaction reduces to the Coulomb type, while at it approximates the properties of a hard sphere interaction.
At = 0, in 3D, the liquid domain is limited to coupling parameter values [20, 21, 22], where the plasma is known to crystallize into a bcc lattice [23, 24]. In 2D, crystallization into hexagonal lattice occurs at a lower value of coupling, at , as found by computer simulations [25, 26] and proven by experiments [27] as well. At 0, 3D systems may crystallize either in a bcc or in a fcc lattice, depending on , as found by Hamaguchi et al. [28] in their calculation of the phase diagram of Yukawa systems. In 2D the crystallized form of the systems is always hexagonal.
The possibility of characterizing a Yukawa liquid with an effective coupling coefficient instead of a () parameter pair has recently been addressed. In 3D was derived by Vaulina, Fortov and coworkers [29, 30, 31] on the basis of the frequency of dust lattice waves. Subsequently, for 2D Yukawa liquids the relationship was established by prescribing a constant amplitude for the first peak of the pair correlation function for fixed values of [32]. The solidliquid melting line ( vs. ) in both 3D and 2D system was found in these studies to follow closely constant values.
Additional characteristic parameters of the systems investigated here are the WS radius and the plasma frequency , which are given for 3D and 2D systems by:
(4)  
(5) 
where and are the 3D (number) density and the 2D areal (number) density of particles, and
(6) 
and
(7) 
Note that in 2D the nominal plasma frequency may also have different definitions, and some of the authors use the lattice constant instead of the WS radius as a length scale.
Other important frequencies characterizing strongly coupled plasmas are the Einstein frequencies, which are the normal modes of oscillation of a test charge in the presence of a given (static) distribution of charges. Einstein frequencies are well known for lattice structures, however, there has been relatively little work done on disordered and liquid phase systems [33]. Such systems are being studied through the combination of theoretical and simulation approaches [34, 35].
As to the possibilities of theoretical description, manybody systems can be treated theoretically in a straightforward way in the extreme limits of both weak interaction and very strong interaction. In the first case, one is faced with a gaseous system, or a Vlasov plasma, where correlation effects can be treated perturbatively (). In the case of very strong interaction, the system crystallizes, the particles are completely localized and phonons are the principal excitations. In the intermediate regime – in the strongly coupled liquid phase – the localization of the particles in the local minima of the potential surface still prevails, however, due to the diffusion of the particles the time of localization is finite [34]. The localization of the particles (which may typically cover a period of several plasma oscillation cycles) serves as the basis of the QuasiLocalized Charge Approximation (QLCA) method [36, 37]. Besides the theoretical approaches computer simulations have proven to be invaluable tools for investigations of strongly coupled liquids of charged particles. Monte Carlo (MC) and molecular dynamics (MD) methods have widely been applied in studies of the equilibrium and transport properties, as well as of dynamical effects and collective excitations. The main difference between the two techniques is that in an MC simulation the particle configuration with the lowest energy is searched for, whereas MD simulations provide information about the timedependent phase space coordinates of the particles, this way allowing studies of dynamical properties.
This paper intends to review the dynamical properties and collective behavior of strongly coupled Yukawa systems in the liquid and solid phases, in two and three dimensions. First we describe the numerical as well as theoretical methods used, in Sections 2 and 3, respectively. The analysis of the collective mode behavior of 3D liquids is presented in Section 4.1. In the 2D case we investigate both an ideally narrow particle layer, and a layer having a finite width, where particles are confined by an external parabolic potential. The analysis of these systems is described in Sections 4.2 and 4.3, respectively. Section 5 gives a brief summary of the experimental studies relevant to the theoretical and simulation results reviewed. Section 6 gives a summary of the paper as well as a short outlook on the topics, which may have been additional subjects of this Review.
2 Molecular dynamics simulations
Molecular dynamics simulations follow the motion of particles by integrating their equations of motion while accounting for the pairwise interaction of the particles, as well as for the forces originating from any external field(s), see e.g. [38]. In the plasma / gas background environment friction forces and randomly fluctuating forces also act on the particles in addition to the forces arising from the interaction of the electrically charged particles. The general form of the equation of motion (of a “testÃÃ particle ) is (see e.g. [39]):
(8) 
where is the force originating from the interaction with particle , is the force originating from any external field, is the friction coefficient, and represents a Brownian randomly fluctuating force (Langevin force). The results presented here correspond to “idealized” Yukawa liquids for which and are assumed. Also, in most of our studies we investigate infinite (unconfined) systems, for which , as well, although as an example a quasitwodimensional liquid – confined by an external parabolic field – is also studied. In the case of unconfined systems periodic boundary conditions (PBC) are imposed in the simulations. In the case of confinement along one of the coordinates, PBCs are used in the unconfined directions. It is noted that for charged colloids Brownian molecular dynamics simulation [3, 40] is widely used, which represents the extreme limit of large friction and large . In this case the inertial term () in eq. (8) is neglected.
The calculation of the force acting on a particle of the system, , is relatively simple in the case of shortrange potentials (like the LenardJones potential or the potential). In this case MD methods make use of the truncation of the interaction potential thereby limiting the need for the summation of pairwise interactions around a test particle to a region of finite size. In the case of longrange interactions (e.g. Coulomb or low Yukawa potentials), which are also of interest here, however, such truncation of the potential is not allowed, and thus special techniques, like Ewald summation [41], have to be used in MD simulations. Besides the Ewald summation technique there exist few additional methods, like the fast multipole method and the particleparticle particlemesh method (PPPM, or P3M), which can be used to handle longrange interaction potentials, see e.g. [42]. The results presented here for Coulomb systems have been obtained from simulations using this latter method [43, 44, 45, 46, 47]. In the PPPM scheme the interparticle force is partitioned into (i) a force component that can be calculated on a mesh (the “mesh force”) and (ii) a shortrange (“correction”) force , which is to be applied to closely separated pairs of particles only. In the mesh part of the calculation charged clouds are used instead of pointlike particles and their interaction is calculated on a computational mesh, taking also into account periodic images (for more details see [46, 47]). This way the PPPM method makes it possible to take into account periodic images of the system (in the PM part), without truncating the long range Coulomb or low Yukawa potentials. For screening values the PP part alone provides sufficient accuracy. In these cases the mesh part of the calculation is not used, the interaction forces are summed for particles situated within a (dependent) cutoff radius around the test particle. Identification of these “neighboring” particles is aided by the “chaining mesh technique”.
In the simulations presented here usually a spatially random particle configuration is set up at the initialization, with particle velocities sampled from a Maxwellian distribution of temperature , which corresponds to the desired value of the coupling parameter [see eq. (2)]. The equations of motion of the particles are integrated using the leapfrog scheme or the velocityVerlet scheme. The desired system temperature is reached by rescaling the particle momenta during an initialization phase of the simulation. In equilibrium MD simulations measurements on the system are taken following this phase, in the state of thermodynamic equilibrium. During this phase thermostation is usually no longer applied. If thermostation is necessary, rescaling of particle velocities is to be avoided, algorithms like the NoséHoover thermostat can be applied (see e.g. [38, 48, 49]).
In our studies measurements on the system are taken at constant volume (), particle number () and total energy (). The MD simulations directly provide the pair correlation function (PCF) of the system, which is the basis for the calculation of thermodynamic quantities (not detailed here, see e.g. [32]), and is also required as input to the QLCA equations for the calculation of the dispersion relations and other quantities (see later).
In the MD simulation information about the (thermally excited) collective modes and their dispersion is obtained from the Fourier analysis of the correlation spectra of the density fluctuations
(9) 
yielding the dynamical structure function as [50]:
(10) 
where is the length of data recording period and is the Fourier transform of (9).
Similarly, the spectra of the longitudinal and transverse current fluctuations, and , respectively, can be obtained from Fourier analysis of the microscopic quantities
(11) 
where and are the position and velocity of the th particle. Here we assume that is directed along the axis (the system is isotropic) and accordingly omit the vector notation of the wave number. The way described above for the derivation of the spectra provides information for a series of wave numbers, which are multiples of , where is the edge length of the simulation box. The collective modes are identified as peaks in the fluctuation spectra. The widths of the peaks provide additional information about the lifetimes of the excitations: narrow peaks correspond to longer lifetimes, while broad features are signals for short lived excitations.
3 Theoretical approaches
The Molecular Dynamics calculations compute the dynamical density–density and current–current correlations (dynamical structure functions), from whose behavior the dispersion relations for the collective modes can be inferred. Following the same route in a theoretical analysis would be an extremely ambitious undertaking. Calculating the dynamical structure functions is not an easy task and not much progress has been achieved so far along this line. The singleparticle and collective microscopic dynamics of a classical 3D Yukawa fluid was first analyzed by Barrat et al. [51], on the basis of memory function and modecoupling theories. They have found that the longitudinal current fluctuations and the velocity autocorrelation function cross over continuously from the behavior characteristic of classical fluids with shortrange interactions to the dynamics of a onecomponent plasma as the screening parameter of the Yukawa potential is reduced.
2D Yukawa systems in the liquid phase were considered by Löwen [40] and Murillo and Gericke [52]. In this latter work radial distribution functions have been computed with the hypernetted chain equations and were compared with those obtained from molecular dynamics simulations. The dynamical structure function obtained from the RPA approach extended by local field corrections was shown to be inadequate to reproduce the features of the structure function obtained from molecular dynamics. Ref. [40] focused mostly on the static properties and Brownian dynamics of the system, while also considering some features of the dynamical fluctuations. Applying the viscoelastic approximation Murillo [53] analyzed some aspects of the transverse current fluctuations.
Fortunately, for the determination of the collective mode spectrum a much more direct approach, via the analysis of the dielectric response (tensor) function, is available. Thus the primary goals of the analytical methods discussed below are the determination of the dielectric function and the derivation of the ensuing dispersion relation for the collective modes.
The dielectric tensor in the spatially homogeneous liquid phase is diagonal in the coordinate system, where is along one of the coordinate axes. Accordingly, the collective modes can be classified by their polarization into longitudinal and transverse modes. In the crystalline solid phase the rotational symmetry is broken, the structure of the dielectric tensor is more intricate and the longitudinal and transverse polarizations do not, in general, represent eigenpolarizations anymore. In this Review we are concerned with the collective mode structure of the liquid phase, but the understanding of the behavior of the collective modes in the solid phase has a bearing, as we will discuss, on the formation of the collective modes in the strongly couple liquid phase as well.
3.1 Fluctuation–Disspation Theorem
The link between the , , spectra measured in the simulations and the dielectric function is provided by the Fluctuation–Dissipation Theorem
(12)  
where , is the susceptibility tensor, and is the proper (or total) susceptibility tensor.
At the value where the dispersion relation is satisfied, vanishes. This, in general, happens only at a complex frequency, the imaginary part of which being characteristic of the damping of the mode. Since the dynamical structure functions are plotted and analyzed for real frequencies only, reaches only a minimum at some value of the real , which can be expected to be in the vicinity of the actual complex : this is the frequency value that can be identified at which the peak of the fluctuation spectrum occurs.
3.2 Dielectric Response Function
The tensorial dielectric response function can be expressed either in terms of the susceptibility tensor or the proper (or total) susceptibility tensor and the Fourier transform of the interaction potential (1). This latter depends on the dimensionality of the system. In 3D
(13) 
and in 2D
(14) 
Then
(15) 
In the coordinate system where is along the axis in 3D and along the axis in 2D, the isotropic liquid has the structure
(16)  
and the dispersion relations for the collective modes are given by
(17)  
The longitudinal dielectric function has the immediate physical significance that it relates the externally imposed electric field to the total (external+polarization) field by . In contrast, the transverse dielectric function has welldefined physical meaning only in terms of the full electrodynamics of the system [54]. Here, there is a certain degree of arbitrariness in the definition of . A useful alternative formulation of the dispersion relations is in terms of the external susceptibility
(18)  
This allows expressing the condition for the collective excitation in the universal form
(19) 
embodies all the dynamical properties of the system, which stem partly from interparticle correlations, partly from the random motion of the particles. Over the past half century an immense effort has gone into the calculation of this quantity for Coulomb systems, both classical and quantum. Most of the work focused on weakly coupled () or moderately coupled () systems. Interest in strongly coupled Coulomb and Yukawa systems is more recent [55, 56]. In the strongly coupled domain the dynamics is dominated by correlations. Here and in the sequel we will mostly ignore the effect of thermal motions on ; some comments on how to abandon this simplification will be made later in this Section.
While our focus in this Review is on the strongly coupled liquid phase, it will be instructive and of interest to begin with an orientation based on the weakly coupled Random Phase Approximation (RPA) theory. The RPA or Vlasov description is based on the assumption that the mean field dominates the particleparticle interaction and correlations can be ignored. This is tantamount to taking as that of the noninteracting gas (although perhaps not quite obviously: for a discussion see e.g. [57]): , which, with the neglect of thermal motion is
(20) 
This leads to the simple expressions for the elements of the dielectric tensor
(21)  
where and are the respective 3D plasma frequency and the 2D nominal plasma frequency defined in Eqs. (6) and (7); .
The dispersion relations (and their small approximations) for the 3D and 2D longitudinal modes follow immediately from (17):
(22)  
For the longitudinal mode is acoustic, i.e. , with the 3D and 2D acoustic velocities and :
(23)  
Note that if we compare 3D and 2D systems with the same average interparticle distance, the acoustic speed in 2D is different by a factor than in 3D. The acoustic behavior in 2D is of course at complete variance with the corresponding of an unscreened Coulomb plasma, i. e. the limit , where .
It is clear that there is no mode satisfying the transverse dispersion relation (17b): the mean field RPA model, which is devoid of correlations, cannot support a transverse shear wave, since shear is a fundamentally correlational phenomenon.
3.3 Quasi Localized Charge Approximation
While the RPA provides a description of the weakly coupled gas, the strongly coupled liquid state of a Coulombic or Yukawa system requires a different approach. There have been various attempts over the years to calculate dispersion relations and related quantities for such systems. Noteworthy approaches include the high frequency sum rule method [50], the application of the STLS (Singwi, Tosi, Land and Sjolander) technique originally developed for the electron gas in metals [58, 59], the memory function approach [60, 61] and the viscoelastic model [53, 62, 63, 64].
In the long run, from a practical perspective most of these methods have turned out to be problematic. The problems that occur vary: they range from weak theoretical foundation through being more appropriate for static than dynamical processes, to resulting in an unwieldy formalism. On the other hand, a method originally proposed by Kalman and Golden in [65] that has become known as the Quasilocalized Charge Approximation (QLCA) (for a review see [36, 37]) has led to quite a successful history of accomplishments. The measure of success in this context is (a) the ability to calculate from available static data dynamical quantities that lend themselves to comparison with numerical or laboratory experiments; (b) solid agreement with the outcomes of MD simulations; and (c) a good accord with the newly available laboratory experiments (still in a rather limited number) on complex plasma wave propagation. The following is a concise description of the QLCA method; for more details the reader is referred to [36].
The conceptual basis for the QLCA has been a model that implies the following assumptions about the behavior of a strongly coupled Coulomb or Yukawa liquid: (i) in the potential landscape within the manybody system deep potential minima form that are capable of trapping (caging) charged particles; (ii) a caged charge oscillates with a frequency that is determined both by the local potential and the interaction with the other (caged) particles in their instantaneously frozen positions; (iii) the potential landscape changes slowly to allow the charges to execute a fair number of oscillations; (iv) the escape from the cages of the particles is caused by the gradual disintegration of the caging environment; the time scale of this process is governed by the coupling strength ; (v) the (time and velocity dependent) correlation between a selected pair of particles is well approximated by the (time and velocity independent) equilibrium pair correlation; (vi) the frequency spectrum calculated from the averaged (correlated) distribution of particles represents, in a good approximation, the average over the distribution of frequencies originating from the actual ensemble.
Hypotheses (i)–(iv) have undergone careful testing by a series of MD simulation experiments both for Coulomb and Yukawa systems, and both for 2D and 3D configurations [34, 35], which will be discussed in Section 4. The validity of hypothesis (v) has recently been called into question in relation to multicomponent systems. The short time evolution of the pair correlation function in the vicinity of a particle moving with respect to its environment can certainly be velocity dependent and anisotropic: it is now believed that it is this behavior that is responsible for some discrepancies between MD simulation results and QLCA predictions occurring in binary Coulomb and Yukawa systems. It is not believed, however, that this behavior would be problematic in a single component system. As to item (vi), the question of the dynamical frequency distribution in a liquid has received very little attention, either theoretically or experimentally (the record is better in relation to disordered crystals, where the problem has been posed and approximation schemes have been proposed, although in a language where the central role of the dispersion relation is obscured – see, e.g. [66]). The extension of the QLCA in this direction, while less than pressing, would be desirable.
The central quantity in the QLCA is the dynamical matrix, either in three dimensions () or in two dimensions ():
(24) 
which is formally similar to the eponymous quantity in the harmonic theory of lattice phonons and is derived from the equation of motion of properly constructed collective coordinates. is the dipoledipole interaction potential associated with . is a functional of the equilibrium pair correlation function (PCF) , or of its Fourier transform .
The longitudinal and a transverse elements of the dielectric tensor are now expressed in terms of corresponding elements of :
(25) 
Thus the and local field functions are the respective projections of [36]. is the 3D or 2D longitudinal mode frequency, found in (22). One should keep in mind that in spite of the universality of the expression (24) the explicit forms of the 3D and 2D s are quite different.
We note that the input required in the calculations is the static pair correlation function (PCF). In earlier works PCFs generated by the HNC (hypernetted chain [67]) technique have been used as input data of the QLCA formulae to calculate the dispersion relations. With the advent of computer simulation techniques it turned out to be both more expedient and more accurate to import simulation generated PCFs in the theoretical calculations. The results of the theoretical calculations presented in this paper use PCFs derived from molecular dynamics (MD) computations.
We now can examine the dispersion relations that emerge from (24) and (25) in conjunction with (17a) and (17b). We will consider both the 2D and 3D cases with the corresponding results for the dispersion relations displayed in figures 1 and 2, respectively. Figure 3 will compare the sound velocities and Einstein frequencies for these two cases.
Turning first to the 2D case, the longitudinal dispersion relation becomes
where , and
(27) 
and are Bessel functions of the first kind.
In the mode frequency in (3.3) the RPA solution and the additional correlational part expressed in terms of the pair correlation function are clearly separated. The result can, however, be transformed into an alternate form, expressed entirely in terms of the pair distribution function . By introducing the extended dynamical matrix :
This result shows that the RPA contribution can be interpreted in terms of the same physical model as the QLCA: in this unified formulation the RPA force experienced by the oscillating particle is due to the mean field [] only. Moreover, a reflection on the origin of the “1” term in the integrand identifies it as the generator of the Einstein frequency, the frequency of oscillation of a single particle in the frozen immobile environment of the other particles (see more discussion below):
(29) 
The behavior of the longitudinal mode is still acoustic, but the correlations reduce the acoustic (sound) speed below its RPA value. For the mode frequency approaches the Einstein frequency . This limiting behavior is a remarkable feature of strongly coupled Coulomb and Yukawa liquids [68, 69].
In contrast to the weakly coupled gas described through the RPA, the strongly coupled liquid supports a shear maintained transverse mode. This is reflected in the QLCA through the transverse dispersion relation
where
(31) 
The last step in (3.3) follows from a simple algebraic identity and it reflects the absence of a mean transverse field in the medium.
The longitudinal and transverse dispersion curves for selected and values are displayed in figure 1. The behavior of the transverse mode is also acoustic, but the correlation maintained acoustic speed is substantially below its longitudinal counterpart. However, the result that the transverse mode extends all the way to is spurious: the liquid is unable to support a shear wave in the uniform limit. The reason for this flaw is well understood: it has to be sought in the neglect of the migrationaldiffusional damping. The introduction of a phenomenological collision frequency [70] or of a semiphenomenological extension of the QLCA (studied sofar only for the Coulomb case  see below) [68] provide an acceptable remedy.
For the transverse mode also approaches the same Einstein frequency as the longitudinal mode, as dictated by the isotropy of the liquid.
Also shown are in figure 3 the and dependences of the acoustic speeds and of the Einstein frequency. For moderate values the sound velocities can also be obtained from the semianalytic formulae [71]:
(32) 
where is the correlation energy per particle.
Turning now to the 3D case, the formal results of the previous derivation can, mutatis mutandis, be taken over, with the understanding that the explicit forms of the RPA frequency and the kernels and are different from their 2D counterparts.
(33) 
(34)  
and
(35) 
There are now 2 degenerate transverse modes.
The expression for the Einstein frequency is also modified:
(36) 
The 3D longitudinal and transverse dispersion curves for selected and values are displayed in figure 2. The and dependences of the 3D acoustic speeds and of the Einstein frequency are shown together with their 2D counterparts in figure 3. Finally, for moderate values the sound velocities can again be obtained from the semianalytic formulae [70, 72]
(37) 
and
(38) 
The results of a comparison between the dispersion properties of the collective modes in the 2D and 3D systems (assuming the same interparticle distance) may be summarized as follows.

While for finite values the qualitative behaviors of the two systems are very much the same, there is the well known fundamental difference in the Coulomb limit between the 2D and 3D systems as to the small dispersion of the longitudinal mode: for 2D, but in the 3D case , the 3D plasma frequency.

Not unrelated to this difference is the behavior of the longitudinal acoustic speeds at finite values: since in 2D and in 3D , for small the latter exceeds the former by the factor .

In contrast, the transverse acoustic speeds exhibit only a mild dependence and it is the 2D speed that is slightly higher than its 3D counterpart.

A similar 2D dominance prevails for the respective Einstein frequencies that govern the behavior of the modes: in 2D the Einstein frequency assumes, for any , a somewhat higher value than in 3D.
While (i) and (ii) are effects originating from the basic difference caused by the long range behavior of the Coulomb potential in a 2D vs. a 3D geometry and are already reflected in the RPA description, (iii) and (iv) are correlational phenomena and they point at the more important role the correlations play in 2D than in 3D.
As a closing comment, it should be reemphasized that the QLCA ignores possible damping mechanisms and Doppler shift, phase mixing, etc., due to the migrationaldiffusive motion of the particles and of velocity dispersion. A method has been recently proposed [73, 74, 75, 76] for the extension of the QLCA to take some of the neglected effects into account by combining the , as local field factors with the Vlasov densitydensity response. In this approximation
(39) 
where is the longitudinal (transverse) Vlasov density response function of noninteracting particles. The application of this formalism to Yukawa systems has not been done yet, but in an early work [68] it was shown that in a 2D Coulomb system the combined effect of phase mixing and Landaudamping leads to the elimination of oscillations in the dispersion curve. The effect of Landau damping, which is not expected to play a major role at high values, is probably overestimated in this work.
3.4 Lattice phonons
With the caveat that the present Review addresses primarily the strongly coupled liquid state, it will be still useful to provide an overview of the phonon dispersion in a 2D or 3D Yukawa crystal. Such an overview will help to understand the structure of the liquid state in terms of a model resembling a disordered lattice and to view the collective modes in the liquid as being akin to the phonon excitations in the lattice.
The phonon dispersion is traditionally calculated in terms of the lattice dynamical matrix defined as
(40) 
with a summation over all the lattice sites , keeping fixed . The resemblance to the extended QLCA dynamical matrix is not accidental. In contrast to the QLCA equivalent, however, the lattice dynamical matrix reflects the symmetry of the underlying lattice and not the rotational invariance of the isotropic liquid. Nevertheless, a dielectric tensor can be constructed along the same line in terms of the matrix , which is defined now as with its mean field contribution removed
(41) 
This leads to a structure analogous to (25):
(42) 
The diagonalization of (or of ) is now possible in the coordinate system of the eigenvectors, whose orientations, in general, do not coincide either with the direction of or with the crystallographic axes.
To find the eigenmodes one can follow the traditional method (see e. g. [77, 78]) of solving the secular equation
(43) 
or continue to follow the path of working with the dielectric tensor. This latter approach ensures that continuity with the liquid and RPA formalism is maintained. The dispersion relation in terms of becomes
(44) 
an obvious generalization of (17a). In fact, except in the degenerate isotropic case, it also includes the transverse relation (17b).
The 2D Yukawa system crystallizes in a triangular (hexagonal) lattice. The phonon spectrum was first calculated by Peeters and Wu [80], followed by Wang et al. [81]; a definitive calculations of the dispersion and the polarization for all propagation angles are given in [79, 82]. These results are shown in figure 4. The mode polarizations are purely longitudinal or transverse for propagation along the crystallographic axes ( and ) only, otherwise they are mixed as shown in the figure. is the propagation angle measured from the axis pointing towards the nearest neighbor. The angle indicated in figure 4 is the polarization angle measured with respect to the propagation vector . The dispersion curves are periodic in , but the period is simply the reciprocal lattice constant only along and ; for intermediate angles it is much longer, given by the formula
(45) 
where and are minimal integers satisfying
(46) 
The dispersion curves of the lattice and those of the strongly coupled liquid do not show much resemblance. Yet, if one views the liquid as an aggregate of locally ordered domains whose symmetry axes are randomly distributed, then the similarity to the liquid dispersion should be sought in a suitably averaged dispersion of the lattice. The strong angular dependence of the period suggests that an angular average should generate through phase mixing a smooth dispersion. This was carried out by projecting out the longitudinal and transverse components of the eigenmodes and comparing their respective angular averages with the longitudinal and transverse liquid modes [82]. Figure 5 shows that the agreement is quite reassuring. In principle, of course, one has to distinguish between the spectrum of an average of configurations and the average of the spectra of each of the configurations: that this observation notwithstanding the similarity persists can be taken as an indication that the sizes of the ordered domains in the liquid state are sufficiently large to diminish the effect of interaction between the domains.
The 3D Yukawa system crystallizes in a bcc or a fcc lattice (depending on the value of ). A phase diagram has been given by Hamaguchi et al. [28]. Due to the existence of 3 rather than 2 eigenmodes and to their dependence both on the azimuthal and the polar angles of propagation a much more complex phonon spectrum is expected than in 2D. So far no systematic published study of this spectrum seems to exist; in an unpublished work, however, Sullivan, Kalman and Kyrkos [83] have generated a series of dispersion and polarization diagrams. A sample of these, for a number of and angles, including propagations along the principal crystallographic axes is given in figure 6. ( is the polarization angle measured with respect to the propagation vector , is the polar angle in the plane and is the azimuthal angle measured from the axis.) Since no averaging has been performed, their comparison with the liquid spectrum at the present time is difficult.
3.5 Einstein frequencies
In addition to the collective excitations, Einstein frequencies represent a dynamical manifestation of the strong interaction in Yukawa systems. Einstein frequencies, as noted above, are the frequencies of oscillation of a single particle of the system (the “test particle”) around its equilibrium position in the immobilized frozen environment of the other particles of the system. For obvious reasons, from the experimental point of view the “freezing” of the system but one particle is not a realistic proposition.
Thus until the advent of dusty plasma experiments Einstein frequencies were considered more of a theoretical construct than an observable quantity. The realization, however, that in the strongly coupled liquid state (but not in the crystalline solid) they represent the asymptotic limit of the mode dispersion has promoted the Einstein frequency to the rank of observable quantities [84, 85].
In the crystalline solid state, where the test particle occupies a lattice site, the assumption that the potential experienced by the test particle is a quadratic function of the coordinates with a positive definite second derivative is in accord with the basic model of the harmonic theory of phonons. The maximum number of eigenfrequencies of oscillation is equal to the dimensionality of the system ( or ); because of the lattice symmetry induced degeneracy the actual number may be less than . In a disordered lattice the degeneracy is removed and the frequencies depend on the actual realization of the disorder. In this case one has to distinguish between the “microscopic” Einstein frequencies () each of which is generated by a particular realization of the disorder and characterized by a distribution over the ensemble, and their ensemble average . It is this latter that will be continued to be referred to as “Einstein frequency” in the rest of this paper.
In addition, it is useful to consider the quantity
(47) 
i.e. the sum of the squared eigenfrequencies in a particular realization. Obviously, but the distributions of and can be quite different.
In a strongly coupled liquid the very notion of “equilibrium position” is questionable. Nevertheless, the “quasilocalization” condition, the basic tenet of the QLCA, is well satisfied for high values, as demonstrated by MD experiments [34] the details of which will be discussed below. It is in this sense that the notion of the Einstein frequency and its distribution can be extended to the case of the strongly coupled liquid.
In a 3D Coulomb crystal the Einstein frequency is determined solely by the background, unaffected by the distribution of the (frozen) particles. This is the consequence of Gauss Theorem which, in turn, follows from the Poisson Equation that the 3D Coulomb potential satisfies. In this case
(48) 
In a disordered lattice or in a liquid (48) is not valid anymore; it is replaced by the weaker statement
(49) 
the socalled Kohn Sum Rule [86], which also follows from the Poisson Equation.
Thus, in a disordered system while has a spread, does not. As to the average, (49) is of course also tantamount to
(50) 
For a genuine Yukawa potential the situation is quite different. The Yukawa potential satisfies the screened Poisson Equation rather than the Poisson Equation. A useful statement can be made now only for , which now can be expressed in terms the average of the Yukawa potential as experienced by the test particle at [33]:
(3.5) is in agreement with (36), the result obtained from the QLCA. The third line clearly shows that, remarkably, in the Coulomb limit the Yukawa Einstein frequency reduces to the background induced (50), even though the Yukawa system exists without any background. It can also be noted that is the interaction energy density of the system (with [positive] Hartree plus [negative] correlation contributions). Since the energy is the lowest in the ordered state, the Einstein frequency must increase with increasing disorder. According to the known phase diagram of the 3D Yukawa system [28] – as already mentioned – the system crystallizes in a bcc or a fcc lattice. The corresponding Einstein frequencies [87]
(52)  
constitute an absolute lower bound.
In the 2D Coulomb system Gauss Theorem does not apply, the background plays no role and neither the Poisson Equation nor its screened variant is satisfied. Consequently, the Einstein frequency is determined by the distribution of the surrounding particles, both for Yukawa and Coulomb systems. In general
(53) 
in agreement with the QLCA result (29).
An argument similar to the one discussed in relation to the 3D case leads to the conclusion that here also the ordered state exhibits the lowest Einstein frequency. The lattice structure is now hexagonal, for which
(54)  
These values constitute then the lowest bound for the 2D Einstein frequencies.
In addition to the frequencies, the Einstein oscillations are also characterized by their eigenpolarizations. It is the distribution of the polarization angles which is of interest; this question has been investigated, however, only for the 2D case [82]. In the perfect hexagonal lattice the degeneracy of the eigenmodes renders this distribution isotropic. It is also isotropic in the liquid phase. However, in the intermediate range where the lattice disorder develops the degeneracy for the microscopic eigenmodes is removed and the rotational invariance of the distribution is reduced to the sixfold symmetry of the underlying lattice. More will be shown about this remarkable effect in Section 4.2.
4 Simulation results
In this Section we review the results of the extensive MD simulation work carried out since the beginning of this decade on the dynamical properties of Yukawa liquids. Most of the work was motivated by the QLCA theory and accordingly a great portion of the results pertaining to the areas where QLCA predictions are available are accompanied by comparisons with the theoretical predictions. However, the information generated by the simulations goes well beyond those predictions: this is eminently true for the frequency spectra of the dynamical densitydensity and currentcurrent correlation functions [dynamical structure functions , , ]. Beyond predictions pertaining to the peak positions of the spectra, identified as the frequencies of the collective excitations, the QLCA does not provide, apart from some qualitative estimates, any basis for comparison in this respect. While other works, based mostly on the memory function formalism [51, 80, 81, 88], have presented theoretical descriptions of some of the features of the structure functions, we have made no attempt to relate to these, rather scant, results for the purpose of comparison with simulations.
As noted in Section 3.3, the basic hypotheses (i)–(iv) of the QLCA theory have undergone careful testing by a series of MD simulation experiments both for Coulomb and Yukawa systems and both for 2D and 3D configurations [34, 35]. With increasing values, a visual inspection of the potential landscape clearly indicates the formation of potential wells [35]. Examination of the phase space trajectories reveals a clear morphological difference between low and high situations: in the first case the trajectories are open, interrupted by propagating oscillatory portions, while in the second case the trajectories are mostly closed and exhibit a loop structure characteristic of localized oscillatory motion [34]. An example of this behavior is illustrated for a 3D Coulomb liquid in figures 7(a) and (b) for = 2.5 and = 160, respectively.
The quantification of the relationship between localization and the strength of the coupling has been carried out by invoking a technique due to Rabani et al. [89]. Here a “cage correlation function” was introduced to characterize the gradual disintegration of the cage of the nearest neighbors and the escape of the caged particle. The main results shown in figures 7 (c) and (d) for 3D and 2D Coulomb and Yukawa systems illustrate the duration (in terms of plasma oscillation cycles) of the caging (decorrelation time, ) as a function of and . In the case of the 3D system, at = 0 and = 160 the cages decorrelate during 50 plasma oscillation cycles. The decorrelation time is reduced to a single cycle at 7. In the case of the 2D system it takes about 100 cycles for the cages to decorrelate at = 0 and = 120, and we reach = 1 at 2.5. In the high domain we observe a strong dependence of the decorrelation time on , both in 3D and 2D systems. At low values of , however, depends only slightly on . The decrease of the decorrelation time for increasing can be compensated by increasing , as it can be seen in figure 7(c) and (d) [35]. It is noted that the data shown in figure 7 convey information about the “average behavior” of the particles, it is however, recognized [35] that the surrounding of individual particles may change in a different way, due to e.g. avalanche type excitation and migration [90]. Finally we note that the caging of the particles at high values determines many of the systems properties as it has been discussed by Daligault for 3D Coulomb liquids [91].
4.1 Threedimensional Yukawa liquids
The first molecular dynamics simulations on the wave dispersion relations in the fluid phase of 3D Yukawa systems were reported by Hamaguchi and Ohta [92, 93]. Their results confirmed the earlier theoretical predictions of Rosenberg and Kalman [72] on the longitudinal wave dispersion and were mostly in agreement with the simultaneously published full QLCA calculations of Kalman et al. [70]. They also demonstrated that the transverse wave dispersion has a cutoff at a long wavelength even in the case of weak screening.
This work was followed by a series of MD simulation for the collective excitations in 3D Yukawa liquids to provide further comparison with the predictions of the QLCA theory. The simulations – of which the results are presented here for the first time – have been carried out using = 12 800 – 15 625 particles.
To illustrate qualitatively the features of the behavior of the collectives excitations the spectral decomposition of the longitudinal and transverse current fluctuations is plotted in figure 8 for 3dimensional Coulomb and Yukawa liquids. In the case of the Coulomb plasma, at low wave numbers the frequency of the longitudinal () mode is concentrated within a narrow frequency range [see figure 8(a)] near the plasma frequency. With increasing wave number the frequency of the mode gradually spreads over a wider domain and shows a slightly decreasing tendency. In sharp contrast with this behavior the mode frequency is spread over a wide domain, as illustrated in figure 8(b). The mode of the Yukawa system is quite different from that in the Coulomb case, the wave frequency approaches zero at wave number. The frequency increases with increasing wave number up to about , and then starts to decrease slightly. Meanwhile the frequency distribution gets gradually wider. The mode in the Yukawa case appears to be similar to the corresponding mode in the Coulomb system, although the frequency is lower, due to the weaker interaction of the particles, as a consequence of the screened potential.
For a better quantitative analysis representative dynamical structure functions (density fluctuation spectra) and spectra of the longitudinal and transverse current fluctuations, and , are plotted in figures 9 and 10, respectively, for a high and a medium case. The obtained for the Coulomb case ( = 160, = 0, see figure 9(a)) peaks at nearly the same frequency for the different values of the wave numbers plotted, which are multiples of (determined by the size of the simulation box). In the presence of screening (Yukawa potential), as shown in figure 9(d), the behavior of changes significantly: at 0 the wave frequency 0 [ is defined by (6)]. The contrast between the = 0 and the 0 cases is also well seen in figure 11(a), where the dispersion curves derived from the fluctuation spectra are displayed. The dispersion curves for 0 are quasiacoustic (), with a linear portion near = 0, which gradually extends when is increased. The (,) pairs for which the dispersion graphs are plotted in figure 11 have been selected to represent a constant effective coupling = 160. This definition of relies on the constancy of the first peak amplitude of the pair correlation function , similarly to the case of 2D Yukawa liquids [32].
Peaks in the spectra of the compressional mode [plotted in panels (b) and (e) of figures 9 and 10] appear at the same frequency as those in the corresponding functions, as these functions are linked via the relation
(55) 
Compared to those characterizing the mode, peaks in the mode spectra are rather broad, as it can be seen in panels (c) and (f) of figures 9 and 10. In the case of this mode there is no dramatic change between the behavior when changes from zero to a nonzero value, only the mode frequency decreases, as can be observed in figure 11(b).
Comparison of the dispersion relations obtained from the MD simulations [via ] and QLCA calculations [see equations (33)–(35)] is presented in figure 11. Here, in the calculations of the QLCA results, we have made use of the functions obtained from the MD simulation. The agreement between the two sets of data is excellent for the mode, while some difference in the frequency of the waves can be seen in figure 11(b). This latter may originate from the inaccurate determination of the peak positions of the rather broad spectra. It should be noted though that while the theoretical calculations provide an oscillatory dispersion curve for (see figure 2), simulations provide reliable results (for collective excitations) for typical liquidphase conditions for . (At higher values the thermal contribution in apparently masks the collective mode peak.). The simulation results here resemble the measured 2D dispersion curves in the liquid phase [15]. Another difference is the cutoff of the mode dispersion curve at finite wave numbers. This disappearance of the shear modes for is a well known feature of the liquid state [50, 95, 96], and the sharp cutoff for a finite has also been observed in simulations of Yukawa systems [53, 93]. It has been already noted that this cutoff is not accounted for by the QLCA, as it does not include damping effects.
The sound velocities, derived in (37), are plotted in figure 12(a), while figure 12(b) displays the Einstein frequency, which is defined in (36). In the limit (36) gives , and decreases with increasing . For comparison, the Einstein frequency data of Ohta and Hamaguchi for a fcc lattice [94] are also plotted in figure 12(b). We find an excellent agreement between the two sets of data.
Numerical experiments were also performed to determine the distribution of the microscopic Einstein frequency . To accomplish this, frequency histograms based on a few hundred, temporally uncorrelated particle configurations have been constructed. For the raw (particle position) data the harmonic matrix for every particle has to be generated:
(56) 
where is the equilibrium position of the th particle (local minimum of the potential surface), is the interaction potential, and represent the Cartesian coordinates. The eigenvalues of are the squared Einstein frequencies (3 for every particle), while the eigenvectors provide the polarization of the oscillation.
A series of frequency histograms for an effective coupling parameter and different values of are shown in figures 13(a)(d). The frequency distributions exhibit three peaks (although this is less visible in the case). With increasing screening the distribution of frequencies becomes wider and its mean value is shifted towards lower frequency. The QLCA results for the Einstein frequency [obtained from Eq.(36) using pair correlation functions generated in the MD simulation], corresponding to the different values of are also indicated in figures 13(a)(d). The values are in good agreement with the simulation results. The effect of at fixed () screening is illustrated in figure 13(a). A six time decrease of the coupling parameter results in approximately doubled width of the Einstein frequency distribution.
Figure 13(e) shows the histograms for , sums of the 3 microscopic squared Einstein frequencies, for different values of : there is a qualitative difference between the Coulomb case where there is only a single frequency (a narrow peak) and the cases, where a distribution of frequencies is apparent. The reason for this difference has been discussed in Section 3.
Further information on the collective behavior is contained in the velocity autocorrelation function (VACF)
(57) 
where the average is taken over the particles and different initial times.
The behavior of the velocity autocorrelation functions of 3D Yukawa liquids obtained at several values of the and parameters is illustrated in figure 14. Analyzing the behavior of at constant , we find a transition from monotonically decreasing into an oscillating type when is increased [figure 14(a)]. Similarly, the shape of changes drastically when is varied at constant , as shown in figure 14(b). For more detailed analysis see [94].
Using the Einstein frequency for the normalization of time, instead of the plasma frequency (as in figure 14) the functions belonging to the same = 160 for a series of values are displayed in figure 15(a). Using this normalization of the timescale the curves exhibit a nearly universal behavior, where at least the first few peaks of nearly overlap. This observation emphasizes the importance of the Einstein frequency in the dynamical behavior of the system.