Brownian motion of a plasma crystal

The dynamics of a plasma crystal under the action of random external forces is considered. The pair correlation functions of the particle displacements are calculated in the harmonic approximation by using the Langevin equations. The case of a planar hexagonal lattice is analyzed in more detail. Analogues of the Van Hove singularities in the spectral densities are discovered.


INTRODUCTION
Dusty plasma is an ensemble of microscopic particles (dust grains) immerged in the ambient ionized gas consisting of weakly interacting electrons, ions, and neutral atoms. Over the past three decades, the physics of dusty plasmas has become a separate research area. Reviews on the main results obtained in this field can be found, e.g., in [1][2][3].
Dust grains are centers of recombination of the ambient plasma; therefore, they usually acquire a large negative charge. At a sufficiently high dust number density, the ensemble of dust grains turns into the liquid or crystalline state. In addition to regular forces exerted on the dust grains by the ambient gas, there are also random fluctuation forces, which result in the Brownian motion of dust grains.
Observation of the Brownian motion of dust grains is an efficient method of dusty plasma diagnostics (see, e.g., [4,5]). The behavior of dusty plasma is often simulated by numerically solving the equations of motion of dust grains under the action of random forces (see, e.g., [6]). This paper is dealing with the effect of interaction between dust grains on the spectrum of their Brownian motion. An ensemble of interacting grains with coordinates of is considered, where the index labels individual grains subject to external random forces. The origin of random forces in actual dusty plasmas is not yet completely understood. Therefore, for simplicity, we will consider the case of white noise, when the random forces acting on separate dust grains are uncorrelated. The interaction between dust grains is assumed to be fairly strong. We also assume that, in the absence of external forces, the grains are at the points of a perfect lattice. When such a crystal is affected by stochastic forces, the grains are randomly displaced from their equilibrium position. The problem is to calculate the correlation functions and their spectral densities.
We note that the term "Brownian motion" is often understood as the motion of an individual particle in the field of random forces resulting in a linear (or more complicated) time dependence of the mean square particle displacement. In this paper, we consider the term "Brownian motion" in a wider sense. A wellknown example of such motion is the dynamics of a harmonic oscillator under the action of random forces, in which case the mean square displacement is limited. In fact, in this work, we consider the Brownian motion of an infinite ensemble of oscillators.
Besides dusty plasmas, similar problems also arise in the physics of colloidal crystals (see, e.g., [7]). However, the motion of colloidal particles in the ambient medium is characterized by high friction forces, due to which the particle inertia can be neglected. For dusty plasmas, the opposite limit of weak friction, when the particle inertia plays a significant role, is of greater interest.
The paper is organized as follows. In Section 2, the Langevin equations for an infinite number of particles are written out and a general expression for the spectral density of displacement correlators in terms of phonon spectra is obtained in the harmonic approximation. Then, the application of the general theory to the case of a two-dimensional hexagonal lattice is analyzed in detail. Oscillations of a hexagonal plasma crystal have been studied sufficiently well (see, e.g., [1,3]); therefore, in Section 3, we only briefly list their properties required for our further analysis. In Sections 4 and 5, we discussed typical features of the spectra of the correlation functions.

DUSTY PLASMA
When describing random processes, we follow the terminology and notation of [8]. We also use some concepts of the solid-state theory that are unusual for plasma physics. As soon as is possible, we explain them. A more detailed discussion can be found, e.g., in [9].

