Structural, dynamical and melting properties of two-dimensional clusters of complex plasmas

The static and dynamical properties of two-dimensional mesoscopic clusters of equally charged classical particles are investigated through the Monte Carlo simulation technique. The particles are confined by an external harmonic potential. The ground-state configuration and the position of the geometry induced defects are investigated as a function of the inter-particle interaction (Coulomb, dipole, logarithmic and screened Coulomb). The eigenmodes are investigated and the corresponding divergence and rotor are calculated which describe the `shearlike' and `compression-like' modes, respectively. The melting behaviour is found to be strongly influenced by the inter-particle interaction potential: a small cluster with a short-range interaction melts earlier than one with long-range interaction. The melting temperature is related to the energy barriers between the ground state and the metastable states. For larger clusters, the melting scenario changes and is strongly influenced by the location of the topological defects.


23.3
with the theoretical prediction [13]. The mode integrated spectrum shows two broad maxima which are found to be due to 'shearlike' and 'compression-like' modes [20,22].
In the first part of the present paper, the ground state and the metastable state configurations for the Yukawa system are obtained and compared with the results of a Coulomb, dipole and logarithmic interacting system. Topological defects which are induced by the symmetry of the confinement potential will be investigated in larger systems. The excitation properties, (i.e. the normal modes) are also discussed. In most papers, one concentrates on the study of the minima in the energy landscape while the barriers between those minima are left out of the discussion. Here the saddle point states, which are the critical states in the transition between the ground state and the metastable states, are studied, from which the potential barriers are obtained. In the second part of this paper, the Monte Carlo (MC) method is used to study in detail the order-disorder transition ('melting') of the cluster.
The paper is organized as follows. In section 2, we describe the model system and the numerical approach. Section 3 is devoted to the structural properties at zero temperature, and both small and large systems are discussed. The eigenmode spectrum of Yukawa systems is shown in section 4. The discussion on the phase transitions and energy barriers are given in section 5. Our results are summarized in section 6.

Numerical approach
The particles interact through a screened Coulomb potential which, depending on the physical parameters and the background plasma, can have different forms. The average interparticle potential is often assumed to be isotropic and purely repulsive and approximated by the Yukawa (Debye-Hückel) potential [2]. The Hamiltonian for such a system is given by where m is the mass of the particles, q is the particle charge, ω 0 is the radial confinement frequency and ε is the dielectric constant of the medium. r i denotes the position of the ith particle, and λ D is the screening length. As the range of the Yukawa system can vary from an extremely short-range hard-sphere-like potential for small λ D to the long-range one-component plasma for λ D = ∞, it is applicable to a wide variety of systems.
To exhibit the scaling of the system, we introduce the characteristic scales in the problem: r 0 = (2q 2 /m ω 2 0 ) 1/3 for the length, E 0 = (mω 2 0 q 4 /2 2 ) 1/3 for the energy and T 0 = (mω 2 0 q 4 /2 2 ) 1/3 k −1 B for temperature. After the scaling transformations (r → r/r 0 , E → E/E 0 , T → T /T 0 ), the Hamiltonian can be rewritten in a simple dimensionless form as [12] where κ = r 0 /λ D is the inverse dimensionless screening length. When κ = 0, the interaction between particles is a pure Coulomb potential. Here, we will also consider other dimensionless interactions as 1/r 3 ij and − ln r ij , for the cases of dipole and logarithmic inter-particle interaction, respectively.
To obtain the ground state, we employed the Newton optimization technique after the standard MC routine. This procedure was outlined and compared with the standard MC technique Table 1. The configurations for the particle packing sequences for N = 4-20 and different screening parameters κ. Here κ is varied from 0 to 8.0. in [13]. At zero temperature, the behaviour of the Yukawa system is completely characterized by the number of particles N and κ. To study the phase transitions in the Yukawa system, we use the standard Metropolis algorithm [23]. The static properties of the classical Yukawa system (κ = 0) at T = 0 are characterized by three dimensionless parameters: N, T and κ. The eigenmode frequencies are obtained from the eigenvalues of the dynamical matrix [13] E αβ,ij = ∂ 2 E ∂r α,i ∂r β,j .
Saddle point configurations are the key configurations for transitions between different stable states. The technique we used to find the saddle point is explained in more detail in [24]. Further, the geometric properties of the energy landscape can be obtained from which the potential barriers are derived.

