Dispersion and damping of two-dimensional dust acoustic waves: Theory and Simulation

A two-dimensional generalized hydrodynamics (GH) model is developed to study the full spectrum of both longitudinal and transverse dust acoustic waves (DAW) in strongly coupled complex (dusty) plasmas, with memory-function-formalism being implemented to enforce high-frequency sum rules. Results are compared with earlier theories (such as quasi-localized charge approximation and its extended version) and with a self-consistent Brownian dynamics simulation. It is found that the GH approach provides good account, not only for dispersion relations, but also for damping rates of the DAW modes in a wide range of coupling strengths, an issue hitherto not fully addressed for dusty plasmas.


Introduction
Dusty, or complex plasmas have grown into a mature research field with a surprisingly broad range of interdisciplinary facets [1, 2,3]. There has been significant interest in the study of dynamics of collective excitations of dusty plasmas, both in experiments [4,5,6,7,8,9,10] and through theory/numerical simulation of strongly-coupled Yukawa systems [11,12,13,14,15,16,17,18,19,20,21,22,23]. Of particular interest in these studies is theoretical modeling of longitudinal and transverse dust acoustic wave (DAW) modes [11], and consequently several assessments of how the existing analytical theories compare with the results from computer experiments have appeared recently [19,20,22,23,24]. Specifically, by comparing the Brownian dynamics (BD) simulation results for the DAW spectra with various theoretical models, it was shown that the so-called Quasi-localized charge approximation (QLCA) [15,20] provides an overall good account of the DAW dispersion relations, i.e., the peak locations in those spectra, for both three-dimensional (3D) and twodimensional (2D) dusty plasma states characterized by the coupling strength Γ in the range 10 < Γ < 1000 [22,24]. However, since the QLCA provides no information about the profile of such spectra and, consequently, gives no account of damping processes for the collective modes [15], it is desirable to explore alternate theoretical approaches to strongly coupled systems. In that context, generalized hydrodynamics (GH) offers possibility to describe full wave spectra by introducing viscoelastic effects in the dynamics of dusty plasmas in a phenomenological manner, as was demonstrated in studying DAWs in 3D dusty plasmas [12,13,14].
In contrast to ordinary hydrodynamics (OH), which is only valid for weakly coupled fluids with Γ ≪ 1 in the long wavelength limit, GH covers wider ranges of coupling strengths and wavelengths, where structural effects begin to show up. Particularly suitable and physically intuitive framework for implementing the GH model is provided by descriptions of the Lennard-Jones fluids, pioneered by Ailawadi et al. [25], Boon and Yip [26], and Hansen and McDonald [27]. This framework may be applied to a dusty plasma by asserting that dust particles form a quasi-neutral fluid with the inter-particle interactions described via a Debye-Hückel, or Yukawa potential, which arises due to screening of the charge accumulated on dust particles by the background electron and ion fluids. Successful applications of the GH model to dusty plasma were conducted by Kaw and Sen [13] and Murillo [14] to study the wave dispersion and the shear wave cutoff in 3D dust liquids, respectively.
However, there is still great demand for a GH model for 2D strongly coupled dusty plasmas (SCDPs), especially because such systems have become particularly favored in recent laboratory experiments [3,28]. Therefore, our principal motivation in this work is to develop and implement a GH model for studying the collective dynamics in 2D SCDPs. In addition, it is necessary to generalize the GH model to include the effects of collisions of dust particles with the neutral gas in the background plasma in a manner similar to that used for colloidal suspensions in describing the collisions of macroions with the solvent molecules [29]. While such collisions give rise to the classical Brownian motion of macroparticles in both these systems, it should be stressed that the friction forces on dust particles due to the neutrals in dusty plasmas are generally much weaker than the friction forces on the colloidal macroions. Accordingly, effects of the neutral gas on the DAW dispersion relations are found to be negligible for the range of parameters of interest in experiments with dusty plasmas, thus validating the use of molecular dynamics (MD) simulations for such systems [24]. However, it is still unclear as to what extent can damping due to the neutral gas compete with the viscous damping of longitudinal DAWs at finite frequencies and finite wavelengths in dusty plasmas [30]. While this difficult issue naturally arises in this work, our primary motivation in using BD simulation lies in the fact that treating the dust particles as Brownian particles provides a natural and convenient way to eliminate the need for thermostation that arises in the MD simulation of dusty plasmas [22,31].
In this paper we perform a BD simulation of strongly coupled 2D Yukawa liquids, and use the results to evaluate the equilibrium radial distribution function (RDF) and the static structure factor, as well as the (power) spectral densities for both the longitudinal and transverse current densities in such systems. The former two quantities enable us to use the GH approach within the memory function formalism to evaluate both the dispersion relation and the damping rate of the DAW modes in 2D dusty plasmas by enforcing the low-order, high-frequency sum rules upon the theoretical spectral densities. It is found that the GH approach provides dispersion relations that compare well with those resulting from the simulation spectra over broad ranges of wavelengths and coupling strengths. Additional comparison with the dispersion relations from the QLCA model shows that the GH approach provides a good account of the direct thermal effect at lower coupling strengths and shorter wavelengths, which is absent in the QLCA approach, but is seen in the simulation data. On the other hand, a simple extension of the QLCA dispersion relation, which was recently proposed in Ref. [22], is found to be in a much better agreement with both the GH results and simulation data. Most importantly, the GH results are shown to yield a good fit to the spectral density profiles from our BD simulations, thus providing semi-analytical modeling of the wave-number dependent damping rates of the DAW modes in strongly coupled dusty plasmas in a broad range of coupling strengths, an issue that has not been fully addressed so far [22]. Finally, we also find that the effect of collisions with the neutrals on the DAW damping rates is well accounted for by introducing a local in time friction force into the GH equations.
The manuscript is organized as follows. The details of the BD simulation are given in Section 2. Theoretical development of the generalized hydrodynamics is discussed in Section 3. Results based on two simple models for the memory function are presented in Section 4, while we conclude in Section 5.

