Two-dimensional Wigner crystals of classical Lennard-Jones particles

# Two-dimensional Wigner crystals of classical Lennard-Jones particles

Igor Travěnec and Ladislav Šamaj Institute of Physics, Slovak Academy of Sciences, Dúbravská cesta 9, 84511 Bratislava, Slovakia
###### Abstract

The ground-state of two-dimensional (2D) systems of classical particles interacting pairwisely by the generalized Lennard-Jones potential is studied. Taking the surface area per particle A as a free parameter and restricting oneself to periodic Bravais lattices with one particle per unit cell, Bétermin L [2018 Nonlinearity 31 3973] proved that the hexagonal, rhombic, square and rectangular structures minimize successively the interaction energy per particle as A increases. We show here that the second-order transitions between the rhombic/square and square/rectangular phases are of mean-field type. The aim of this paper is to extend Bétermin’s analysis to periodic 2D lattices with more than one particle per elementary cell. Being motivated by previous works dealing with other kinds of models, we propose as the ground-state the extensions of the 2D rectangular (1-chain) lattice, namely the “zig-zag” (2-chain), 3-chain, 4-chain, etc. structures possessing 2, 3, 4, etc. particles per unit cell, respectively. By using a recent technique of lattice summation we find for the standard Lennard-Jones potential that their ground-state energy per particle approaches systematically as the number of particles per unit cell increases to the one of a phase separated state (the optimal hexagonal lattice). We analyze analytically the low-density limit A\to\infty and the limiting hard-core case of the generalized Lennard-Jones potential.

###### pacs:
64.60.-i,61.50.Ah,64.70.K-,68.35.Rh
: J. Phys. A: Math. Gen.

Keywords: structural phase transitions, Wigner crystal, Lennard-Jones potential, hard disks

## 1 Introduction

Systems of classical pointlike particles interacting pairwisely via Lennard-Jones potential (LJ) [1] have been studied intensively. For two particles at distance r, the generalized version of this potential

is often referred to as the Mie potential. Here, n and m are positive integers, the inequality n>m prevents the collapse of particle pairs in the ground state. The first repulsive term is dominant at small distances, the second attractive term prevails at large distances between the particles. Since the particle distance with the minimum energy is given by r_{\min}=1 and V(1)=-1, the length and energy can be taken in arbitrary units. The potential (1) is denoted as the n-m LJ potential.

We will concentrate on the choice m=6, i.e. the n-6 model. This combination is considered very often, see e.g. [2], because 1/r^{6} is the attractive long-range van der Waals potential between molecules with electric dipole moments [3]. The original LJ value for the repulsive n=12 potential applies to argon [4], but some other rare gases have an effective value of n close to 11. For n finite, we speak about the soft potential since the particles can approach to one another arbitrarily close. In the n\to\infty limit, we get the hard-core potential of diameter 1 around each particle with an attractive interaction added,

Values of n up to n=36 were studied in Ref. [2]. Numerous other constraints for m and n were studied. The case n=2m seems to be motivated by the most common 12-6 potential rather than by a convincing physical situation. Nevertheless the 24-12 and 36-18 cases mimic, at least locally, colloidal systems during phase separation [5]. For systems in higher d-dimensions, the choice n=d+9 and m=d+3 was suggested for 3D gases. Sometimes the 9-3 potential is used for the solid-fluid interaction model [6]. The special case n=6k-3, m=3k-3 was found to be adequate for solids satisfying the spinodal condition [7]. The treatment of dynamics and pattern formation in systems with competing interactions or an external force can be found in Ref. [8].

As has been already mentioned, the potential (1) is attractive at large distances. Hence, at zero temperature, if we take a finite number of particles and gradually enlarge the space at their disposal, starting from some volume the particles will occupy only a restricted space region and the rest remains empty [9, 10]. In the two-dimensional (2D) space, particles create a phase separated state (the “optimal” hexagonal lattice); the hexagonal structure is known to minimize the energy for majority of purely repulsive pair interactions, ranging from the long-range Coulomb potential to hard disks. This phenomenon does not occur if external fields are applied, e.g. a periodic field with maxima at the vertices of the square lattice [11, 12, 13, 14], random pinning [14] or particles are restricted to a space domain by geometry, e.g. Wigner bilayers on parallel planes [15, 16].

Another possibility is to consider periodic structures with a fixed area A per particle (i.e., the particle density equals to 1/A) and to minimize the interaction energy per particle with respect to those structures. Bétermin [17] restricted himself to the standard 2D Bravais lattices with one particle per unit cell and applied this approach to the original 12-6 LJ potential. He has found that for small A, where the repulsive potential dominates at short distances, the hexagonal lattice provides the lowest energy. But as A rises, the rhombic, square and rectangular lattices become successively the energy minimizers.

It is natural to extend the ground-state problem to periodic 2D lattices with more than one particle per elementary cell. Our primary motivation comes from the work [12] about very distinct 2D systems of particles, interacting via the dipole 1/r^{3} potential and being in an external square-periodic potential. For such systems, the structural transition from the square to hexagonal lattices is realized via solitonic structures and the double-periodic 2-chain zig-zag lattice with two particles per unit cell. It is shown that this zig-zag phase provides a substantially lower ground-state energy than simple Bravais structures for intermediate and large values A. A further generalization comes from the study of quasi-one-dimensional systems of repulsive Yukawa particles in a confining potential [18]. In dependence on the parametric form of the potential, the particles form 2-, 3-, etc. chain configurations. We introduce 2D counterparts of these quasi-one-dimensional structures with respectively 2, 3, etc. particles per unit cell and find for the standard 12-6 LJ potential that their ground-state energy approaches systematically to the one of the optimal hexagonal lattice.