Ground-state configurations
The configuration of the cluster is determined by two competing effects, namely the circular symmetry of the confinement potential and the triangular structure of the Wigner lattice. For a small cluster (with small number of particles N), particles are packed into concentric shells.
Depending on the total number of particles, the functional form of the confining potential and the mutual repulsion potential, a complicated structure with inner triangular cores surrounded by outer circular shells can be observed [19]. Table 1 shows the packing sequences of the ground state for different Yukawa parameters κ for N = 4-20 particles. The parameter κ was varied from 0 to 8.0. Comparing the configurations for κ = 0.125-0.5 with the pure Coulomb potential 23.5 configurations [16], one notices the similarity of the ground-state configuration. However, a slight increase of κ beyond this value can already change the ground-state configuration. For small N ∈ {4-9}, the parameter κ does not influence the packing of the particles since the confinement dictates the favourable configuration. Table 1 clearly shows that the structure of the Yukawa clusters with N < 10 is independent of the screening strength κ. Even when κ becomes as large as 8.0, the structure of the cluster is exactly the same as that for the case of κ = 0, only the size of the system is smaller. However, for N = 10 particles the ground-state configuration (2,8) will change into (3, 7) with increasing κ. This is not the case for N ∈ {11, 13-18, 20}, and these configurations appear to be very stable with respect to the screening parameter κ.
To summarize the results of table 1: the screening parameter κ can influence and change the configuration for a particular number of particles. With increasing κ the ground state changes from a shell-like structure into a hexagonal lattice, and the radius of the system decreases.
In the following discussion, we focus on the N = 19 system as an interesting example. The energy together with the ground-state configurations for different values of κ are plotted in the insets of figure 1. One can clearly see that the shell structure for small κ transfers into a hexagonal structure for increased κ. In figure 1, we plotted the energy of the N = 19 cluster as a function of κ for two states, i.e. with the (1,6,12) and (1,7,11) configuration. Notice that the energy decreases with increasing screening parameter κ because of the decreased radius of the system which decreases the confinement potential energy. When κ < 0.66, the configuration (1, 6, 12) is the ground-state configuration (see inset figure 1(a)), while for 0.66 < κ < 4.58, the configuration (1, 7, 11) (see inset figure 1(b)) becomes the ground state. However, for κ larger than 4.58, a re-entrant behaviour of the first ground state is observed (see inset figure 1(c)). Notice the perfect hexagonal structure for κ = 16 in agreement with [9]. From figure 1(d), one can see that the radius of the system monotonically decreases with increased κ. For a large value of κ the particles are located in the central region where the influence of the confinement is very small, and they form a hexagonal structure as in a 2D infinite system.
For larger Yukawa systems, we use defect theory to get a better understanding of the topological nature of the cluster structure. When the parameter κ is small, the configuration is similar to the one of a pure Coulomb system (κ = 0), and large clusters form a hexagonal lattice in the central region with a ring structure near the edge [12]. We make use of the Voronoi construction [25] to picture the cluster structure (for a detailed description see [20]) and to show the coordination number of each of the particles. Two kinds of defect appear in the Yukawa system, i.e. dislocations and disclinations. A dislocation is a pair of two disclinations, namely, a defect with fivefold (−) and a defect with sevenfold (+) coordination number [26]. When those geometry induced defects are present, the net topological charge N − − N + always equals six as already demonstrated in [10].
In figure 2, the Voronoi structures are shown as insets for 300 particles with different κ values. The location of the defects is strongly influenced by the screening parameter κ. One can see that the defects can move from the edge (figure 2(a)) to the transition region (between the ring structures and the central region, figure 2(b)) when κ changes from κ = 0.5 to 2.0. For κ = 2.0, the defects are situated at the six corners of a hexagon, each corner with a single net topological charge. This situation is similar to the case for pure Coulomb inter-particle interaction [20]. Beyond κ = 4.0 (see figure 2(c)) the defect region will move back from the transition area to the edge. Meanwhile the system is reduced in size and the number of particles in the outmost ring changes from 43 for κ = 0.25 to 40 for κ = 2.0, and to 41 for κ = 16.0. This result is a consequence of the balance between the two competing forces. The cluster patterns are determined by the need to balance the tendency to form a triangular lattice against the formation of a compact circular shape. When increasing κ, the interaction potential will change into a shorter-range potential, the defect region is pushed to the centre; meanwhile the system radius is reduced. Further increasing the screening parameter κ, the particles move further to the centre of the confinement potential reducing its effect, which will squeeze the defect region back to the edge of the system. The same behaviour follows from the number of particles in the outmost ring. This quantity can decrease firstly and increase again while increasing the κ value. In order to measure the relative position of the defects (RPD), we introduce the following quantity: where r − and r + denote the position of '−' and '+' defects, respectively, R is the radius of the system and N −(+) is the number of '−' ('+') defects. The RPD for 300 particles with different inter-particle interaction is shown in figure 2. For the purpose of comparison, the RPD of the long-range logarithmic potential, the dipole potential and the Coulomb potential between particles are also shown [28]. Notice that initially the positions of the defects move from the edge of the cluster towards the inner part of the cluster with decreasing range of the inter-particle interaction. A local minimum is reached for κ ≈ 0.125. In the range κ = 0 → 0.25 the RPD varies by at most 8%. A pronounced minimum in RPD is found for κ ≈ 2 where the defects are situated around (2/3)R. This figure nicely illustrates the effect of the interplay between the confinement potential and the inter-particle potential on the position of the geometry induced topological defects.