Simulation
We perform BD simulation of a 2D system consisting of N = 4000 dust particles, each carrying a constant charge q d , which are initially placed at random positions within a square simulation cell in the r = (x, y) plane with periodic boundary conditions. Evolution of the system is governed by equations that may be regarded as a stochastic generalization of the usual MD method, where the effect of collisions with the neutral particles in the background plasma is modeled by a Langevin generalization of the usual Newton's equations. Therefore, equations for the velocity and the position vectors of the i-th dust particle are given by where m d is the mass of each dust particle and F(r) = −(r/r)dU(r)/dr is the pair-wise force between two dust particles a distance r = x 2 + y 2 apart, interacting via the Yukawa potential, U(r) = (q 2 d /r) exp (−r/λ D ), with λ D being the Debye screening length of electrons and ions in the background plasma. Collisions of the i-th dust particle with the neutrals in the plasma give rise to a systematic drag force, which we model by the Epstein drag coefficient γ n , and to a cumulative random force, which we model by a delta-correlated Gaussian white noise. When the system is in thermal equilibrium, the friction coefficient γ n and the stochastic acceleration of the i-th particle A i (t) are related to the ambient temperature T d via the fluctuation-dissipation theorem. Equations (1) are solved for our system by the Gear-like predictor-corrector algorithm for BD simulation [31], which was previously used in modeling both the equilibrium structure [32] and collective modes in 2D dusty plasmas [22].
Equations (1) can be expressed in terms of just three parameters [22]: the coupling strength Γ = q 2 d / (ak B T d ), the reduced neutral drag coefficient γ = γ n /ω pd , and the screening parameter κ = a/λ D , where a ≡ 1/ √ πρ is the average inter-particle separation with ρ being the average surface density of dust particles, and ω pd = 2q 2 d / (m d a 3 ) is the characteristic dust-plasma frequency. In this work, we perform BD simulations using γ = 0.06 (which is a value chosen as a matter of convenience only) and κ = 1 as standard values, in a range of coupling strengths 20 ≤ Γ ≤ 1000 covering liquid-to-crystalline states of dusty plasma [22].
Further, the current density of dust particles is defined as with the Fourier transform of its cartesian component α given by Assuming that the dust collective modes propagate along the x-direction, i.e., k = (k, 0), we let α = x or y, allowing us to define the corresponding longitudinal or transverse current density auto-correlation functions as where the angular brackets denote ensemble averaging over the initial time. The longitudinal/transverse spectral densities are then obtained as the Fourier transform of the respective current density auto-correlation functions as, The Fourier transforms defined in Eq. (5) are evaluated using Fast Fourier transform from simulation data. An equivalent and useful definition of the spectral densities, to be used later, is given as the real part of the Laplace transform [27] of the current density auto-correlation functions,

