Spontaneous disordering of a two-dimensional (2D) plasma crystal

Spontaneous disordering plays an important role in the physics of highly ordered complex plasmas. In this paper, an analytical theory is proposed for the process of ‘cold amorphization’, which has been observed. This consists of splitting a plasma crystal into sub-domains, followed by disordering. The results of recent simulations and experiments showing such spontaneous disordering have been reviewed and interpreted in this paper. Complex plasmas can serve as a powerful tool providing fundamental insight into this process generically.


Introduction
In this paper, we address 'amorphization' (through the appearance of a network of sub-domains) of a cold plasma crystal associated with spontaneous disordering of a highly ordered twodimensional (2D) complex plasma.
Low-pressure, low-temperature plasmas are called complex plasmas if they contain microparticles as an additional thermodynamically active component. In the size domain of 1-10 µm (normally used in experiments with complex plasmas), these particles can be visualized individually, hence providing an atomistic (kinetic) level of investigations [1,2]. The interparticle spacing is typically of the order of 0.1-1 mm and the characteristic time scales are of the order of 0.01-0.1 s. These unique characteristics allow investigations of many dynamical processes at the kinetic level.
One of the interesting aspects of complex plasma research is the possibility to levitate a monolayer of particles under gravity conditions. In this case, the particle suspension has a flat practically 2D structure. Depending on the discharge conditions, the monolayer can have crystalline or liquid order. 2D configurations of dust particles either in crystalline or liquid state were successfully used to study phase transitions, dynamics of waves and many transport phenomena in complex plasmas (see e.g. [2]- [4]). A symmetry disordering accompanying the crystalline-liquid phase transition has been investigated experimentally in [5]- [7].
The microparticles, collecting electrons and ions from the plasma background, become charged (most often negatively [1]) and hence can be confined by external electric fields. The configuration of the confining forces affects the geometry and actual structure of the microparticle cloud. In rf discharge complex plasmas, the particles are self-trapped inside the plasma because of a favorable configuration of the electric fields [8].
3 2D configurations constitute an attractive simplification, significantly lowering the description difficulties. It is not surprising then that there are many simulations of 2D plasma crystals in the literature; it is impossible to cite here all related studies. A comprehensive study of 2D classical clusters was started by Bedanov and Peeters [9] and then continued by many other authors. The first step in understanding the effects of the confinement and the anisotropy of confinement was taken by Cândidoy et al [10]. Defect excitation, melting of a 2D Coulomb cluster with charged particles interacting with r −1 force, was studied numerically by Lai and Lin [11], perhaps the first discussion of the symmetry breaking by packing a triangular lattice in a circular boundary. (For a review of the dynamical heterogeneity in 2D dusty plasma liquids, see also [12].) The ground states of the particle configurations at different interaction ranges and the number of particles, and the statics and dynamics of 2D mesoscopic clusters have been investigated by Kong et al [13]. The global structure of large finite 2D Yukawa clusters confined by the parabolic potential was simulated and analyzed by Totsuji et al [14,15] and by Zhdanov et al [16]. Controversial remains the problem of the so-called hexatic state in plasma crystal studies [17]- [19]. Experiments and simulations of the hexatic state have recently been performed by Sheridan [18,19]. Spontaneous disordering, splitting of the lattice layer and assembling of a system of sub-domains are of principal importance in this problem. The actual influence of the shape of the confinement well on the plasma crystal dynamical properties has been discussed by Durniak et al [20].
Here we would like to highlight and report on the physics of spontaneous disordering of a plasma crystal. The results of recent experimental observations and simulations cannot be properly addressed without a deep understanding of this important issue. However, to the best of our knowledge, such a detailed analysis is still outstanding. This paper is organized as follows. Following this introduction (section 1), we formulate the main features of the model (section 2) and then start the discussion with the crystal fluctuation spectra and introduce the correlation length, which is directly related to the confinement strength (section 3). The influence of the confinement is discussed in section 4. This gives us an opportunity to link qualitatively spontaneous disordering with 'amorphization' -the splitting of a crystal into sub-domains. A simple criterion of such 'amorphization' is proposed in section 5. In section 6, the results obtained are summarized and the main conclusion is presented.