As a computational tool we apply a recent technique of lattice summations which provides a representation of the energy per particle as quickly converging series of Misra functions [15, 16]. In contrast to computer simulation methods, the technique requires ad-hoc lattice candidates for the ground-state. On the other hand, close to second-order phase transitions the method allows to derive the Landau expansion of the ground-state energy in the order parameter. The critical point is then given as the nullity equation for an expansion coefficient and can be specified with a prescribed accuracy by exploring a short computer time. Critical exponents can be obtained analytically. In the ground-state, critical exponents are usually of mean-field type, but there are exceptions from this rule [19]. In the context of Bétermin’s analysis [17] we show analytically that the second-order transitions between the rhombic/square and square/rectangular phases are of mean-field type. Furthermore, being motivated by previous works [12, 18], for the standard 12-6 LJ potential we consider the 2D zig-zag (2-chain), 3-chain and 4-chain structures which provide a lower ground-state energy than simple Bravais lattices for intermediate and large values of A. It is suggested that by increasing the number of particles in elementary cell the energy goes to the one of the optimal hexagonal lattice. Our technique of lattice summation is also very appropriate to study specific limiting cases where tiny energy differences among different structures are not accessible by simulation methods. In particular, for the zig-zag phase we analyze the n-6 potential up to the hard-disk limit n\to\infty or the small particle-density limit A\to\infty.

The paper is organized as follows. In section 2 we summarize basic facts about Bravais lattices, the double-periodic (2-chain) zig-zag, 3-chain and 4-chain structures and write down Misra series representations of their energies per particle. The phase diagram for the standard 12-6 LJ potential is constructed in section 3. The general n-6 case of the LJ potential with n going to large values, up to the limit of hard disks is the subject of section 4. A short recapitulation is presented in section 5.

## 2 Periodic structures and their energy

To calculate the energy per particle for periodic lattice structures, we use a recent method applied already to long-range Coulomb [15, 20] and Yukawa systems [16]; for those who would like to reproduce our present formulas, the cited papers describe the method in detail. The method involves the application of the gamma identity to each interaction term; note that each pair of particles shares the corresponding interaction energy which therefore must be taken with factor 1/2 when calculating the energy per particle. Subsequently, the lattice sum is expressed in terms of an integral over certain products of the Jacobi theta functions with zero argument [21]. Repeating twice the Poisson summation formula leads to the representation of the lattice sum as an infinite series of (generalized) Misra functions [22]. The series is converging very quickly; for instance, the truncation of the series at the 1st, 2nd, 3rd and 4th term reproduces the Coulomb Madelung constant for the hexagonal lattice up to 2,5,10,17 decimal digits, respectively. The computation of one energy value, e.g. by using MATHEMATICA, requires around one second of the CPU time on the standard PC. It turns out that the Misra series representation is technically simpler to derive for our potential (1) if n and m are even integers, and we shall restrict ourselves to this case.

We start with four Bravais lattices with one particle per elementary cell (see Fig. 1) which were treated in Ref. [17].

In the case of the rhombic lattice, the size of each side is equal to a, the smaller angle is \varphi and the area A=a^{2}\sin{\varphi}. We introduce the parameter \delta=\tan(\varphi/2) whose special values correspond to two rigid Bravais lattices. The choice \varphi=\pi/3 (i.e. \delta=1/\sqrt{3}), or equivalently \varphi=2\pi/3 (i.e. \delta=\sqrt{3}), corresponds to the hexagonal lattice, sometimes referred to as the equilateral triangular lattice; note that adjacent triangles form a rhomb with the area A=a^{2}\sqrt{3}/2. For \varphi=\pi/2 (i.e. \delta=1) one has the square lattice. The energy per particle of a rhombic structure with the n-m potential (1) is derived in the form

 E_{n,m}^{\rm rh}(\delta,A)=\frac{1}{n-m}\left[m\frac{v_{n}^{\rm rh}(\delta)}{A% ^{n/2}}-n\frac{v_{m}^{\rm rh}(\delta)}{A^{m/2}}\right], (3)

the Misra series representation of {v_{n}^{\rm rh}(\delta)} is given (for even values of n or m) in equation (A.4) of Appendix. Note the obvious symmetry \delta\to 1/\delta of the rhombic energy which corresponds to the equivalence of choosing either of the angles \varphi or \pi-\varphi.

Another structure of interest is the rectangular lattice with aspect ratio \Delta and A=a^{2}\Delta. The energy per particle of this lattice with the n-m potential (1) reads as

 E_{n,m}^{\rm rec}(\Delta,A)=\frac{1}{n-m}\left[m\frac{v_{n}^{\rm rec}(\Delta)}% {A^{n/2}}-n\frac{v_{m}^{\rm rec}(\Delta)}{A^{m/2}}\right]. (4)

The Misra series representation of {v_{n}^{\rm rec}(\Delta)} is given in equation (A.5) of Appendix. It stands to reason that this formula gives the same energy as (3) for the square lattice, i.e. for \Delta=\delta=1. Note the obvious symmetry \Delta\to 1/\Delta for the rectangular energy.

Besides simple Bravais lattices there exist other periodic structures with more particles per unit cell which can be potential energy minimizers. Here, we consider the successive more-chain extensions of the Bravais rectangular (1-chain) lattice pictured in the left panel of Fig. 2.

The zig-zag (2-chain) structure, adopted from Ref. [12], is presented with the notation of lattice parameters in the central panel of Fig. 2. The elementary cell, given by the rectangle defined by two vertical solid lines whose ends are connected by two horizontal dashed lines, contains two particles. Due to the lattice symmetry, all particles are equivalent in the sense that they see the same lattice arrangements of other particles. Since the density of particles is defined as 1/A, it holds that

 2A=\psi a^{2}. (5)

The length a^{\prime} is given by the relation

 \cos\varphi=\frac{\psi a}{2a^{\prime}}. (6)