Generalized Hydrodynamics
The basic approach we will take in the development of a GH model for a 2D dusty plasma is to treat the dust system as a compressible, viscous neutral fluid of particles interacting via the Yukawa potential, analogous to the way how molecules in a Lennard-Jones fluid interact [27,30]. In that context, we note that the OH approach provides a satisfactory account of the long-time (or low-frequency) response of a fluid in the weak coupling regime by means of the familiar Naiver-Stokes equation, in which viscous effects are treated as instantaneous (or local in time) internal friction force, related to fast thermalization of dust particles due to their mutual collisions. However, as the coupling strength increases, the effect of "caging" sets in, so that individual dust particles are temporarily trapped in potential wells that migrate through the system on a slow time scale, similar to the picture invoked in the QLCA model. Therefore, the short-time (or high-frequency) response of the system is dominated by elastic effects due to the restoring forces on dust particles in the itinerant potential wells. In the GH approach to the Naiver-Stokes equation, a transition from the purely viscous, fluid-like behavior to the elastic response of a solid-like system is described by postulating a non-local (in time) friction with a memory function characterized by relaxation time τ , such that the limit of a viscous fluid is recovered at times t ≫ τ , whereas the solid-like elastic effects are dominant at times t ≪ τ . In addition to the collisions among dust particles, it is necessary to include in the GH approach also the effect of their repeated collisions with the neutral molecules in the background plasma. These collisions may be treated as a Gaussian white noise, giving rise to a picture of Brownian motion, where each dust particle is subject to a local (in time) friction force with the neutral drag coefficient γ n . We note that for typical dusty plasmas one expects γ n τ ≪ 1. On the other hand, at times t ≪ γ −1 n , the Brownian motion of dust particles due to the neutrals may be considered as ballistic. Therefore, it is justified to ignore the effect of neutral drag at the time scale t ∼ τ ≪ γ −1 n where viscoelastic effects dominate due to interactions among the dust particles [29]. On the other hand, at long times, such that t ∼ γ −1 n ≫ τ , both the viscous drag and the neutral drag on dust particles may be treated by means of two additive, local in time friction forces in a generalized Navier-Stokes equation because the underlying physical mechanisms of energy dissipation are statistically independent.
Therefore, we first ignore the neutral drag and introduce a memory function formalism, which incorporates the effects of viscous relaxation into the short-time dynamics of the dust collective modes by enforcing high-frequency sum rules upon the power spectral densities of such modes [26]. The even-order sum rules are defined in terms of the corresponding frequency moments as 1 2π whereas all odd moments of frequency vanish since P l,t (k, ω) are even functions of ω. For example, the zeroth-order sum rule is given by [26] 1 2π where v th = k B T d /m d is the dust thermal speed. A connection with microscopic properties of the system is accomplished by expressing, e.g., the second moments of frequency for the longitudinal and transverse collective modes in Eq. (7a) with n = 1 in terms of the RDF of the dust layer, g(r), respectively as [26] Once the initial-value properties of the relevant memory functions are fixed by enforcing the sum rules in Eqs. (7b) and (7a), one may incorporate the effect of the neutral drag by simply adding a local dissipative force into the fluid equations of motion [13,30], as explained in section 3.2.

Collective Modes
The starting point for the extension of OH to GH is given by the linearized continuity equation and the Navier-Stokes equation, which may be written for the longitudinal and transverse dust current densities as [26] ∂ ∂t where χ T is the isothermal compressibility, ν l = 2ν 1 + ν 2 is defined as the longitudinal viscosity, with ν 1 being related to the shear viscosity and ν 2 to the bulk viscosity, ρ is the average surface density, δρ is a small perturbation in density, and j l,t (k, t) is a small perturbation to the longitudinal/transverse current density, assumed to be zero on average. Inserting the integral of Eq. (10a) into Eqs. (10b) and (10c), one can obtain integrodifferential equations for the current density auto-correlation functions C l,t (k, t) defined in Eq. (4), which must be subject to the initial condition C l,t (k, 0) = v 2 th . In order to model the short-time correlation effects, the resulting equations can be modified to satisfy the frequency sum rules in the following way [25,26]. In the first step, longitudinal spectral density is forced to satisfy the zeroth frequency sum rule, Eq. (7b), by replacing the isothermal compressibility with the static structure factor, S(k), as which generalizes a standard statistical-mechanical relation in the limit k → 0. In the next step, generalizations of the longitudinal and shear viscosities are introduced via the as yet undefined memory functions, φ l (k, t) and K t (k, t), respectively, giving integro-differential equations of the form where It can be shown that the second frequency sum rule, Eq. (7a) with n = 1, can be satisfied if the initial values of the longitudinal and transverse viscosity memory functions are expressed in terms of the corresponding second frequency moments, Eqs. (8) and (9), respectively, as [26]