General aspects
Mathematical models (see e.g. [6], [14]- [16], [19]- [23]) used in analytical theories and molecular dynamics simulations of a 2D complex plasma at the kinetic level are traditionally based on a system of N interacting particles described by the Newtonian equations where M is the particle mass and r i the coordinate of the i-particle. The first term on the right-hand side describes symbolically the mutual interparticle repulsive interactions (via e.g. the Yukawa forces (2)), whereas the second one represents a sum of all possible external 4 interactions, e.g. friction (external dissipation) −mγ˙ r i due to collisions with neutrals (Epstein drag or neutral damping with a damping rate γ ), excitation, 'heating' and confinement forces. Actually, this is a 'standard' model, and modifications are only in the concrete realization of friction, confinement, excitation (for numerical realization of the wave excitations, see e.g. [20] and references therein) or heating applied.
For example, in [6,22], randomization (heating) was achieved by applying a random force (produced by a standard random number generator) to each individual particle; the particle kinetic temperature was dissipated through neutral gas drag, and no thermostat was used to control it. In [19], the kinetic particle temperature of the lattice layer was introduced, implementing the Metropolis algorithm, and there is formally no neutral gas drag involved. The co-existence of these two approaches reflects two different scenarios of complex plasma dynamics: either the particles are strongly coupled with a background, providing Brownian dynamics (diffusion), or the frictional coupling is relatively weak, mutual particle interactions are dominating and a complex plasma can be treated as a dynamical system with background friction [24].
Concerning particle confinement, a parabolic confinement well is often suggested in simulations. In this case, the confinement force −m 2 r i increases linearly with distance [16,19,20]. Higher-order polynomials |f conf | ∝ r n−1 , n > 2 are also a promising tool, allowing a more homogeneous highly ordered crystal structure with a weaker density gradient in the center. At fixed crystal size, such confinement is steeper toward the crystal edges but softer in the center of the configuration compared to the parabolic well, and the defect-free central region is larger [20]. This is in good qualitative agreement with the conclusions drawn in section 5: 'softening' of the confinement might help us to simulate crystals that are practically free of defects.

Confinement-influenced acoustic fluctuation spectra
It is well known that two broken symmetries distinguish the crystalline state from the liquid: the broken translational order and the broken orientational order. In 2D for ordinary crystals, it is also well known that even at low temperatures the translational order is broken by spontaneous disordering mediated by thermal fluctuations [25]. As a result, the fluctuation deflections (disordering) grow with distance and translational correlations decay (algebraically, see [26]).

Particulars of the analytical model
2D plasma crystals also obey this common rule. Nonetheless, there is an aspect that is worth noting: the character of disordering may be strongly affected by the confinement forces. The new approach developed below enables us to introduce a specific correlation length (see equation (23) below) defined by the confinement forces.
Usually such an 'in-plane' confinement is due to the bowl-shaped potential well selfmaintained inside the discharge chamber, which to first order is approximately parabolic (see e.g. [19,27]), that is, U conf ∝ r 2 , where r is the distance. We introduce U conf = 1 2 M 2 r 2 , where M is the particle mass and is the confinement parameter [16]. (The 'out-of-plane' confining forces, controlling the position of the entire lattice, are normally much stronger; below we consider the 'pure' 2D case, assuming an absolutely stiff out-of-plane confinement.) We consider a 2D cluster of equally charged particles interacting pairwise via a Yukawa (the screened Coulomb) force, where R i, j = r i − r j is the relative coordinate and R = |R i, j | is the distance between the particles i, j with the coordinates r i , r j ; q is the particle charge and λ is the screening length.

Phonon contribution to the energy
The key point is to analyze the fluctuation spectra. The fluctuation spectra can be calculated in the following manner. (We follow the direct expansion method developed in [21].) The longrange phonon contribution to the free energy of a 2D system of particles interacting via Yukawa forces and confined by a shallow isotropic parabolic well can be conveniently represented as where c tr,l are the transverse (shear wave) and the longitudinal (compressional wave) sound speed, respectively, V k , D k are the Fourier components of the vorticity V = curl z u and the divergency D = div u of the particle displacements u = k u k exp(ikr), and k is the wave vector (k = |k|). The unperturbed crystal structure is supposed to be hexagonal. For a definition of the phonon speed, see e.g. [21,28,29]. The ratio of the phonon speeds c 2 l /c 2 tr depends only on the interaction range and not on the particle charge; typically, c 2 l /c 2 tr ∼ 10-30 1 at the values of the interaction range a/λ = 0.5-2 most interesting in complex plasma experiments [2,3,21,30,31]. Here, a is the interparticle separation or the lattice constant of the ideal plasma crystal.