BASIC EQUATIONS
We consider an ensemble of particles with coordinates and momenta that interact via the force potential . The particles move in a neutral gas, and each of them is subject to the friction force and the random force . The dynamics of the particle ensemble can be described by the set of Langevin equations, (1) where is the particle label and is the coefficient of friction in the ambient gas. The potential energy of the particle ensemble described by Eqs. (1) is (2) The generalization of Eqs. (1) with allowance for external forces, the nonreciprocity and nonpotentiality of interaction between particles, and other complicating and refining factors is quite obvious.
In what follows, we assume that the external forces are caused by a stationary random process with the correlation function (3) where and is intensity of the source of external forces. Here, the angle brackets denote averaging over the ensemble and the indices correspond to the Cartesian components. The forces affecting different particles are assumed to be uncorrelated.
If the random forces are caused the effect of the ambient plasma, then, generally speaking, the random process described by Eq. (3) is not -correlated. Moreover, there are no reasons to consider that the fluctuations of the forces acting on different particles are independent. The calculation of the spectral density of plasma fluctuations is a separate problem. As will be shown below, the explicit form of the correlation function of the external noise is of minor importance; therefore, for simplicity, we will use expression (3).
Let us assume that the interaction forces between particles are so strong that, in the absence of external fields ( ), the particles are at some equilibrium positions with zero velocities. We designate the equi- t f librium coordinates as and small displacements from the equilibrium positions as . In other words, the particles are numerated using the vector index , which takes discrete values. Making a change of variables in Eqs. (1) and (2) and expanding them in powers of , we obtain the following linearized set of Langevin equations: (4) where the components of the force matrix are (5) and . Hereinafter, the dot indicates the scalar product of two vectors or the product of a matrix and a vector.
Our aim is to calculate the correlation functions of the particle displacements and their spectral densities . Turning to the Fourier transforms in linearized equations (4), we find that the function satisfies the following equation: (6) where the matrix is (7) and is the unity matrix. Note that matrices (5) and (7) are symmetric. We multiply Eq. (6) by the complex conjugate equation and make the change . Then, we average it over the ensemble and take into account that, for a stationary random process, the following relation is true: . As a result, we arrive at the equation for the spectral density, (8) We suppose that the particles in their equilibrium positions form a crystalline structure with a Bravais lattice, i.e., there is one particle per cell. Then   For the three-dimensional lattice, we have , whereas for the twodimensional lattice considered below, we have . Due to the homogeneity of the problem, the correlation functions related to different particles should depend only on the interparticle distance, i.e., we can assume that the relation is true. In this case, it is convenient to solve Eq. (8) by means of the discrete Fourier transformation. For any function defined on a periodic lattice, we can write (9) The function is periodic on the reciprocal lattice, i.e., , where are arbitrary integers and are the primitive vectors of the reciprocal lattice, defined as . The function is completely determined by its values in the primitive cell of the reciprocal lattice. It is convenient to choose the first Brillouin zone as a primitive cell of the reciprocal lattice. The transformation inverse to transformation (9) can be written as (10) where integration is performed over the Brillouin zone and is the volume or area of this zone.
After performing transformation (9), Eq. (8) for the correlation function takes the form . The solution to this equation is , where (11) It follows from definition (7) that the Fourier transform of the matrix is , where the components of the force matrix with allowance for expressions (5) are (12) The physical meaning of these formulas is quite simple and clear. In the absence of damping and external action, the eigenvalues of force matrix (12) determine the spectrum of eigenmodes of the crystal. The matrix in expression (11) describes the collective response of the crystal to a random external action in the presence of damping. In fact, this matrix acts as a filter selecting certain frequencies from the spectrum of external noise. It is obvious that, for a sufficiently small damping ( ), the response of the crystal is resonant. In a certain frequency range, the components of the correlation matrix obtained from expression (11) by means of transformation (10) can be = , , The frequency dependence of this matrix possesses some interesting features, which will be discussed below by using a particular example.

