# Molecular dynamics study of a classical two-dimensional electron system: Positional and orientational orders

## 1 Introduction

More than 60 years ago, Wigner pointed out that an electron system will crystallize due to the Coulomb repulsion for low enough densities (Wigner crystallization) . Although quantum effects play an essential role in a degenerate electron system, the concept of Wigner crystallization can be generalized to a classical case where the Fermi energy is much smaller than the thermal energy. A classical two-dimensional (2D) electron system is wholly specified by the dimensionless coupling constant $`\mathrm{\Gamma }`$, the ratio of the Coulomb energy to the kinetic energy. Here $`\mathrm{\Gamma }(e^2/4\pi ϵa)/k_BT`$, where $`e`$ is the charge of an electron, $`ϵ`$ the dielectric constant of the substrate, $`a`$ the mean distance between electrons and $`T`$ the temperature. For $`\mathrm{\Gamma }1`$ the system will behave as a gas while for $`\mathrm{\Gamma }1`$ as a solid. Experimentally, Grimes and Adams  succeeded in observing a transition from a liquid to a triangular lattice in a classical 2D electron system on a liquid-helium surface around $`\mathrm{\Gamma }_c=137\pm 15`$, which is in good agreement with numerical simulations .

On the theoretical side, two conspicuous points have been known for 2D systems: (i) Mermin’s theorem dictates that no true long-range crystalline order is possible at finite $`T`$ in the thermodynamic limit . To be precise, the $`1/r`$ Coulomb interaction is too long ranged to apply Mermin’s arguments directly. Although there have been some theoretical attempts  to extend the theorem to the Coulomb case, no rigorous proof has been attained. (ii) A theory due to Kosterlitz, Thouless, Halperin, Nelson, and Young (KTHNY) predicts that a “hexatic” phase, characterized by a short-range positional order and a quasi-long-range orientational order, exists between a liquid phase and a solid phase . Because the KTHNY theory is based on various assumptions and approximations, its validity should be tested by numerical methods such as a molecular dynamics (MD) simulation. Several authors have applied numerical methods to classical 2D electron systems, but they arrived at different conclusions on the KTHNY prediction .

In order to address both of the above problems, the most direct way is to calculate the positional and the orientational correlation functions, which is exactly the motivation of the present study .

## 2 Numerical Method

A detailed description of the simulation is given elsewhere , so we only recapitulate it. We consider a rectangular area with a rigid uniform neutralizing positive background in periodic boundary conditions. The aspect ratio of the rectangle is taken to be $`L_y/L_x=2/\sqrt{3}`$, which can accommodate a perfect triangular lattice . The Ewald summation method is used to take care of the long-range nature of the $`1/r`$ interaction. We have employed Nosé-Hoover’s canonical MD method  to incorporate temperature accurately.

The system is cooled or heated across the transition with a simulated annealing method. The results presented here are for $`N=900`$ electrons with MD runs with 30,000–110,000 time-steps for each value of $`\mathrm{\Gamma }`$. The correlation functions and other quantities are calculated for the last $``$20,000 time-steps of each run.

Following Cha and Fertig , we define the positional and the orientational correlation functions from which we identify the ordering in each phase. The positional correlation function is defined by

$$C(r)\rho _𝐆^{}(𝐫)\rho _𝐆(\mathrm{𝟎})=\frac{{\displaystyle \underset{i,j}{}}\delta (r|𝐫_i𝐫_j|)\frac{1}{6}{\displaystyle \underset{𝐆}{}}e^{\mathrm{i}𝐆(𝐫_i𝐫_j)}}{{\displaystyle \underset{i,j}{}}\delta (r|𝐫_i𝐫_j|)},$$

where $`𝐆`$ is a reciprocal vector of the triangular lattice with the summation taken over the six $`𝐆`$’s that give the first peaks of the structure factor \[inset (a) of Fig. 1\]. The orientational correlation function is defined by

$$C_6(r)\psi _6^{}(𝐫)\psi _6(\mathrm{𝟎})=\frac{{\displaystyle \underset{i,j}{}}\delta (r|𝐫_i𝐫_j|)\psi _6^{}(𝐫_i)\psi _6(𝐫_j)}{{\displaystyle \underset{i,j}{}}\delta (r|𝐫_i𝐫_j|)},$$