A probability-based algorithm of averaging
The relationship (3) provides the probability of the fluctuations (see e.g. [26]) w ∼ exp(− U/T ), allowing the calculation of all necessary statistical means.
As an example of the averaging procedure, let us calculate the averaged spectral intensity of the fluctuations u k . The calculations are straightforward; therefore we skip most of the routine details. Obtaining the mean values (here of a spectral intensity), we introduce the average by the relationship where dS D,V = ρd ρdϕ | D,V ; as a pair (ρ, ϕ) has to be substituted to (5) either (|D k |, arg D k ) or (|V k |, arg V k ) depending on whether the divergency or vorticity fluctuations are considered. Next, noting that 6 we get where v T is the particle thermal velocity. It is worth noting that the averaged spectral intensity does not depend explicitly on the dissipation rate. Hence, the analytical model is equally applicable in situations with strong or weak friction. Note also that, even if friction is strong, it is not an obstacle to interpreting the fluctuations as phonon excitations. For instance, the phonon excitations are important in a highly ordered system of charged particles immersed in a colloid, which is a system with evidently strong friction [28].
Evaluating the mean values of the functions given by the Fourier arrays, one needs to average the exponential functions exp(ikx) . Performing averaging, it is enough to apply the well-known identity (it is assumed that x = 0) and then follow the procedure given e.g. in [26,32]. Calculating the density-density correlations shown in section 4, we adopt exactly this procedure. Averaging in (8) is carried out by applying rule (5).

Averaging with the help of the fluctuation-dissipation theorem (FDT)
Equivalently, the fluctuation-dissipation theorem (FDT) can be applied to implement the averaging. For the Yukawa lattice layer, the space-time fluctuations of the particle positions induced by a given random force f can be written as where ω is the frequency, G 1,2 = (ω(ω + iν) − 2 1,2 (k) − 2 ) −1 with the dissipation rate ν are the propagators corresponding to the eigenvectors e 1,2 = e 1,2 (k) and the eigenvalues 1,2 (k) of the dynamical matrix introduced in [22].
It is easy to verify that the inverse Fourier transform (over ω) yields The cores K 1,2 (τ ) give the response functions. The spectrum of the fluctuations is stationary if the noise in the origin of the random force is stationary. Assuming a white noise delta-correlated in time, we obtain Further, substituting |e 1,2 A| 2 = 2νv 2 T and taking into account that in the long wavelength approximation 1,2 (k) = c tr,l k, we again obtain (7). Compared to (7), the relationship (12) is valid for any value of k.

p-parameter
The relationship (7) (or (12)) defines the spectrum of fluctuations at any given . It is evident, however, that this parameter must be estimated for a given cluster directly from the confinement conditions [14,16]. There are different ways of establishing a necessary condition. Here we introduce a parameter Based on the results published earlier [14,16], we conclude that a lattice layer of finite size R is stably confined if roughly p const. It is worth noting that p is equivalent (but not identical) to the parameter α introduced in [14]. The parameter p turned out to be a weaker function on N , the number of particles in the cluster, and on a/λ, the interaction range, making it more advantageous for subsequent developments.
Here, r = r 1 − r 2 , b is the reciprocal lattice vector, and n = 2π (small terms ∝ c 2 tr /c 2 l omitted). It is assumed that r is large compared to the interparticle separation a.
Relationship (14) can be validated in the following manner. Note that in the case of = 0, the phonon energy (3) can be readily rewritten in matrix form as where δ is the unit matrix. By using this convenient representation, the averaged displacement cross-products, for example, can be calculated as Next it is straightforward to check that the correlation function can be represented through a sum running over the lattice Burgers vectors (for more details see [26]), Relationship (14) can be obtained in a standard way [26] by replacing the sum over k by an integral. This operation leads to where S is the area of the elementary cell, Sb 2 4π = 2π √ 3 and k max ∼ a −1 . Substitution of (18) back into (17) yields (14).