The special case of the square lattice corresponds to \varphi=0, a=a^{\prime} and \psi=2. The hexagonal lattice is identified with the following parameters \varphi=\pi/6, a=a^{\prime} and \psi=2\cos{\varphi}=\sqrt{3}. Another possibility for the hexagonal lattice is the parameter choice \varphi=\pi/3, a^{\prime}=\psi a and \psi=1/\sqrt{3} which we shall consider in what follows. The energy per particle of the zig-zag structure reads as

 \displaystyle E^{\rm zz}_{n,m}(\varphi,\psi,A)=E_{n,m}^{\rm rec}(\psi,2A)+E_{n% ,m}^{\rm shift}\left(\psi,2A,\frac{1}{2},\frac{\psi}{2}\tan{\varphi}\right), (7)

where the term E_{n,m}^{\rm shift}(\psi,2A,c_{1},c_{2}) sums the energy contributions over an infinite rectangular lattice of sides a and \psi a, the sum being made with respect to a reference point shifted by the vector (c_{1}\psi a,c_{2}a) from a lattice point, see equations (A.6) and (A.7) of Appendix. For a given value of A, the energy (7) is minimized with respect to the free lattice variables \varphi and \psi, the lengths a and a^{\prime} are fixed by the relations (5) and (6).

The 3-chain structure, adopted from Ref. [18], is presented in the right panel of Fig. 2. The elementary cell, as before given by the central rectangle defined by two vertical solid lines whose ends are connected by two horizontal dashed lines, contains three particles. There are two non-equivalent sets of particles: the first set is given by two particles inside the unit cell while the second set consists in four particles forming the ends of the elementary rectangle (each weighted by the factor 1/4). Due to the reflection lattice symmetry, the number of lattice parameters is the same as in the zig-zag structure and they are defined analogously. The zig-zag relation (5) is modified to

 3A=\psi a^{2}, (8)

while the relation (6) remains unchanged. The energy of the 3-chain structure is obtained in the form

 \displaystyle E^{\rm 3ch}_{n,m}(\varphi,\psi,A) \displaystyle= \displaystyle E_{n,m}^{\rm rec}(\psi,3A)+\frac{4}{3}\ E_{n,m}^{\rm shift}\left% (\psi,3A,\frac{1}{2},\frac{\psi}{2}\tan{\varphi}\right) (9) \displaystyle+\frac{2}{3}\ E_{n,m}^{\rm shift}\left(\psi,3A,0,\psi\tan{\varphi% }\right).

We consider also a 4-chain structure (not pictured) with two non-equivalent sets of particles. For symmetry reason the distances between rows 1-2 and 3-4 must be the same and equal to (\psi/2)\tan\varphi_{1} while the distance between rows 2-3 (\psi/2)\tan\varphi_{2} may be different. The optimization of the energy for A=4 implies the angle values \varphi_{1}=0.33385\pi and \varphi_{2}=0.33328\pi which are very close to one another. To simplify numerical calculations, we set \varphi_{1}=\varphi_{2}=\varphi; the change in energies appears at the fourth decimal digit. Respecting that

 4A=\psi a^{2}, (10)

the energy of the 4-chain structure reads as

 \displaystyle E^{\rm 4ch}_{n,m}(\varphi,\psi,A) \displaystyle= \displaystyle E_{n,m}^{\rm rec}(\psi,4A)+E_{n,m}^{\rm shift}\left(\psi,4A,0,% \psi\tan{\varphi}\right) (11) \displaystyle+\frac{3}{2}E_{n,m}^{\rm shift}\left(\psi,4A,\frac{1}{2},\frac{% \psi}{2}\tan{\varphi}\right) \displaystyle+\frac{1}{2}E_{n,m}^{\rm shift}\left(\psi,4A,\frac{1}{2},\frac{3% \psi}{2}\tan{\varphi}\right).

One can proceed to simplified equidistant more-chain structures by considering only two independent variables \psi and \varphi, but increasing the number of chains the derivation of the energy per particle becomes more cumbersome.

## 3 Standard 12-6 LJ potential

In this section, we study the standard 12-6 case of the LJ potential (1).

### 3.1 Simple Bravais lattices

Restricting himself to simple Bravais lattices, Bétermin [17] showed that by increasing the value of A, the hexagonal, rhombic, square and rectangular lattices become successively energy minimizers. The transition between the hexagonal and rhombic lattices was found at A_{\rm BZ}\approx 1.1378475, the transition between the rhombic and square lattices at A_{1}\approx 1.1430032 and finally the transition between the square and rectangular lattices at A_{2}\approx 1.2679987. Our high-precision calculations give A_{\rm BZ}\approx 1.13784740849, A_{1}\approx 1.14300316307 and A_{2}\approx 1.26799868096, meaning a minor change in the last digit of A_{\rm BZ} and fully confirming the rest. The transitions at A_{1} and A_{2} are of second order, i.e. the corresponding parameters \delta and \Delta change continuously from the value 1 at the critical point. The change of \delta at A_{\rm BZ} is discontinuous, which means the first-order transition between the hexagonal and rhombic structures. The A-dependence of the ground-state energies for simple Bravais lattices is pictured in Fig. 3, a zoom around the short interval (A_{\rm BZ},A_{1}) where the rhombic structure becomes energy minimizer is given in Fig. 4.

As is seen in Fig. 3, the energy of the hexagonal lattice exhibits a minimum at A_{\rm opt} which corresponds to the previously defined optimal hexagonal lattice. Recalling that the hexagonal lattice is a special case of the rhombic one, namely E^{\rm hex}_{12,6}(A)=E_{12,6}^{\rm rh}(\delta=\sqrt{3},A), and using the representation (3), one has

 E^{\rm hex}_{12,6}(A)=\frac{v_{12}^{\rm rh}(\sqrt{3})}{A^{6}}-2\frac{v_{6}^{% \rm rh}(\sqrt{3})}{A^{3}}. (12)