Effects of neutral drag
For typical conditions in 2D dusty plasmas, the effect of neutral drag on the dispersion relations of DAW modes is expected to show up at very low frequencies only, ω < γ n , and we have indeed found in our BD simulations that the actual value of the drag coefficient γ n has very little effect on the dispersion relations in comparison to the viscoelastic effects. However, the spectral widths in Eqs. (5) or (6) are found to be affected by the neutral drag, as discussed in section 4. With the short time dynamics of both the longitudinal and transverse dust collective modes fixed by enforcing the sum rules upon the initial-value properties of the relevant memory functions, we may now follow suggestions of Refs. [13,30] and introduce the long time relaxation effect due to neutral drag by simply adding a local dissipative term in Eq. (12), We proceed with solving Eq. (15) by means of Laplace transform in order to derive expressions for the spectral density by invoking the definition Eq. (6). For the longitudinal mode we obtain where φ ′ l (k, ω) and φ ′′ l (k, ω) are the real and imaginary parts, respectively, of the Laplace transformed longitudinal viscosity memory function, L[φ l (k, t)] s=iω , and where we have explicitly defined Similarly, we obtain for the transverse mode with K ′ t (k, ω) and K ′′ t (k, ω) being the real and imaginary parts, respectively, of the Laplace transformed transverse viscosity memory function,

Model Memory Functions
We note that no approximations were used in the development of the GH approach so far. Specific forms of Eqs. (16) and (18) may now be obtained by introducing simple phenomenological models for the corresponding memory functions.

Exponential model Choosing an exponential longitudinal viscosity memory function
to be used in Eq. (16). We first discuss the limits of very short and very long viscous relaxation times τ l . In the limit of small relaxation time, ωτ l ≪ 1, we recover the so-called delta function model for the memory function be letting φ l (k, 0) = ν l /τ l , where ν l is the longitudinal viscosity. This situation corresponds to a fluid where viscous drag is described by a local friction force, so that Eq. (16) is reduced to [26] giving the dispersion relation as ω = ω 0 (k) with ω 0 (k) defined in Eq. (17), and having the full width at half maximum of γ n + k 2 ν l . This result provides a relatively simple account of the interplay between the neutral drag and viscous damping, assuming that the quasi-static longitudinal viscosity ν l can be defined properly [33].
On the other hand, it is important to study the opposite limit, ωτ l ≫ 1, for the sake of comparison with the QLCA model, which is inherently valid in this regime [35]. In particular, we find that the viscous damping vanishes in this limit and the elastic effects in the system's response prevail, giving the dispersion relation ω = ω l with the second frequency moment given in terms of the RDF via Eq. (8). We note that this situation is analogous to that in the QLCA model with two important differences: the dispersion relation ω = ω l ∞ (k) is different from that found in the QLCA model, and the neutral drag still provides a mechanism for damping of the DAW modes.
For finite values of τ l , a dispersion relation in the exponential model for the longitudinal viscosity memory function may be obtained from the peak positions in the spectral density given by Eq. (16) with Eqs. (19a) and (19b). This spectral density is fully determined by the functions S(k) and ω 2 l (k) via Eqs. (17) and (14a), respectively, and by the longitudinal viscosity relaxation time τ l , which may be a function of k as well. We note that the former two functions can be calculated, at least in principle, from the RDF of the dust layer by using the usual definition for S(k), and Eq. (8) for ω 2 l (k) . However, while g(r) can be obtained directly from BD simulations, or may be available from first principles, there is no simple way to determine τ l from first principles. One possibility to proceed is to ignore its k dependence and treat τ l as a free parameter, possibly dependent on Γ, that can be determined from, e.g., fitting the peak positions of Eq. (16) to the peak positions of the spectral densities obtained from experiments or computer simulations. It is expected that such a procedure would yield a dispersion relation lying somewhere between ω = ω 0 (k) and ω = ω l ∞ (k), defined via Eqs. (17) and (21), respectively.
Implementation of the exponential model for transverse viscosity memory function K t (k, t) is a straightforward repetition of the procedure outlined above, and it produces a result for the spectral density that is entirely equivalent to that derived by Murillo for transverse modes in 3D strongly coupled Yukawa liquids [14]. In brief, one assumes K t (k, t) = K t (k, 0)e −t/τt with the initial value given in Eq. (14b), where the second moment of frequency for the transverse mode may also be obtained from the RDF by using Eq. (9). One further obtains the real and imaginary parts of the Laplace transform L[K t (k, t)] s=iω by using expressions analogous to Eqs. (19a) and (19b), which, upon substitution into Eq. (18) yield a spectral density for the transverse mode in the exponential model.
It should be reiterated that, in the case of vanishing relaxation time, τ t = 0, the transverse mode is purely diffusive and the resulting spectral density in Eq. (18) exhibits no dispersion [26]. In the opposite limit of τ t → ∞, viscous damping vanishes, while the dispersion relation, given by can be shown to exhibit a quasi-acoustic behavior in the long wavelength limit. On the other hand, when using the peak positions of the spectral density in Eq. (18) for the exponential model with finite τ t , one encounters a cutoff wavenumber, k c , in the dispersion relation for the transverse mode, which can be estimated from the condition ω 2 t (k c ) = (v th /τ t ) 2 in the limit of vanishing neutral drag [14,22]. However, like in the case of the longitudinal mode, since there is no simple way to determine τ t from first principles, one can again treat it as a free parameter that may be determined from an appropriate fitting, e.g., by using the cutoff values k c observed in the experiments on shear waves in 2D Yukawa liquids [8].