Strong confinement
In the experiments, = 0, i.e. it is always finite (although noticeably small, one or two orders of magnitude less than the frequency of the local 'caged' oscillations of the individual particles [19,27]). From (7) at non-vanishing , it immediately follows that the fluctuations remain finite at k → 0. Instead, the new singularities-the poles-appear at complex k. This absence of a singularity alters the character of disordering from algebraic (14) to exponential at a scale depending on the confinement parameter, It follows e.g. from the asymptotic behavior of the correlation function, where |u k | 2 as defined by relationship (7) is a sum of two terms. Actually, since c l c tr , the second term is most important, giving asymptotically Here, K 0 is the Macdonald function. The fact that the scale length (19) is determined dominantly by c l and not e.g. by c tr is actually not surprising, because (i) compression is at the origin of disordering and therefore any distortion of the lattice layer has to be scaled at a length dictated by the compressibility, and (ii) the shear modulus of the plasma crystals is typically much less than the compressibility modulus. It is essential that both asymptotes-algebraic and exponential-are treated as 'near-field' (r r conf ) and 'far-field' (r ∼ r conf ) approximations. Hence, it would be logical to assume that the ordering decay alternates with distance from algebraic to exponential. This is indeed in qualitative agreement with observations [17,18,33].
Remarkably, (3) and (7) are formally similar to the equations describing director fluctuations in nematic crystals in the presence of a magnetic field [26,34]. The action of the magnetic field is known to suppress the large-scale director fluctuations in liquid crystals.
There is another interesting interpretational aspect, which we discuss quantitatively below.

'One-plus' correlation length
The length scale ∼r conf seems to be of fundamental importance and hence the correlation length ∼r cor must be elucidated a bit differently. Let us take into account that the particles experiencing a horizontal confinement are distributed non-uniformly. The steady-state displacements of the particles u in the plasma crystal from their ideal locations in a uniform 2D crystal represent a growing function with distance u = r 3 /8r 2 conf [16]. Taking this into account and exploiting the Lindemann rule [35], we conclude that the lattice breaks up when (The density was obtained by computing the bounded Voronoi diagram using Delaunay triangulation, a standard algorithm provided by Mathematica 5.1. The density is normalized to n max = 1.63 mm −2 .) On average, the cluster density systematically decreases toward the edge. The local density maxima (minima) are easily identified as the positions of 5 fold (7 fold) cells. Note that the first circular row of the defects happens to appear at 0.46R from the center. This agrees with the theoretical estimate ∼(0.44-0.48)R obtained by using (23). Note also that the left part and the right part of the simulated cluster image are perfectly chiral (mirror-symmetric).
where L is the Lindemann parameter. Since L = 0.16-0.18 (see e.g. [36]), it follows that the first row of defects most probably appears at r r conf √ 8L (1.1-1.2)r conf . Making use of (13), (19), substituting r = r cor in (22) and replacing the inequality by the equality, one can estimate the size of the domains (or the equivalent correlation length) as We call the scale length (23) the 'one-plus' correlation length and refer to it as supplementary to the bond-order correlation length and translational correlation length, controlling ordering in homogeneous crystals [37]. The correlation length (23) does not depend explicitly on the temperature. In other words, for purely topological reasons, the big crystal spontaneously splits up into an array of subdomains, even at zero temperature. The estimated values of r cor agree well with those obtained in the simulation (see figures 1 and 2) and in experiments. For instance, it has been observed in  [13]. The radius of the defect-free region obtained in simulations with particles interacting with ∝ r −1 force [11,12] is shown for comparison (star). a laboratory experiment and a computational 'experiment' by Sheridan [18,19] that the crystal orientational order has a power-law decay at distances r/R < 1/4, which are in fairly good agreement with the value ≈0.3 obtained from (23). (This point is represented by the triangle in figure 2, which was calculated using parameters from [18,19].) The appearance of a network of growing domains as a 2D plasma crystal freezes was recently observed by Knapek et al [6] and by Hartmann et al [38].
The correlation length (23) depends only weakly on a/λ, on the interaction range and on N , the number of the particles inside the cluster. For instance, we observe that r corr /R ∝ (a/R) 1/2 ∝ N −1/4 because a/R ∝ N −1/2 according to [14]. The interaction-range dependence can also be explored. It is sufficient to note that relationship (23) is equivalent to the following ansatz, where α is the parameter introduced in [14]. At all other parameters fixed α ∝ (a/λ) 4 according to [14], and we conclude that In the Yukawa model, ξ depends only on the interaction range [21]; it is a rather weak function, showing, at a/λ = 0.5-2, approximately 17% variation.
It is worth validating also that the phonon dispersion effects are small, and hence the long-wavelength approximation we based our consideration on is well justified. It is important because usually the simulated clusters are limited in size: the ratio R/a ≡ the cluster size/interparticle spacing rarely exceeds 100 in simulations (see examples shown in figure 2). The dispersion of plasma crystal phonons is due to discreteness of the lattice and can be introduced by the parameter ς = (l disp /L) 2 , where l disp ∼ a is the dispersion length and L is the space scale of the wave field [21]. In our case, L = c l / , see (19), and the dispersion correction is small because For instance, there is R/a = 20, and expected ς ∼ 10 −3 for the cluster shown in figure 1. Note also that if necessary the dispersion corrections can be obtained more precisely with the help of relationship (12), which is valid at all wavelengths.