The energy minimization with respect to A implies

 A_{\rm opt}=\left[\frac{v_{12}^{\rm rh}(\sqrt{3})}{v_{6}^{\rm rh}(\sqrt{3})}% \right]^{1/3}\approx 0.8491235647. (13)

Using the hexagonal relation A=a^{2}\sqrt{3}/2, the optimal lattice spacing a_{\rm opt}\approx 0.990193636287 is very close to 1 which is the distance between two isolated particles corresponding to the energy minimum. The optimal energy per particle

 E_{\rm opt}=E_{12,6}^{\rm rh}(\sqrt{3},A_{\rm opt})\approx-3.38212347861 (14)

is plotted in Fig. 3 as a horizontal dashed line. It applies for A>A_{\rm opt} where only a fraction of the available space is occupied by the optimal hexagonal lattice. For A\leq A_{\rm opt}, the space is fully occupied by the hexagonal lattice whose energy per particle depends on A (see the solid curve in Fig. 3).

Returning to periodic Bravais lattices, for very large A the rectangles become extremely narrow [17] which means that also the aspect ratio \Delta\to\infty, see the illustration in the left panel of Fig. 2. It can be readily shown that in that limit the most relevant terms from the Misra series representation of {v_{n}^{\rm rec}(\Delta)} in equation (A.5) correspond to the sum \sum_{j=1}^{\infty}z_{n/2+1}(j^{2}/\Delta)/(\frac{n}{2}-1)!. Consequently,

 \displaystyle E_{12,6}^{\rm rec}(\Delta,A) \displaystyle= \displaystyle\frac{v_{12}^{\rm rec}(\Delta)}{A^{6}}-2\frac{v_{6}^{\rm rec}(% \Delta)}{A^{3}} (15) \displaystyle\mathop{\sim}_{A\to\infty} \displaystyle\frac{1}{120A^{6}}\sum_{j=1}^{\infty}z_{7}\left(\frac{j^{2}}{% \Delta}\right)-\frac{1}{A^{3}}\sum_{j=1}^{\infty}z_{4}\left(\frac{j^{2}}{% \Delta}\right).

Since

 z_{4}(y)\mathop{\sim}_{y\to 0}\frac{2}{y^{3}}+O(1),\qquad z_{7}(y)\mathop{\sim% }_{y\to 0}\frac{120}{y^{6}}+O(1), (16)

we end up with

 E_{12,6}^{\rm rec}(\Delta,A)\mathop{\sim}_{A\to\infty}\left(\frac{\Delta}{A}% \right)^{6}\zeta(12)-2\left(\frac{\Delta}{A}\right)^{3}\zeta(6), (17)

where \zeta(s)=\sum_{j=1}^{\infty}1/j^{s} is the Riemann’s zeta function. Recalling that \Delta/A=1/a^{2} the energy per particle (17) corresponds to the one-dimensional array of particles with the 12-6 LJ potential; the point is that for \Delta\to\infty all other particles are at infinite distance from the reference one. For fixed A, the energy (17) as a function of \Delta is minimal for

 \Delta^{*}\mathop{\sim}_{A\to\infty}\left[\frac{\zeta(6)}{\zeta(12)}\right]^{1% /3}A=1.00566397A. (18)

The corresponding energy per particle is finite:

 E_{12,6}^{\rm rec}(\Delta^{*},A)\mathop{\sim}_{A\to\infty}-\frac{\zeta(6)^{2}}% {\zeta(12)}=-1.034732272. (19)

### 3.2 Second-order transitions between rhombic, square and rectangular phases

First we study the second-order transition between the rhombic and square phases. Parameterizing \delta=\exp({\epsilon}), the symmetry \delta\to 1/\delta of the energy (3) is equivalent to the transformation \epsilon\to-\epsilon and the energy is an even function of \epsilon. The Ginsburg-Landau form of its expansion around \epsilon=0 reads as

 E_{12,6}^{\rm rh}(\epsilon,A)=E_{12,6}^{\rm rh}(0,A)+g_{2}(A)\epsilon^{2}+g_{4% }(A)\epsilon^{4}+\ldots, (20)

where the Misra representation of the expansion functions g_{2}(A) and g_{4}(A) is at our disposal (we do not write them explicitly because formulas are too lengthy). The critical point A_{1} is identified with the condition g_{2}(A_{1})=0. The functions g_{2}(A) and g_{4}(A) behave in the vicinity of the critical point A_{1} as follows

 \displaystyle g_{2}(A) \displaystyle= \displaystyle g_{21}(A_{1}-A)+{\cal O}[(A_{1}-A)^{2}], \displaystyle g_{4}(A) \displaystyle= \displaystyle g_{40}+{\cal O}(A_{1}-A), (21)

where the constant prefactors g_{21}<0 and g_{40}>0. The minimum energy is reached at \epsilon^{*}\approx\delta^{*}-1 given by

 \frac{\partial}{\partial\epsilon}E_{12,6}^{\rm rh}(\epsilon,A)\big{|}_{% \epsilon=\epsilon^{*}}\approx 2g_{2}(A)\epsilon^{*}+4g_{4}(A)(\epsilon^{*})^{3% }=0. (22)

For A>A_{1}, there is only one trivial solution \epsilon^{*}=0 which corresponds to the square lattice. For A<A_{1}, besides the trivial solution \epsilon^{*}=0 one gets also two non-trivial conjugate solutions \pm\epsilon^{*} with

 \epsilon^{*}=\left(-\frac{g_{2}(A)}{2g_{4}(A)}\right)^{1/2}\approx\left(-\frac% {g_{21}}{2g_{40}}\right)^{1/2}\sqrt{A_{1}-A}, (23)