The eigenmode spectrum
In this section, we will discuss the excitation spectrum corresponding to the ground-state configuration of the Yukawa system. In the pure Coulomb interaction system, it is well known that there are three eigenfrequencies which are independent of N [13]: ω/ω 0 = 0, 1 and √ 3, which correspond to the rotation of the system as a whole, the centre of mass mode and the breathing mode, respectively. In the Yukawa system, the modes with frequency ω = 0 and ω/ω 0 = 1 do not depend on the interaction between the particles [27]. But the frequency of the breathing mode depends on the screening parameter κ and the number of particles. In figure 3 we show the eigenfrequencies for the N = 19 particles as a function of κ where the three previous eigenmodes are indicated by different symbols. The eigenfrequencies of a uniform rotation and the centre of mass mode are independent of the value of κ. The eigenfrequency of the breathing mode increases with increasing κ. This agrees with previous theoretical work [9] and recent experimental results [21]. Note also that the spectrum exhibits a discontinuous behaviour of the value of the eigenmodes as function of κ at κ = 0.66 and 4.58 where the ground-state configuration changes between the (1, 6, 12) and (1,7,11) configurations (see figure 1).
The lowest non-zero normal mode frequency (LNF) is a measure for the stability of the ground state. It tells us how easy (or difficult) it is to deform this state. In general, the LNF will increase when increasing the value of κ for the same configuration, since the radius of the system decreases with increasing κ, and the system becomes stiffer [9]. Notice that the LNF changes sharply when the configuration changes. For example, the LNF as a function of κ for N = 19 particles is shown in figure 4. When κ < 0.66, the LNF increases with increasing κ. Its groundstate configuration is (1,6,12) and the LNF eigenmode (see inset (a) of figure 4) corresponds to a vortex-antivortex excitation [13]. At the critical value κ = 0.66, the LNF suddenly drops almost to zero. When 0.66 < κ < 4.58 the (1,7,11) is the ground state configuration and the LNF slightly increases with increasing κ. The LNF excitation now corresponds to an intershell   (1,6,12) configuration is a so-called 'magic number' [13] configuration which is more stable against deformations than the (1, 7, 11) configuration. In the former the positions of the particles on the two rings are commensurable which makes their positions strongly inter-locked, which is not the case for the (1, 7, 11) configuration.
Note that it was previously found [13] that for small Coulomb systems, the lowest non-zero frequency mode corresponds to intershell rotation, whereas for larger numbers of particles it corresponds to the excitation of a vortex-antivortex pair [13]. This is not always so for the Yukawa system as we just found. For a fixed value of N and depending on the value of κ, inter-shell as well as vortex-antivortex pair excitations, and even the centre of mass mode, can be the LNF eigenmode.
The modes can have a shearlike or a compression-like character. The compressional and shear properties can be extracted from the divergence and rotor of the velocity field. Recently, we calculated the spatially resolved divergence and rotor of the different modes [20] for a Coulomb interacting 2D cluster. In this paper, we will associate a single number for the shearlike and compression-like character of the different modes by calculating the spatial average of the divergence ∇ · v and the vorticity ( ∇ × v) z of the velocity field, following the approach of [13].
The z component of the rotor ψ r (k) = e z · rot ψ(k) and the divergence ψ d (k) = divψ(k)  of the field of eigenvectors of mode k are where the values of ψ d,i (k), and ψ r,i (k) for the ith particle are given by Here, m and M denote the index and number of neighbouring particles of particle i, r m is the positional coordinate of a neighbouring particle and A i (k) is the eigenvector of particle i for mode k. Note that in contrast to [13] where only the nearest-neighbour particle was used instead of the sum in equations (6a), (6b), we sum over all the neighbours as was recently also done in [22]. This turns out to be the correct approach because the position and direction of the nearest neighbour is strongly influenced by small variations of the local environment of particle i. Note also that we calculate the squared average over all the particles because the simple spatial average is of course zero. For small Yukawa clusters (i.e. N < 40), the properties of the rotor and divergence were discussed in [22]. Here we extend those results to large clusters. In figure 5 we plot ψ d (k) and ψ r (k) as a function of the excitation frequency for N = 200 particles for different values of κ. In general, the lower eigenfrequency spectrum corresponds to a rotational type of excitations which are vortex-antivortex-like motions of the particles which lead to practically no density fluctuations. Note that ψ r (k) exhibits a peak with a shoulder at intermediate frequencies.
The shoulder and peak merge into a broad peak in the limit of large κ. In the second half of the spectrum the divergence ψ d (k), which corresponds to compression-like motion, is appreciably different from zero. Notice that for relatively small κ values there is a clear separation between the shear and the compressional modes. However, for large values of κ, the divergence increases more monotonically with frequency before it reaches a maximum. For smaller N the separation between rotational and compressional modes is less strict and there is a large range of intermediate frequencies where eigenmodes are partly compressional and partly rotational [22].