Relative position of defects
A symmetric situation like that shown in figure 1 is unlikely in general. Even after the multistep equilibration process, the final configuration of the cluster might be still far from the minimal energy state. The ground state might be difficult to reach practically due to, e.g., an enormously long equilibration time. Note that usually the ground state of a cluster is approached by annealing-repeating a heating-cooling procedure. It is a quite effective way to equilibrate the clusters consisting of a few tens of particles (see e.g. [13] and references therein). There are no 'typical' recommendations on how to achieve the ground state of a large multi-particle cluster with a few hundreds of particles or more. If the cluster is believed to be close enough to equilibrium, it is useful to estimate the RPD [13]-the relative position of defects effectively introducing the size of the defect-free central region by the relationship where r def is the mean position of defects, indicating the distance from the center of the cluster. A few cluster RPDs obtained in simulations are shown in figure 2.

Relation to the hexatic phase
Relationships (3) and (7), introducing the energy and the spectral intensity of the fluctuations, are universal in the sense that they do not depend explicitly on the dissipation rate. Hence, both relationships are equally applicable at any rate of friction, either strong or weak. The fluctuation-dissipation theorem implemented allowed us to prove this fact directly (see section 3.4). The necessity to confine the particles constituting the plasma crystal gives rise to a specific correlation length (23). The consequence of this 'one-plus' correlation length is that it unavoidably introduces a network of sub-domains to a confined lattice layer even at low temperatures. As a consequence (see section 4.3), high ordering is expected only at the shorter sub-domain scales, while at the larger super-domain scales the particle ordering decays exponentially, as our simple model predicts.
This might be of crucial importance e.g. for the interpretation of observations of the so called hexatic state in plasma crystals-still an outstanding and controversial issue in complex plasma studies [17]- [19].
An opportunity to observe the hexatic state [37], an intermediate state between the solid and liquid states, is of great interest. A hexatic phase has to appear in a 2D system with increasing temperature as the result of the second-order phase transition, and it is attributed to the clear difference in the character of translational ordering and orientational ordering.
We could not apply our findings directly to the study of the hexatic phase: the mean kinetic energy of the particles considered in our analytical model is in any case suggested to be small. Still, as it follows from (19)-(23), the confinement-induced self-ordering might be sensitive to the temperature increase only through softening of the elastic moduli. This is a rather weak effect for the plasma crystals, and a network of the 'confinement-induced' domains might exist at higher temperatures as well.
The existence of such sub-domains may eventually affect hexatic ordering due to the effective lowering of the global crystalline order. For instance, Sheridan explained the steeper decay of the bond-orientation order in the hexatic phase by clustering of dislocations observed in both the experiment [18] and the model [19]. It was suggested that the boundary effects stabilized the hexatic phase; otherwise it is only metastable.
For this reason it is of great importance to perform experiments and simulations with bigger clusters of microparticles, enabling the observation of bigger 'hexatic droplets'.
The logical development of this study is to discover a 'confined hexatic state', a natural generalization of the present results. We reserve this important but formidable subject for our later works.

Conclusion
We have demonstrated that the important features of the spatial structure of plasma crystals can be explained as the result of spontaneous disordering. Disordering is established by a fluctuation-activity-based self-organizing process. We characterize the spatial structure of the model crystals (with Yukawa-type interparticle interactions) by Fourier transform of the particle displacements and correlation functions. A probability-averaging algorithm is supplemented by a fluctuation-dissipation analysis. There are two new results: (i) an analytical approach was developed, which introduces a correlation length with respect to spontaneous 'cold' splitting of a crystal into an array of sub-domains, and (ii) a simple criterion for the appearance of such 'amorphization' is proposed. This could help us to explore and quantify the topology of microparticle clusters. The results were compared with the data obtained from numerical simulations for systems with strong and weak damping, respectively. Our analytical model describes the spontaneous crystal disordering observed recently in simulations and in experiments quite well-even quantitatively. Plasma crystals are perfectly suited to studying order-disorder transitions and symmetry alternations.