PLANAR HEXAGONAL LATTICE
Let the particles in their equilibrium positions form a planar hexagonal lattice and the distance between the neighboring particles be . Bearing in mind applications of this theory to plasma dust crystals, we assume that the particle charge is and, at short distances, the particles interact via the Coulomb repulsion potential, i.e., at . We assume that the potential is isotropic and depends only on the distance. It is convenient to introduce dimensionless variables in which the distance is normalized to the lattice constant and the frequency is normalized to . We also assume that all particles move only in the plane. In the harmonic approximation, the displacements of particles along the axis do not depend on the particle motion in the plane; therefore, they can be analyzed separately.
We choose the following primitive vectors of the hexagonal lattice in Cartesian coordinates (see Fig. 1a): and . The primitive vectors of the reciprocal lattice are and (Fig. 1b). The Brillouin zone is a regular hexagon with a side of and an area of . It is shown by the dashed line in Fig. 1b. One vertex of the Brillouin zone is at the point , while the others are obtained by successive rotations by an angle of .

Lattice Sums
As was mentioned above, the lattice sums of form (12) determine periodic functions of the wave vector . The hexagonal lattice is also invariant with No. 6 2017 IGNATOV respect to rotations, forming the group . For an isotropic interaction potential , it is expedient to take into account additional restrictions on force matrix (12) that arise due to this symmetry. It is easy to prove that, when the lattice rotates by an angle of , force matrix (12) is transformed as (13) where is an orthogonal matrix.
Any symmetric matrix having symmetry property (13) can be expressed in terms of one scalar function (14) where is the exterior product of two vectors, i.e., it is a matrix with the components . Thus, the function can be expressed in terms of the initial matrix, The validity of relations (14) and (15) can be verified by direct substitution into expression (13). With the help of relations (15) and (12), the function can be written as the lattice sum, where and the coefficients and depend only on the distance from the lattice center, .
We will use some general properties of sum (16). Obviously, this even function ( ) is periodic on the reciprocal lattice, i.e., , where are arbitrary integers. In a general case, it has no symmetry with respect to rotations, except for certain values of the wave vector. First, these are the vertices of the Brillouin zone (Fig. 1b). As follows from Eq. (16), the following equality is true at these points: Second, the gradient of function (16) vanishes ( ) at the midpoints of the edges of the Brillouin zone, i.e., at the points , , and so on (Fig. 1b).
When numerically calculating force matrix (12) or sum (16), it is necessary to take into account the symmetry of the problem. We divide all points of the hexagonal lattice into separate shells, which are concentric circles of radii . The number of the lattice points in each shell are multiples of , and this set of points remains the same when the lattice rotates by an angle of . The symmetry of the problem remains unchanged if we rewrite function (16) in the form , where each term of the series corresponds to the sum calculated over a certain shell. For an interaction potential sufficiently rapidly decreasing at large distances, we can limit ourselves to a small number of shells. The first two terms of such a shell expansion of Eq. (16) are (18) where .
All the plots presented below are depicted using the screened Coulomb potential in the approximation of the nearest neighbor interaction.