Melting and energy barriers
In this section we will add thermal energy to the particles and we investigate how the ordered state is lost with increasing temperature. Previously several groups used the molecular dynamics [18,29] simulation technique to investigate melting of a finite 2D system. Here we apply the MC method and relate the melting to the dynamical properties of the system. To characterize the melting behaviour of the cluster, we calculated the mean square displacement u 2 R which was introduced in [12] as

23.11
with a = 2R/ √ N the average distance between the particles. Figure 6 shows u 2 R as a function of the scaled temperature T /T 0 for N = 18, 19 and 300 for different κ. At low temperatures the particles exhibit harmonic oscillations around their T = 0 equilibrium position, and the oscillation amplitude increases linearly but slowly with temperature: the particles are well localized and still display an ordered structure. Melting occurs when u 2 R increases very 23.12 sharply with T . Because of the finite number of particles one has a region rather than a well defined melting temperature. After the melting point, the particles exhibit liquid-like behaviour. Following [30], we can 'define' a melting temperature at the point where u 2 R ≈ 0.10, which in fact is the Linderman melting criterion. The melting temperature as defined in this way is depicted in figure 7. Notice that the ground-state configuration for N = 19 changes from (1,6,12) to (1,7,11) in the same range of κ. This may influence the melting behaviour. To avoid this, we calculated the melting properties for a cluster with N = 18 particles (see figure 6(b)) which has the same structure (only different radius) for all κ. Again we found that the melting temperature decreases with increasing screening parameter κ. Therefore we conclude that the melting temperature of a small Yukawa cluster will decrease with increasing screening parameter κ. However, from figure 6(c), we see that this conclusion is not valid for a large system, in fact the opposite κ dependence is found as is clearly shown in figure 7.
In the first part, we mentioned that the system becomes stiffer with increasing κ. Therefore we expected a higher melting temperature with increasing κ which is, as just shown, not the case except for large clusters. To understand this unexpected melting behaviour, we investigated the energy landscape. It is well known that the number of saddle points increases rapidly with the particle number N which leads to a very complex energy landscape. Therefore, we considered first the case with six particles. We obtained the barrier height between the ground state and the metastable state. For the six-particle system, there are always two saddle point states, one ground state (1, 5) and one metastable state (6) for every value of κ. One of these two saddle point states is between the ground state and the metastable state. The other saddle point is between the same metastable state and another one which are simply rotated with respect to each other. The potential barriers and the melting temperature for N = 6 particles as functions of the parameter κ are plotted in figure 8. One can see that the melting temperature has a maximum in the range 1 < κ < 2. Notice that the potential barrier between the ground state and the metastable state is closely related to the melting temperature. The melting occurs through thermal excitation of one of the particles which hops from the centre of the system to the shell. We made a similar investigation for N = 8 particles which is a more complex system having several saddle points and the results are shown in the inset of figure 8. Notice that the lowest two potential barriers The melting scenario of a small laterally confined 2D system was found earlier [12] to be a two-step process. Upon increasing the temperature, first intershell rotation becomes possible where orientational order between adjacent shells is lost while retaining their internal order and the shell structure. At even higher temperatures, the growth of thermal fluctuations leads to radial diffusion between the shells, which finally destroys positional order. To characterize the relative angular intrashell and the relative angular intershell displacement, we consider the functions defined in [12], for the relative angular intrashell square deviation and the relative angular intershell square deviation with ϕ i the angle of particle i with respect to the x-axis; i1 indicates the nearest particle to particle i in the same shell, i2 refers to the nearest particle to particle i in the adjacent shell and ϕ 0 = 2π/N R is the angle between adjacent particles in a shell for a uniform shell with N R particles. u 2 a1 characterizes the angular motion in a particular shell, while u 2 a2 shows if the motion of the two adjacent shells is correlated or not.
A clear two-step melting process is also found for the Yukawa system. The radially dependent mean square displacement and the angular square deviation for N = 19 particles with different κ are plotted in figure 9. From figures 9(a)-(c), one can see that the radial melting always appears at higher temperature than the angular melting. In figure 9(a), the u 2 a1 of the inner ring and the outer ring are almost the same as u 2 a2 which is the relative angular intershell square deviation. This means that when the inner ring loses its order, the relative order is lost simultaneously, and the outer ring has almost the same melting temperature as the inner ring. For κ = 2.0 (see figure 9(b)), the u 2 a2 displacement starts for much lower temperature than the u 2 a1 . It should be noted that when κ = 2.0, the ground-state configuration for N = 19 is (1, 7, 11), which is more unstable against intershell rotation than the configuration (1, 6, 12). Therefore, the relative intershell order disappears for very low temperature, and the two shell structures still exist and they melt at higher temperature. The radial melting will start at even higher temperature. Figure 9(c) is quite similar to figure 9(a), due to the same structure of the groundstate configuration. Using the definition u 2 i = 0.1 as the melting criterion we found for the loss of angular intershell correlations the melting temperature T a2 /T 0 = 0.0042, 0.0012, 0.0033, and for the loss of intrashell correlation T a1 /T 0 = 0.0061, 0.0055, 0.0045 and for the radial melting T R /T 0 = 0.0097, 0.011, 0.0067 for κ = 0.5, 2.0 and 8.0, respectively.
In order to better describe the spatial dependence of the melting process in large clusters, we separate the large system, for example N = 300, into three regions as we did earlier in [20]: region I is comprised of the defect-free hexagonal centre, region II is a transition region between the two outer rings and the centre region and region III consists of the outermost two rings. The outer two rings have the same number of particles which varies between N R = 40 and 44 depending on the value of κ. We plot the radial dependent mean square displacement for the three different regions and the angular square deviation of the two outer rings for three different values of κ in figure 10. We find that the melting always starts from the defect region as in [20]. From figures 10(a)-(c) we find that for a large system the melting temperature obtained from u 2 a1 23.16 of the two outer rings is almost the same as the melting temperature obtained from u 2 a2 . This means that the outer rings are always locked which is different for small clusters. Comparing figures 10(a), (b) with 9, we see that in the case of a large system with small κ the radial and angular displacements start to increase rapidly at approximately the same temperature. Thus for large clusters with small κ value angular melting of the outer rings appears at the same temperature as the radial melting. But this is not the case for N = 300 and κ = 16 as shown in figure 10(c). The angular melting of the outer rings starts before the radial melting. The reason is that for large κ the defects are situated at the edge (the Voronoi structures are shown in figure 2). This results in a radial melting of the outer rings which starts before the radial melting of the transition region and the centre region. Moreover, the angular motion of the outer part is even easier to excite than the radial motion in the outer ring. This results in a lower angular than radial melting temperature at the edge. Thus also here the melting is a two-step process: first angular melting occurs then the radial melting.