Gaussian model
Instead of using τ l and τ t as free parameters, one may attempt enforcing higher-order frequency sum rules to see if a refinement of the above exponential model can be achieved. However, since the exponential model for the viscosity memory functions does not support moments higher than the second order when used in Eqs. (16) or (18), one needs a model with different time dependence. Besides having to satisfy the second-frequency sum rule via Eq. (14a), it is shown in the Appendix that such model must additionally satisfy both the third-order sum rule giving ∂ ∂t φ l (k, t) t=0 = 0, and the fourthorder sum rule giving We see now that a Gaussian model for the longitudinal viscosity memory function of the form φ l (k, t) = φ l (k, 0) exp [−t 2 /σ 2 l (k)] will satisfy all moments up to and including the fourth if the Gaussian relaxation time σ l (k) is given by with the initial value of the memory function, φ l (k, 0), given by Eq. (14a). Finally, by using the Laplace transform of the Gaussian longitudinal viscosity memory function with s = iω in Eq. (16), we obtain a spectral density which is fully determined by three functions: S(k), ω 2 l (k) , and ω 4 l (k) via Eqs. (17) and (14a) without the need for free parameters. Unfortunately, analytical expressions that can be used for calculation of the fourth frequency moment are rather cumbersome and require a three-particle distribution function, which is difficult to obtain from simulations [26]. Therefore, since calculation of ω 4 l (k) from first principles is impractical, we may use the Gaussian model by computing the fourth moment numerically from Eq. (7a) with n = 2 where P l (k, ω) is obtained from simulation. For consistency, it is then desirable to also compute the second moment from Eq. (7a) with n = 1 using the same P l (k, ω) based on simulation data. In this way, we may use the Gaussian model as a simulation-based, parameter-free test of the quality of the approximation achieved in using the exponential model with a suitable choice of finite relaxation time τ l .
We finally note that a completely analogous development of the Gaussian model can be applied to the transverse mode.

Results and Discussion
We tested several values for the reduced neutral drag coefficient γ in our BD simulation and found no noticeable dependence of the resulting spectra on γ. Moreover, all the quantities used in the GH modeling of the collective modes (e.g., the RDF and the relaxation times in the exponential model of the memory functions for the longitudinal and transverse modes) were also found to be robustly independent of the (small) values of γ used in simulation. Therefore, all results will be shown for a standard value of γ = 0.06, along with the standard screening parameter κ = 1.