which provide a lower energy than the trivial solution. The order parameter \epsilon^{*}\propto\sqrt{A_{1}-A} is thus associated with the mean-field critical index \beta=1/2. The fitting of the numerical dependence of \delta^{*}-1 on A_{1}-A gives \beta=0.5003, confirming the mean-field character of the phase transition.

Analogously one can find the same type of transition between square and rectangular phases, using the energy expression (4) and the parameterization \Delta=\exp({\epsilon^{\prime}}) ensuring the symmetry \epsilon^{\prime}\to-\epsilon^{\prime} of the energy. The numerical value \beta=0.5006 confirms the same scenario and the mean-field character of this transition for A\to A_{2}^{+}.

### 3.3 Beyond simple Bravais lattices

Including the double-periodic zig-zag (2-chain) structure modifies the ground-state phase diagram significantly. As one can see in Fig. 3, the hexagonal lattice remains the energy minimizer for low values of A, namely 0<A<A_{t}=1.0650263. But for A>A_{t} that is the zig-zag structure which minimizes the energy and all other Bravais lattices become suppressed. The phase transition between the hexagonal and zig-zag structures is of first order, as the parameters \varphi and \psi undergo a stepwise change at A=A_{t}, namely from the hexagonal values \varphi=\pi/3=1.0471975511966 and \psi=1/\sqrt{3}=0.57735026918963 to the zig-zag values \varphi=1.03912932433 and \psi=0.5313649234.

Let us study the limiting case A\to\infty of the zig-zag phase. Numerical calculations up to A=5 reveal that A\psi\to{\rm const}\approx 0.429294, i.e. \psi\to 0 when A\to\infty. This corresponds to elongated zig-zag strips as counterparts of the elongated rectangles, see Fig. 2. In analogy with the elongated rectangles it can be shown that in the limit \psi\to 0 the leading terms of the zig-zag energy (7) correspond to the triangle row adjacent to the reference site:

 \displaystyle E^{\rm zz}_{12,6}(\varphi,\psi,A) \displaystyle\mathop{\sim}_{A\to\infty} \displaystyle\frac{\zeta(12)}{(2A\psi)^{6}}-\frac{2\zeta(6)}{(2A\psi)^{3}}-% \frac{1}{(2A\psi)^{3}}\sum_{j=-\infty}^{\infty}\frac{1}{(j^{2}+j+\frac{1}{4% \cos^{2}{\varphi}})^{3}} (24) \displaystyle+\frac{1}{2(2A\psi)^{6}}\sum_{j=-\infty}^{\infty}\frac{1}{(j^{2}+% j+\frac{1}{4\cos^{2}{\varphi}})^{6}}.

The first two terms on the r.h.s. sum over the line on which the reference site belongs, in analogy with the rectangular equation (17), with the substitution A\to 2A and \Delta\to 1/\psi. The last two terms sum over the upper line of the triangle array. They are given by the k=0 case of last sum in equation (A.7) with c_{1}=1/2 and c_{2}=(\psi\tan{\varphi})/2. The sums in equation (24) are slowly converging, but they can be summed up analytically by using a generating sum, see formulas (A.8)-(A.10) in Appendix. Minimizing the energy (24) with respect to \varphi and \psi yields

 \displaystyle\varphi^{*} \displaystyle\mathop{\sim}_{A\to\infty} \displaystyle 1.0500805037=0.33425100567\pi, \displaystyle\psi^{*} \displaystyle\mathop{\sim}_{A\to\infty} \displaystyle\frac{0.42929422368}{A}. (25)

Note that the angle \varphi^{*} is close to \pi/3 for the present 12-6 potential which means that the rows are composed of almost equilateral triangles. Within the calculations presented in the next section, for the n-6 potential with very large n=42 we get \varphi^{*}=0.33349\pi. Therefore, it is likely in the limit n\to\infty that \varphi^{*}=\pi/3 holds exactly. In other words, in the low-density limit A\to\infty the system seems to produce the zig-zag structure with distant rows of optimal equilateral triangles. The corresponding energy per particle

 E_{12,6}^{\rm zz}(\varphi^{*},\psi^{*},A)\mathop{\sim}_{A\to\infty}-2.1163 (26)

is well below the one for elongated rectangles (19).

As is clear from Fig. 3, the 3-chain and 4-chain structures become dominant at intermediate and large values of A and provide even lower energies than the zig-zag phase. The phase transitions to the hexagonal phase are of first order. Performing the large-A analysis for 3-chains and 4-chains, the corresponding energies

 E_{12,6}^{\rm 3ch}(A)\mathop{\sim}_{A\to\infty}-2.5279,\qquad E_{12,6}^{\rm 4% ch}(A)\mathop{\sim}_{A\to\infty}-2.7392 (27)

are lower than the one for the zig-zag structure, see also Fig. 3. By increasing the number n of chains to infinity one might expect that the energy goes to the optimal hexagonal one. Applying the [1/1] Padé extrapolation

 E_{12,6}^{\rm n-chain}(A\to\infty)=\frac{a_{0}+a_{1}n}{b_{0}+n} (28)

to n=2 (26) and n=3,4 (27) data, one obtains the value

 E_{12,6}^{\infty}=a_{1}=-3.3961 (29)

which is indeed very close to the optimal hexagonal energy (14).

## 4 General n-6 case

For the considered Bravais lattices, we analyze the n-6 interaction with even n ranging from 8 to 36. The results are presented graphically in Fig. 5. As concerns the transition between the hexagonal and rhombic structures at A_{\rm BZ} (diamonds), there was a possibility that starting from some n the first-order phase transition becomes of second order, but this is not the case. In the hard-core n\to\infty limit, A_{\rm BZ}(n) approaches to the close-packed hexagonal value A_{h}=\sqrt{3}/2\approx 0.8660254 and the rhombic parameter \delta goes to its hexagonal value 1/\sqrt{3}\approx 0.57735 which was chosen as the bottom value of the y-axis in Fig. 5. The rhombic-square transition value A_{1}(n) (squares) reaches asymptotically the closed-packed square value 1.