Phonon Spectra
We designate the eigenvalues of the force matrix as . These quantities determine the spectrum of natural vibrations of the crystal and can be expressed in terms of the function by using representation (14), Since the equality is true, the proper frequencies are obviously invariant with respect to rotations ( ).
The above-discussed properties of sum (16) also manifest themselves in the eigenfrequencies. In view of relations (17), the spectrum becomes degenerate at the vertices of the Brillouin zone, i.e., , and matrix (12) is proportional to the unity matrix. Using sum (16), we can show that, in the vicinity of the vertices of the Brillouin zone at , the frequencies depend on the wave vector as , where is a certain positive factor. By analogy with graphene (see, e.g., = , , , , 1 3 2 7 3... l 6  [10]) and other similar substances, we will call these singularities the Dirac points. 1 The group velocities of both vibration branches vanish at points on the edge of the Brillouin zone (Fig. 1b). The positions of these singular points are determined only by the lattice symmetry. In addition, there may be other critical points the positions of which are determined by the type of the interaction potential.
In the long-wavelength acoustic limit, a crystal with the symmetry group has the same properties as a solid body that is isotropic in the plane (see, e.g., [11], Section 10). Therefore, at , matrix (12) depends only on the wave vector, i.e., the relation is true, where are the velocities of transverse and longitudinal acoustic vibrations. In this limit, the frequencies given by expression (19) are equal to and function (16) can be written as (20) For an arbitrary interaction potential, the velocities of acoustic waves can be expressed as sums (21) Near the edges of the Brillouin zone, the oscillations cease to be purely longitudinal or transverse. In particular, at the Dirac points, there is no difference between longitudinal and transverse waves. Neverthe-1 It would be unreasonable to draw direct analogy between graphene and plasma crystal; however, in both cases, singular points of this type appear due to the symmetry of the system. less, for simplicity, we will refer to the low-frequency and higher-frequency modes as transverse and longitudinal modes, respectively. Figure 2 shows the contour lines of the frequency in the Brillouin zone for (a) transverse and (b) longitudinal oscillations, calculated for the screened Coulomb potential in the approximation of the nearest neighbor interaction. Figure 2 also shows the directions of oscillation polarization, i.e., the eigenvectors of force matrix (12) corresponding to the frequencies , as well as the critical points at which the group velocity vanishes. The Dirac points at the vertices of the Brillouin zone are clearly seen in Fig. 3, which shows three-dimensional plots of both oscillation branches.
It is seen from Figs. 2b and 3 that, in addition to the critical points at the boundary of the Brillouin zone, there are also six saddle points in the spectrum of longitudinal oscillations, one of which has the coordi-   Fig. 2b, these points are designated by crosses. For the screened Coulomb potential in the approximation of the nearest neighbor interaction, we obtain .

CORRELATION FUNCTIONS
Correlation matrix (11) can be written as , where is force matrix (12). It is convenient to write it in an explicit form as (22) where is the unity matrix, are eigenfrequencies (19), and The validity of relation (22) can be ascertained by going over to a basis in which the symmetric matrix is diagonal.
In the case of small damping ( ), function (23) has a sharp maximum at the frequency and, in the limiting case, transforms into the Dirac function (24) Thus, in the case of small damping, expression (22) as a function of the wave vector is localized in the vicinity of the contour lines of the frequency shown by the solid lines in Fig. 2.
We are interested in the correlation matrices of particle displacements in the crystal, which can be calculated by applying Fourier transformation (10) to matrix (22), In the general case, integrals over the Brillouin zone in the form of expression (25) can be calculated only numerically; however, some limiting cases can be studied analytically.
Let us consider the limit of weak damping and sufficiently low frequencies, . In this case, we can make substitution (24) and use expression (20). After these simplifications, integrals of form (25) can easily be = , and . It is clear from this expression that, in the case of weak damping, correlations between the displacements of different particles are caused by the exchange of acoustic waves with the frequency and wave vectors . Note that the matrix describing the mean-square displacements of an individual particle is always proportional to the unity matrix, i.e., the displacements are independent of the directions of the crystal axes.
At each point of lattice , the correlation matrix is real and symmetric. We designate the eigenvalues of the matrix and the normalized eigenvectors as ( ) and , respectively. Depending on the frequency, the eigenvalues can be either positive or negative. If a certain eigenvalue is positive, e.g., , then, on the average, random displacements of the central particle with and the particle at the lattice point in the direction of the vector are inphase. However, if , then, on the average, the particles shift in opposite directions. Figure 4 gives an idea of the structure of the correlation matrices calculated by expression (25). The figure is constructed using expression (26) for and the screened Coulomb potential with . Matrices with the same signs of eigenvalues are shown by ellipses with semiaxes proportional to and directed along the eigenvectors . Red and blue ellipses correspond to matrices with positive and negative eigenvalues, respectively. Matrices with different signs of eigenvalues are represented by crossed segments having lengths proportional to and directed along the eigenvectors. Red and blue lines correspond to positive and negative eigenvalues, respectively. Figure 4 shows that, for the nearest neighbors ( ), the correlation matrices are nearly isotropic. On the second shell of the hexagonal lattice ( ), the correlation matrices become more anisotropic and the amplitude of the radial displacement becomes substantially larger than that of the transverse dis- l a a R l a a l a a − 1 2 ( ) a a = , placement. As the distance increases ( ), the radial displacements of particles become, on the average, directed oppositely, while the transverse displacements occur in the same direction. It is clear from expression (26) that, at large distances, the matrix elements oscillate with increasing distance , whereas the displacement amplitudes decrease as .

VAN HOVE SINGULARITIES
As was mentioned above, the dependence of the oscillation eigenfrequencies on the wave vector has singular points at which the group velocity vanishes or the two branches of the spectrum merge. This means that, for a sufficiently weak damping, matrices (25) as functions of the frequency have singularities related to the topological rearrangement of the contour lines . Similar singularities in the spectral density of excitations appear in the solid-state theory, where they are called Van Hove singularities.
For the hexagonal crystal under consideration, there are three types of such singularities. First, these are points (Fig. 1b) at the edge of the Brillouin zone, where the frequency of longitudinal oscillations is maximal. In the proper coordinate system in wave vector space, the frequency near these points can be written as , where and are positive coefficients. In the vicinity of the maximum frequency , integrals of form (25), (where is a certain smooth function), depend on the frequency as , i.e., they vanish abruptly. Near the saddle points, the relation is true. At the fre- . There are two sets of saddle points. For transverse waves, they are located at points (Fig. 1b) and the corresponding frequency is equal to . For longitudinal waves, the saddle points are located inside the Brillouin zone (crosses in Fig. 2b). We will designate their frequency as . Finally, integration over the vicinities of the Dirac points, where the relation is true, makes the contribution of the form . Hence, in this case, there is no singularity; however, we can expect the appearance of a local minimum of the elements of matrix (25).
The matrix determined by expression (25) is always proportional to the unity matrix, i.e., both its eigenvalues are equal. Figure 5 shows the frequency dependence of the eigenvalue of the matrix calculated numerically for the screened Coulomb potential under conditions of finite damping. The above singularities turn out to be smoothed out, but they are still clearly visible even at a sufficiently strong damping. The most pronounced is the Van Hove singularity near the saddle point of transverse oscillations with . The eigenvalues of the correlation matrix for (i.e., for different particles) have similar properties. However, in this case, the eigenvalues are different and their signs can depend on the frequency . The structure of the correlation matrices changes substantially in the vicinity of the frequencies corresponding to the saddle points of the phonon spectrum. Figure 6 shows the matrices obtained by numerical integration of expression (25) for several shells at the frequency . It can be seen from the figure that, for neighboring particles ( ), the polarization of displacements does not change qualitatively; however, at a larger distance ( ), the eigenvalues of some matrices change their signs. Similar properties are also observed near the saddle point of longitudinal oscillations at the frequency (Fig. 7).
As the frequency varies in the vicinity of the Dirac point, the structure of the correlation matrices changes insignificantly. Finally, considerable correlations between the displacements of distant particles are observed at frequencies slightly below the maximum frequency; however, these correlations disap-= 1 l = 2 l ω ≈ Ω sl pear almost completely as the frequency slightly increases (Fig. 8).

SIMULTANEOUS MOMENTS
Using the matrices determined by expression (22), we can calculate the other spectral densities. In particular, the spectral density of the velocity correlator is . Applying the frequency Fourier transformation to expression (25), we can calculate the time dependences of the velocity correlators.
Since, for the Fourier components, the relation is true, we have On the other hand, we have and, in the long-wavelength limit , the integral of expression (22) depends on the wave vector as . Hence, spatial Fourier transform (25) diverges logarithmically. This divergence illustrates the general assertion on the absence of strict long-range order in two-dimensional systems (see [12], Section 137).

CONCLUSIONS
In this work, we have studied the correlation functions of particle displacements by using a planar hexagonal crystal as an example. In particular, we have shown that certain features of the phonon spectra manifest themselves in the frequency dependences of the spectral densities.