Longitudinal Wave Mode
We use coupling strength with values Γ = 20, 60, 100, 200, 600, and 1000 to investigate the longitudinal mode in a broad range of dusty plasma conditions, going from liquid to crystalline states.
As described above, both the static structure factor S(k) and the second frequency moments can be calculated from the equilibrium RDF g(r), which is shown in Fig. 1 for several values of Γ. However, there are significant difficulties related to using S(k) for modeling the dispersion relation of the longitudinal mode at long wavelengths. Namely, when S(k) is calculated from the RDF, its small-k values exhibit a very strong dependence on the fluctuations appearing in the simulation due to the finite number of dust particles [34]. This is expected to give rise to rather noisy dispersion curves, especially in the situations characterized by short relaxation times when the dispersion is dominated by the frequency ω 0 (k) defined in Eq. (17). To remedy the situation to some extent, we resort to calculating S(k) directly from the simulation data by using the long time average of the Fourier transformed particle density in equilibrium, as follows where In Fig. 2, we show the thus obtained results for S(k) for several Γ values, and note that the peak structures in that figure were found to be identical to those arising in S(k) computed from the corresponding RDFs in Fig. 1. In Fig. 3 we show the simulation data for the spectral density of the longitudinal mode, along with the dispersion curves from the exponential model, represented by three solid curves corresponding to, in descending order, infinitely long, finite, and vanishing relaxation times τ l . Also shown in Fig. 3 by a thick black solid line is the dispersion curve obtained from the Gaussian model with the second and fourth frequency moments evaluated from Eq. (7a) with n = 1 and 2, respectively, where P l (k, ω) was taken directly from the simulation data for spectral density. While the three curves for the exponential model extend over the full range of wave-numbers shown in Fig. 3, the results for the Gaussian model are shown over a reduced range of wave-numbers covering the first Brillouin zone only, ka < 6, because computations of the frequency moments via Eq. (7a) were hampered by a significant noise in the simulation spectra at larger wave-numbers and lower coupling strengths. One notices in Fig. 3 that, for the exponential model with vanishing relaxation time (or, equivalently, the delta-function model), the resulting dispersion relation, ω = ω 0 (k) with ω 0 (k) given in Eq. (17), displays some noise coming from the simulation data via Eq. (25). We found this noise to be significantly weaker than the noise arising when S(k) is calculated from the RDF. On the other hand, as the longitudinal relaxation time τ l increases in the exponential model, there is a relative increase in the contribution of the second frequency moment ω 2 l (k) , which, when evaluated from the RDF via Eq. (8), exhibits a much smoother k dependence than ω 0 (k) with S(k) evaluated from Eq. (25).
As mentioned before, the dispersion relation for finite values of the relaxation time τ l in the exponential model for the longitudinal viscosity memory function is obtained from the peak positions in the spectral density given by Eq. (16) with Eqs. (19a), (19b), (17), and (14a). We have chosen the values of τ l which provide the best fit to the peaks in the simulation spectral densities, and we have found τ l to be a relatively weak function of Γ that may be reasonably well approximated by We emphasize that the quality of our choice of the values for τ l in the exponential model is confirmed through a close agreement with the dispersion curves from the Gaussian model for wave-numbers in the first Brillouin zone, as displayed in Fig. 3. It is worth mentioning that the extreme cases of the exponential model, corresponding to the limits τ l → 0 and τ l → ∞, yield two parameter-free dispersion relations for the longitudinal mode, ω = ω 0 (k) and ω = ω l ∞ (k), respectively. Remarkably, it is noticed in Fig. 3 that these two dispersion relations provide good account of the simulation data by straddling the thermal noise, seen in the recorded spectra at the low-to-medium Γ values. Moreover, at the medium-to-high Γ values, one notices that ω = ω 0 (k) closely follows the dispersion relation from the Gaussian model for the wave-numbers in the first Brillouin zone and, in particular, reproduces the near-vanishing of frequency at around ka ≈ 4 for Γ = 600 and 1000. This is somewhat surprising, given that the delta-function model is stretched well into the condensed state at those two coupling strengths.
We further compare in Fig. 4 the dispersion relations ω = ω 0 (k) and ω = ω l ∞ (k) (shown by the upper and lower solid curves, respectively) with the results from the QLCA model [24], shown by the lower dashed curve. It is immediately obvious that the QLCA dispersion does not reproduce the direct thermal effect, which is responsible for a quasi-linear increase in the peak frequencies of the simulation spectra at higher wave-numbers and lower coupling strengths. In a previous work, we have shown that this deficiency can be easily rectified by a simple extension of the QLCA model, labeled as the EQLCA model in Ref. [22], which is shown by the upper dashed curve in Fig. 4. One notices a surprisingly good agreement between the EQLCA dispersion relation and the curve ω = ω l ∞ (k) from the exponential model for τ l → ∞. While an agreement with the QLCA dispersion relation was expected at the higher Γ values and lower k values owing to the fact that the QLCA model inherently assumes ωτ l ≫ 1 [35] and the direct thermal effect is relatively weak at low k and high Γ values [22], it is remarkable how the parameter-free dispersion relation ω = ω l ∞ (k) provides justification for the empirically derived EQLCA model over the broad ranges of wavenumbers and coupling strengths.
We next compare the performances of the exponential model with finite relaxation time τ l from Eq. (27) and the Gaussian model with the frequency moments evaluated from the simulation spectra via Eq. (7a) by using these two models in Eq. (16) to predict the frequency dependent spectral profiles at several fixed wave-numbers. Given that the half-width at halfmaximum (HWHM) of these profiles may be directly related to the wave-number dependent damping rate of the longitudinal DAW mode, in this way we demonstrate the advantage of using the memory function formalism within the GH model in tackling the difficult issue of damping of the collective excitation modes in dusty plasmas.
Specifically, we show in Figs. 5-10 the frequency dependencies obtained from the spectral density in Eq. (16) with the exponential and the Gaussian models for four wavenumbers, ka = 0.61, 1.17, 2.29, and 3.36, and compare them with the corresponding profiles of the simulation spectral density for Γ = 20, 60, 100, 200, 600, and 1000, respectively. One notices in Figs. 5-10 a fair agreement between the simulation data and both GH models, which is especially good for lower Γ values where the hydrodynamic regime is expected to be more pronounced, while for very high coupling strengths, say, Γ > 200, the agreement begins to deteriorate. Comparison of both GH models with the simulation data may be considered fairly reasonable, even at high coupling strengths where the hydrodynamic model has been stretched well into the crystalline state.
Finally, in Fig. 11, we show a comparison between the wave-number dependent HWHM values, obtained from the simulation spectral density profiles and from Eq. (16) where we used both the exponential model with finite τ l values given in Eq. (27), and the Gaussian model with the frequency moments evaluated from the simulation spectra via Eq. (7a). One notices that all three sets of data are noisy due to the noise in the simulation spectra (one recalls that, in the case of the exponential model, the noise stems from ω 0 (k) with S(k) evaluated from Eq. (25)). As expected, the HWHM values from the Gaussian model are systematically somewhat larger than those from the exponential model, but they are both seen to be in reasonably good agreement with the simulation data, even though the agreement becomes blurred by the increased noise at the highest Γ values, say Γ > 200. One also notices in Fig. 11 that all damping rates roughly follow a quadratic dependence on k, with increasing opening of the parabola as Γ gama increases, and with the limiting value as k → 0 being close to the damping rate due to the neutral drag, given by γ/2 = 0.03.