For any finite n, there is a first-order phase transition at A_{t}(n) between the hexagonal lattice, which is energy minimizer for A<A_{t}(n), and the zig-zag structure, which is energy minimizer for A>A_{t}(n), see circles in Fig. 6. In the limit n\to\infty, the close-packed hard disks form a hexagonal lattice with A=A_{h} (horizontal dashed line). There is no phase for A<A_{h} due to the hard-core restrictions among particles. The hexagonal phase, which exists only at A_{h}, converts immediately to the zig-zag phase which is the energy minimizer for all A>A_{h}. The corresponding transition values of the lattice parameters \varphi_{t} (squares) and \psi_{t} (triangles) go to their hexagonal values \varphi_{h}=\pi/3 and \psi_{h}=1/\sqrt{3}.

Now we analyze in detail the limit n\to\infty for the zig-zag structure to describe its transformation to the hexagonal lattice of hard disks with diameter 1 at A_{h}. In that limit, one keeps in the energy (7) only the attractive part of the interaction, whereas the repulsive part of interaction is replaced by a set of hard-core constraints. These constraints can be established by simple geometrical considerations, namely the three nearest neighbors of a chosen particle cannot be at distance smaller than 1 from that particle. This implies the conditions

 \displaystyle\psi \displaystyle\geq \displaystyle\frac{1}{2A}, \displaystyle\frac{\psi}{\cos{\varphi}^{2}} \displaystyle\geq \displaystyle\frac{2}{A}, \displaystyle\frac{\psi^{2}}{4}+\left(1-\frac{\psi}{2}\tan{\varphi}\right)^{2} \displaystyle\geq \displaystyle\frac{\psi}{2A}. (30)

With these constraints, the minimization of the attractive part of the zig-zag energy for values of A slightly above A_{h} showed that only the first constraint remains to be a strict inequality \psi>1/2A, while the other two become the equalities. Thus we can write down the set of two equations for three free parameters A, \varphi and \psi:

 \displaystyle\frac{\psi}{\cos^{2}\varphi} \displaystyle= \displaystyle\frac{2}{A}, \displaystyle\frac{\psi^{2}}{4}+\left(1-\frac{\psi}{2}\tan{\varphi}\right)^{2} \displaystyle= \displaystyle\frac{\psi}{2A}. (31)

Taking \varphi as an independent parameter, one gets

For \varphi=\varphi_{h}=\pi/3, one can check that we get the close-packed hexagonal values A_{h}=\sqrt{3}/2 and \psi_{h}=1/\sqrt{3}. Based on (32), for A close to A_{h} we get the expansions for small deviations of \varphi and \psi from their hexagonal values in powers of A-A_{h},

 \displaystyle\varphi-\varphi_{h} \displaystyle= \displaystyle A-A_{h}+\sqrt{3}(A-A_{h})^{2}-\frac{2}{3}(A-A_{h})^{3}+{\cal O}(% A-A_{h})^{4}, \displaystyle\psi-\psi_{h} \displaystyle= \displaystyle\frac{4}{3}(A-A_{h})+\frac{16}{3\sqrt{3}}(A-A_{h})^{2}+\frac{112}% {9}(A-A_{h})^{3}+{\cal O}(A-A_{h})^{4}.

These expansions were checked also numerically. The fact that the expansions are analytic in the deviation A-A_{h} implies that there is no phase transition in the limit of close-packed hard disks. We conclude that the first-order transition between the hexagonal and zig-zag phases, which exists for any finite n, disappears for n\to\infty.

A similar analysis can be made for 4-chain structures with qualitatively similar results.

## 5 Conclusion

In this paper, we have studied the 2D ground-state Wigner structures of classical particles with standard 12-6 and n-6 LJ pairwise interactions. We have used special lattice summation techniques [15, 16, 20] which permit one to express the energy per particle as a quickly converging series of Misra functions. The main aim was to extend Bétermin’s treatment of simple Bravais lattices [17] to periodic 2D lattices with more than one particle per elementary cell.

Fixing the area per particle A, we have introduced in section 2 four simple Bravais lattices (see Fig. 1) and the extensions of the rectangular (1-chain) lattice to double-periodic zig-zag (2-chain), 3-chain and 4-chain structures (Fig. 2).

Section 3 deals with the standard 12-6 LJ potential. We have confirmed the results of Bétermin [17] that as A increases the hexagonal, rhombic, square and rectangular phases become successively the energy minimizers. The second-order transitions between the rhombic/square and square/rectangular phases are shown to be of mean-field type. As is seen in Fig. 3, the inclusion of the zig-zag (2-chain), 3-chain and 4-chain structures modifies the phase diagram substantially. These phases become dominant at intermediate and large values of A and their phase transitions to the hexagonal phase are of first order. We anticipate that as the number of chains goes to infinity, this transition will be of second order. The limit A\to\infty can be treated analytically and we found that the energy of the zig-zag phase (26) is well below the rectangular one (19); the rectangular phase is the energy minimizer of simple Bravais lattices for A\to\infty. The corresponding results for the 3-chain and 4-chain structures (27) indicate a systematic decrease of the energy to the one of the phase separated state (the optimal hexagonal lattice) as the number of chains increases; the Padé extrapolation (28) of available n=2,3,4 data implies the result (29) which is indeed close to the optimal value (14).