Conclusion
We presented the configurational and melting properties of 2D finite Yukawa clusters. The results are obtained numerically using the MC simulation technique. For a small number of particles, the ground-state configuration consists of a ring structure which transforms into a hexagonal structure with increasing screening parameter κ. For large clusters, geometry induced defects appear and the location of the defects depends on κ. For logarithmic and strongly screened inter-particle interaction those defects occur near, or at the boundary of the cluster. For κ = 2 the defects are found deepest inside the cluster near (2/3)R.
The frequencies of the eigenmodes were investigated as functions of κ. In particular the lowest non-zero eigenfrequency which is a measure for the stability of the cluster was investigated. The shear and compression content of the modes was calculated and it was found that for large N the lower half of the eigenfrequency spectrum is mostly 'shearlike' modes while the upper half of the spectrum contains the 'compression-like' modes.
The melting behaviour is found to be strongly influenced by the screening parameter κ: a small cluster with a short-range interaction melts earlier than the one with long-range interaction. This was explained by the dependence of the energy barriers on κ. The two-step melting processes are also found for small Yukawa clusters. The large clusters have a different melting scenario which is strongly influenced by the location of the defect region. The angular and radial melting temperatures coincide. A two-step melting can appear again for large systems with a high value of κ, if the defect region is situated at the edge of the system.