Transverse Wave Mode
Transverse, or shear waves are expected to only exist for higher coupling strengths Γ [14,8]. In the following we compare the simulation results with the GH approach for a typical case of Γ = 100. Using Eq. (9) to evaluate the second frequency moment from the RDF, we have obtained a dispersion relation from the peak positions of the spectral density in Eq. (18) for the exponential model with τ t = 5.0 (note that typical values of τ t for the transverse mode are generally found to be larger than those for the longitudinal mode, in agreement with the remarks of Boon and Yip [26]). The results are shown in Fig. 12, along with the dispersion relation in the limit of infinite relaxation time, τ t → ∞, given by ω = ω t ∞ (k) with Eq. (22), and are compared with the simulation spectrum for Γ = 100. Both theoretical dispersion curves show good agreement with the simulation data for ka < 4 and, while the simulation spectrum is not conclusive about the cutoff wave-number, the model with finite τ t clearly exhibits a cutoff at k c ≈ 0.5/a [14,8].
Further, several profile plots of the spectral density are shown in Fig. 13 for Γ = 100, where the results from Eq. (18) for both the exponential model with τ t = 5.0 and the Gaussian model with the frequency moments evaluated from the simulation data via Eq. (7a) are compared with the simulation profiles for ka = 0.61, 1.17, 2.29, and 3.36. One notices that the agreement of the exponential model with the simulation spectral profiles is quite satisfactory for the lower two wave-numbers, but the simulation profiles at the two higher wave-numbers are seen to be much broader than the corresponding peaks from the exponential model. On the other hand, the Gaussian model does not reproduce the dispersion of the transverse mode at all, since the corresponding profile curves peak at zero frequency for all wave-numbers, as displayed in Fig. 13. Nevertheless, the broadening of the simulation profiles at the two higher wave-numbers in Fig. 13 seems to be better echoed by the Gaussian than by the exponential spectral profile curves.

Concluding Remarks
We have carried out Brownian Dynamics simulation of a 2D layer of strongly-coupled dusty plasma and extracted the equilibrium radial distribution function, static structure factor, and the spectral densities for both the longitudinal and transverse dust acoustic wave modes.
We have then shown that the memory function formalism of the theory of generalized hydrodynamics provides good semi-analytic results to model these wave modes. In particular, following the approach of Boon and Yip [26], we have developed an exponential model for the viscosity memory functions that satisfies the second frequency sum rule with the relaxation time being used as a fitting parameter. With such exponential model we have obtained good fits to the simulation data for the longitudinal dispersion curves over a wide range of coupling strengths, 20 ≤ Γ ≤ 1000, as well as reasonable match with the spectral density profiles for several values of the wave-number k. Next, we have extended the theory in a straightforward manner to satisfy the third and fourth frequency sum rules, providing a parameter-free Gaussian model for the viscosity memory function. The results using this model provided a successful test for our choice of the relaxation time in the exponential model in modeling both the dispersion relations and the spectral density profiles. Moreover, we pointed out that the limits of an infinitely long and infinitesimally short relaxation times in the exponential model represent two parameter-free models that are fully determined by the equilibrium radial distribution function and the static structure factor of the system. It was then shown that these two limits provide good upper and lower bounds for the thermal noise observed in the simulation spectra for the longitudinal mode. Comparisons of the dispersion relations obtained from these two limits of the exponential model with the dispersion relations obtained from both the Quasi-localized charge approximation (QLCA) and the extended QLCA (EQLCA) showed that the limit of an infinitely long relaxation time reproduces remarkably well all the features of the EQLCA model, including the direct thermal effect. On the other hand, the limit of an infinitesimally short relaxation time was found to reproduce the near-vanishing of the dispersion relation, which is seen in the simulation spectra at certain short wave-lengths for large coupling strengths, but is not reproduced by the QLCA group of results. In addition, we have also tested the memory function formalism within the theory of generalized hydrodynamics against the simulation results obtained for transverse wave modes at the coupling strength of Γ = 100, and found that the exponential model with relaxation time  Figure 11. Half-width at half maximum of the spectral profile curves versus the reduced wavenumber ka for κ = 1.0, γ = 0.06, and Γ = 20, 60, 100, 200, 600, and 1000, with the (blue) dots representing the simulation data, the (red) noisy solid curve representing the exponential model with finite relaxation time, and the squares representing the Gaussian model. used as a fitting parameter can provide a reasonable estimate for the wave dispersion relation, including the appearance of a cutoff wave-number. Both the peak positions and the widths in the spectral density profiles from the simulation were well reproduced by the exponential model at long wavelengths, whereas the broadening of these spectra at short wavelengths was better echoed by the Gaussian model.
The most important finding of the present analysis is that the wave-number dependent damping rates, which were extracted from the simulation spectra, were found in good agreement with the results from both the exponential and Gaussian models, testifying to the strengths and potentials of using the memory-function formalism within the GH theory in describing the damping of the longitudinal dust acoustic wave in strongly-coupled dusty plasmas.
Finally, we have included the effect of neutral drag due to the collisions of dust particle with neutral molecules in the background plasma via a (local in time) drag coefficient. While we did not find that the neutral drag affects our modeling of dispersion relations for both the longitudinal and transverse waves in any significant way, the overall damping rates for the longitudinal wave showed possibly important roles of the neutral drag at long wavelengths and high coupling strengths where viscous damping is expected to be reduced. This aspect of the GH theory shall be studied in future work.
In conclusion, the memory function formalism of the generalized hydrodynamics provides a very promising route to analytical modeling of spectral density of collective modes in strongly coupled 2D dusty plasmas and, in particular it demonstrates a good account of viscoelastic effects which are not well described by other analytical methods.