The general n-6 (n being an even integer) LJ potential is the subject of section 4; here, we restrict ourselves to Bravais and zig-zag structures. In the case of simple Bravais lattices n ranges from 8 to 36 (Fig. 5), while the first-order transition between the hexagonal and zig-zag phases is studied for n ranging from 8 to 42 (Fig. 6). We studied analytically the interesting limit n\to\infty, where the hexagonal closed-packed lattice of hard disks exists only at A_{h}=\sqrt{3}/2 (there is no phase for A<A_{h} due to the hard-core restrictions among the particles) and the zig-zag phase dominates in the whole region A>A_{h}. It was shown that the parameters of the zig-zag phase are analytical functions of the deviation A-A_{h}, so its transformation to the hexagonal phase at A_{h} is smooth and there is no phase transition for hard disks. In other words, the first-order phase transition between the hexagonal and zig-zag structures, which exists for any finite value of n, disappears for n\to\infty.

The support received from VEGA Grant No. 2/0003/18 is acknowledged.

## Appendix

The first step of the method developed in [15, 16, 20] is the replacement of the interaction term of the type 1/r^{n} by using the gamma identity:

 \frac{1}{r^{n}}=\frac{1}{(n-1)!}\int_{0}^{\infty}{\rm d}t\ t^{n-1}\exp(-rt). (A.1)

A series of subsequent transformations leads to the representation of the energy per particle with the aid of Misra functions [22]:

For even powers n or m of r=\sqrt{r_{x}^{2}+r_{y}^{2}} one gets the Misra functions with an integer values of \nu; their representation in terms of special functions, suitable especially in MATHEMATICA, is easy to obtain. We cut the Misra series representation at the eighth term to get the precision of the energy more than 25 digits. In particular, for the case 12-6 one needs the following Misra functions

 \displaystyle z_{-1}(y) \displaystyle= \displaystyle\frac{1}{2}\left[\frac{e^{-\pi{y}}(1-\pi y)}{\pi^{2}}+y^{2}\Gamma% (0,\pi y)\right], \displaystyle z_{-4}(y) \displaystyle= \displaystyle\frac{{\rm e}^{-\pi{y}}}{120\pi^{5}}\left(24-6\pi y+2\pi^{2}y^{2}% -\pi^{3}y^{3}+\pi^{4}y^{4}\right) \displaystyle-\frac{y^{5}}{120}\Gamma(0,\pi y), \displaystyle z_{4}(y) \displaystyle= \displaystyle\frac{{\rm e}^{-\pi{y}}\left(2+2\pi y+\pi^{2}y^{2}\right)}{y^{3}}, \displaystyle z_{7}(y) \displaystyle= \displaystyle\frac{{\rm e}^{-\pi{y}}}{y^{6}}(120+120\pi y+60\pi^{2}y^{2}+20\pi% ^{3}y^{3} (A.3) \displaystyle+5\pi^{4}y^{4}+\pi^{5}y^{5}),

where \Gamma(x,y) stands for the incomplete Gamma function [21].

The function v_{n}^{\rm rh}(\delta) in the energy representation (3) for the rhombic lattice and v_{n}^{\rm rec}(\Delta) in the energy representation (4) for the rectangular lattice read (for even values of n)

 \displaystyle v_{n}^{\rm rh}(\delta) \displaystyle= \displaystyle\frac{1}{2^{n/2+1}(\frac{n}{2}-1)!}\Bigg{\{}4\pi^{n-1}\sum_{j=1}^% {\infty}\left[z_{2-n/2}\left(\frac{4j^{2}}{\delta}\right)+z_{2-n/2}\left(4j^{2% }\delta\right)\right] (A.4) \displaystyle+\frac{\left(\frac{n}{2}+1\right)\pi^{n/2}}{\frac{n}{2}\left(% \frac{n}{2}-1\right)}+4\sum_{j,k=1}^{\infty}z_{n/2+1}\left(\frac{j^{2}}{\delta% }+k^{2}\delta\right) \displaystyle+4\pi^{n-1}\sum_{j,k=1}^{\infty}\left[1+(-1)^{j}(-1)^{k}\right]z_% {2-n/2}\left(\frac{j^{2}}{\delta}+k^{2}\delta\right) \displaystyle+2\sum_{j=1}^{\infty}\left[z_{n/2+1}\left(\frac{j^{2}}{\delta}% \right)+z_{n/2+1}\left(j^{2}\delta\right)\right] \displaystyle+4\sum_{j,k=1}^{\infty}z_{n/2+1}\left[\frac{(j-1/2)^{2}}{\delta}+% (k-1/2)^{2}\delta\right]\Bigg{\}},
 \displaystyle v_{n}^{\rm rec}(\Delta) \displaystyle= \displaystyle\frac{1}{2(\frac{n}{2}-1)!}\Bigg{\{}2\pi^{n-1}\sum_{j=1}^{\infty}% \left[z_{2-n/2}\left(\frac{j^{2}}{\Delta}\right)+z_{2-n/2}\left(j^{2}\Delta% \right)\right] (A.5) \displaystyle+\frac{\pi^{n/2}}{\frac{n}{2}\left(\frac{n}{2}-1\right)}+4\pi^{n-% 1}\sum_{j,k=1}^{\infty}z_{2-n/2}\left(\frac{j^{2}}{\Delta}+k^{2}\Delta\right) \displaystyle+2\sum_{j=1}^{\infty}\left[z_{n/2+1}\left(\frac{j^{2}}{\Delta}% \right)+z_{n/2+1}\left(j^{2}\Delta\right)\right] \displaystyle+4\sum_{j,k=1}^{\infty}z_{n/2+1}\left(\frac{j^{2}}{\Delta}+k^{2}% \Delta\right)\Bigg{\}}.