where $`\psi _6(𝐫)\frac{1}{n_c}_\alpha ^{\mathrm{n}.\mathrm{n}.}e^{6\mathrm{i}\theta _\alpha (𝐫)}`$ with $`\theta _\alpha (𝐫)`$ being the angle of the vector connecting an electron at $`𝐫`$ and the $`\alpha `$-th nearest neighbor with respect to a fixed axis. The summation is taken over $`n_c`$ nearest neighbors that are determined by the Voronoi construction .

## 3 Results and Discussions

Let us first look at the positional and the orientational correlation functions in Fig. 1 for $`\mathrm{\Gamma }=200`$ and $`\mathrm{\Gamma }=160`$, for which the system is well in the solid phase. The positional correlation is seen to decay slowly, indicative of an algebraic (power-law) decay at large distances. The round-off in the correlation function around half of the system size should be an effect of the periodic boundary conditions. The algebraic decay of the positional correlation indicates that the 2D electron solid has only a quasi-long-range positional order. Thus we have obtained a numerical indication that Mermin’s theorem  remains applicable to the $`1/r`$ Coulomb interaction, which is consistent with the analytical but approximate results obtained in .

On the other hand, the orientational correlation rapidly approaches a constant, indicating a long-range orientational order. Therefore, while the 2D electron solid has no true long-range crystalline order, we can say that it has a topological order. From a snapshot of the configuration \[see inset (b) in Fig. 1\], we can see that a long-range orientational order is preserved since defects (5- or 7-fold disclinations, etc) tend to appear as dislocation (5-7 combination of disclinations) pairs, i.e., 5-7-5-7 disclination quartets that only disturb the orientational correlation locally. Here the coordination number is again determined from the Voronoi construction.

Now we move on to the orientational correlation function (inset of Fig. 2) around the crystallization, which is obtained by cooling the system from a liquid to a solid. In between a short-range orientational order for $`\mathrm{\Gamma }=120`$ and a long-range one for $`\mathrm{\Gamma }=140`$, the correlation appears to decay algebraically at $`\mathrm{\Gamma }=130`$ with an exponent approximately equal to unity, which deviates from the upper bound of $`1/4`$ predicted by KTHNY . However, numerical difficulties arising from finite-size and finite-time effects prevent us from drawing any definite conclusion on the existence of the hexatic phase. Namely, we cannot rule out the possibility that the power-law behavior is an artifact of insufficient equilibration. In fact, the solid phase persists down to $`\mathrm{\Gamma }=130`$ when the system is heated from a solid, which is understandable if the solid-hexatic and the hexatic-liquid transitions are of first and second order, respectively.

The KTHNY theory is based on a picture that the hexatic-liquid transition occurs through unbinding of disclination pairs. To examine if this is the case, we have calculated a defect-defect correlation function (Fig. 2), which we define as a distribution of 7-fold coordinated electrons with respect to a 5-fold coordinated electron. The correlation function exhibits no qualitative difference between $`\mathrm{\Gamma }=120`$ and $`\mathrm{\Gamma }=130`$, for which the disclinations are not tightly bound as compared with $`\mathrm{\Gamma }=140`$. If we look at a snapshot for $`\mathrm{\Gamma }=130`$ (Fig. 3), we see some domain structure as far as the present numerical simulation with finite-size and finite-time restrictions are concerned. A finite-size scaling will be an interesting future problem if the transition into the hexatic phase is of second order.

We have also looked at a dynamical property, i.e., the power spectral density of the velocity (Fig. 4), which is related to the vibrational density of states and corresponds to the Fourier transform of the velocity autocorrelation function via Wiener-Khinchin’s theorem. While the difference between the solid and the liquid phases appears in the spectrum around zero frequency, which is proportional to the diffusion constant (finite in the liquid or vanishingly small in the solid), we find a peak around the typical phonon frequency, which, curiously enough, persists even in the liquid phase. This indicates that the liquid has well-defined local configurations despite the short-range positional and orientational correlations.