Let us consider an infinite 2D rectangular lattice of sides \psi a and a, and a reference point shifted with respect to a lattice point by the vector (c_{1}\psi a,c_{2}a). Summing the potential (1) over all points of the rectangular lattice with respect to the reference point in close analogy with [16], we obtain the energy

 E_{n,m}^{\rm shift}\left(\psi,2A,c_{1},c_{2}\right)=\frac{1}{n-m}\left[m\frac{% \Sigma_{n/2}(\psi,c_{1},c_{2})}{(2A)^{n/2}}-n\frac{\Sigma_{m/2}(\psi,c_{1},c_{% 2})}{(2A)^{m/2}}\right], (A.6)

where, considering that n/2 or m/2 is an integer N,

 \displaystyle\Sigma_{N}(\psi,c_{1},c_{2}) \displaystyle= \displaystyle\frac{1}{2(N-1)!}\Bigg{\{}2\pi^{2N-1}\sum_{j=1}^{\infty}\bigg{[}% \cos(2\pi jc_{1})z_{2-N}\left(\frac{j^{2}}{\psi}\right) (A.7) \displaystyle+\cos(2\pi jc_{2})z_{2-N}\left(j^{2}\psi\right)\bigg{]}+\frac{\pi% ^{N}}{N-1} \displaystyle+4\pi^{2N-1}\sum_{j,k=1}^{\infty}\cos(2\pi jc_{1})\cos(2\pi kc_{2% })z_{2-N}\left(\frac{j^{2}}{\psi}+k^{2}\psi\right) \displaystyle+\sum_{j,k=-\infty}^{\infty}z_{N+1}\left[(j+c_{1})^{2}\psi+\frac{% (k+c_{2})^{2}}{\psi}\right]\Bigg{\}}.

To express analytically the sums in (24), one starts with the generating formula

 \sum_{j=-\infty}^{\infty}\frac{1}{j^{2}+j+\kappa}=2\pi\frac{\tan\bigg{(}{\frac% {\pi}{2}\sqrt{1-4\kappa}}\bigg{)}}{\sqrt{1-4\kappa}}. (A.8)

Taking consecutive derivatives with respect to \kappa one gets the required power of the denominator on the l.h.s. A minor complication is that after inserting \kappa=1/(4\cos{\varphi}^{2}) we get

 \sqrt{1-4\kappa}={\rm i}\tan{\varphi}, (A.9)

but applying the equality \tan({\rm i}x)={\rm i}\tanh x in the numerator of the r.h.s. of equation (A.8) we get real expressions as expected, e.g.

 \displaystyle\sum_{j=-\infty}^{\infty}\frac{1}{[j^{2}+j+1/(4\cos{\varphi}^{2})% ]^{3}}=-2\pi\cot^{3}{\varphi}\Bigg{\{}-6\cot^{2}{\varphi}\tanh\left({\frac{\pi% }{2}\tan{\varphi}}\right) \displaystyle\qquad+\pi\cosh^{-2}\left(\frac{\pi}{2}\tan{\varphi}\right)\left[% 3\cot{\varphi}+\pi\tanh\left({\frac{\pi}{2}\tan{\varphi}}\right)\right]\Bigg{% \}}. (A.10)

## References

• [1] Lennard-Jones J E and Devonshire A F 1937 Proc. R. Soc. A 163 53
• [2] Sousa J M G, Ferreira A L and Barroso M A 2012 J. Chem. Phys 136 174502
Delage-Santacreu S, Galliero G, Hoang H, Bazile J P, Bonet C and Fernandez J 2015 J. Chem. Phys 142, 174501
• [3] Milton K A 2001 Casimir Effect (Singapore: World Scientific)
• [4] Verlet L 1967 Phys. Rev. 159, 98
• [5] Lodge J F M and Heyes D M 1999 Molecular Simul. 23, 203
• [6] Smith L S and Lee L L 1979 J. Chem. Phys. 71, 4085
• [7] Sun J X 2005 J. Phys. Cond. Mat. 17, L103
• [8] Zhao H J, Misko V R and Peeters F M 2010 New J. of Phys. 2012 14, 063032; 2013 Phys. Rev. E 88, 022914
• [9] Gardner C S and Radin C 1979 J. Stat. Phys. 20, 719
• [10] Blanc X and Lewin M 2015 EMS Surveys in Math. Sci. 2, 255
• [11] Pogosov V W, Misko V R, Zhao H J and Peeters F M 2009 Phys. Rev. B 79, 014504
• [12] Gränz B, Korshunov S E, Geshkenbein V B and Blatter G 2016 Phys. Rev. B 94, 054110
• [13] Pokrovsky V L and Talapov A L 1980 Zh. Exp. Theor. Phys. 78, 134
• [14] Pogosov V W, Zhao H J, Misko V R and Peeters F M 2010 Phys. Rev. B 81, 024513
• [15] Šamaj L and Trizac E 2012 Europhys. Lett. 98, 36004; 2012 Phys. Rev. B 85, 205131
• [16] Travěnec I and Šamaj L 2015 Phys. Rev. E 92, 022306
• [17] Bétermin L and Zhang P 2015 Commun. Contemp. Math. 17, 1450049
Bétermin L 2016 SIAM J. Math. Anal. 48, 3236; 2018 Nonlinearity 31, 3973
• [18] Galván-Moya J E, Misko V R and Peeters F M 2010 Phys. Rev. B 90, 094111; 2015 Phys. Rev. B 92, 064112
• [19] Antlanger M, Kahl G, Mazars M, Šamaj L and Trizac E 2016 Phys. Rev. Lett. 117, 118002; 2018 J. Chem. Phys. 149, 244904
• [20] Travěnec I and Šamaj L 2016 Phys. Rev. B 93, 104110; 2017 J. Phys. A: Math. Theor. 50, 485001
• [21] Gradshteyn I S and Ryzhik I M 2000 Table of Integrals, Series, and Products, 6th edn. (Academic Press, London)
• [22] Misra R D 1940 Math. Proc. Cambridge Philos. Soc. 36, 173
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters