Twodimensional Wigner crystals of classical LennardJones particles
Abstract
The groundstate of twodimensional (2D) systems of classical particles interacting pairwisely by the generalized LennardJones 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 secondorder transitions between the rhombic/square and square/rectangular phases are of meanfield 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 groundstate the extensions of the 2D rectangular (1chain) lattice, namely the “zigzag” (2chain), 3chain, 4chain, 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 LennardJones potential that their groundstate 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 lowdensity limit A\to\infty and the limiting hardcore case of the generalized LennardJones potential.
pacs:
64.60.i,61.50.Ah,64.70.K,68.35.RhKeywords: structural phase transitions, Wigner crystal, LennardJones potential, hard disks
1 Introduction
Systems of classical pointlike particles interacting pairwisely via LennardJones potential (LJ) [1] have been studied intensively. For two particles at distance r, the generalized version of this potential
V(r)=\frac{1}{nm}\left(\frac{m}{r^{n}}\frac{n}{r^{m}}\right),\quad n>m  (1) 
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 nm LJ potential.
We will concentrate on the choice m=6, i.e. the n6 model. This combination is considered very often, see e.g. [2], because 1/r^{6} is the attractive longrange 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 hardcore potential of diameter 1 around each particle with an attractive interaction added,
V(r)=\cases{\infty\qquad r<1,\cr\frac{1}{r^{m}}\quad r\geq 1.}  (2) 
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 126 potential rather than by a convincing physical situation. Nevertheless the 2412 and 3618 cases mimic, at least locally, colloidal systems during phase separation [5]. For systems in higher ddimensions, the choice n=d+9 and m=d+3 was suggested for 3D gases. Sometimes the 93 potential is used for the solidfluid interaction model [6]. The special case n=6k3, m=3k3 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 twodimensional (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 longrange 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 126 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 groundstate 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 squareperiodic potential. For such systems, the structural transition from the square to hexagonal lattices is realized via solitonic structures and the doubleperiodic 2chain zigzag lattice with two particles per unit cell. It is shown that this zigzag phase provides a substantially lower groundstate energy than simple Bravais structures for intermediate and large values A. A further generalization comes from the study of quasionedimensional 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 quasionedimensional structures with respectively 2, 3, etc. particles per unit cell and find for the standard 126 LJ potential that their groundstate 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 adhoc lattice candidates for the groundstate. On the other hand, close to secondorder phase transitions the method allows to derive the Landau expansion of the groundstate 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 groundstate, critical exponents are usually of meanfield type, but there are exceptions from this rule [19]. In the context of Bétermin’s analysis [17] we show analytically that the secondorder transitions between the rhombic/square and square/rectangular phases are of meanfield type. Furthermore, being motivated by previous works [12, 18], for the standard 126 LJ potential we consider the 2D zigzag (2chain), 3chain and 4chain structures which provide a lower groundstate 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 zigzag phase we analyze the n6 potential up to the harddisk limit n\to\infty or the small particledensity limit A\to\infty.
The paper is organized as follows. In section 2 we summarize basic facts about Bravais lattices, the doubleperiodic (2chain) zigzag, 3chain and 4chain structures and write down Misra series representations of their energies per particle. The phase diagram for the standard 126 LJ potential is constructed in section 3. The general n6 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 longrange 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 nm potential (1) is derived in the form
E_{n,m}^{\rm rh}(\delta,A)=\frac{1}{nm}\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 nm potential (1) reads as
E_{n,m}^{\rm rec}(\Delta,A)=\frac{1}{nm}\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 morechain extensions of the Bravais rectangular (1chain) lattice pictured in the left panel of Fig. 2.
The zigzag (2chain) 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 zigzag 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 3chain 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 nonequivalent 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 zigzag structure and they are defined analogously. The zigzag relation (5) is modified to
3A=\psi a^{2},  (8) 
while the relation (6) remains unchanged. The energy of the 3chain 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 4chain structure (not pictured) with two nonequivalent sets of particles. For symmetry reason the distances between rows 12 and 34 must be the same and equal to (\psi/2)\tan\varphi_{1} while the distance between rows 23 (\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 4chain 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 morechain 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 126 LJ potential
In this section, we study the standard 126 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 highprecision 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 firstorder transition between the hexagonal and rhombic structures. The Adependence of the groundstate 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})\approx3.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 onedimensional array of particles with the 126 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 Secondorder transitions between rhombic, square and rectangular phases
First we study the secondorder 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 GinsburgLandau 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 nontrivial 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 meanfield critical index \beta=1/2. The fitting of the numerical dependence of \delta^{*}1 on A_{1}A gives \beta=0.5003, confirming the meanfield 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 meanfield character of this transition for A\to A_{2}^{+}.
3.3 Beyond simple Bravais lattices
Including the doubleperiodic zigzag (2chain) structure modifies the groundstate 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 zigzag structure which minimizes the energy and all other Bravais lattices become suppressed. The phase transition between the hexagonal and zigzag 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 zigzag values \varphi=1.03912932433 and \psi=0.5313649234.
Let us study the limiting case A\to\infty of the zigzag 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 zigzag 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 zigzag 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 126 potential which means that the rows are composed of almost equilateral triangles. Within the calculations presented in the next section, for the n6 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 lowdensity limit A\to\infty the system seems to produce the zigzag 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 3chain and 4chain structures become dominant at intermediate and large values of A and provide even lower energies than the zigzag phase. The phase transitions to the hexagonal phase are of first order. Performing the largeA analysis for 3chains and 4chains, 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 zigzag 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 nchain}(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 n6 case
For the considered Bravais lattices, we analyze the n6 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 firstorder phase transition becomes of second order, but this is not the case. In the hardcore n\to\infty limit, A_{\rm BZ}(n) approaches to the closepacked 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 yaxis in Fig. 5. The rhombicsquare transition value A_{1}(n) (squares) reaches asymptotically the closedpacked square value 1.
For any finite n, there is a firstorder phase transition at A_{t}(n) between the hexagonal lattice, which is energy minimizer for A<A_{t}(n), and the zigzag structure, which is energy minimizer for A>A_{t}(n), see circles in Fig. 6. In the limit n\to\infty, the closepacked hard disks form a hexagonal lattice with A=A_{h} (horizontal dashed line). There is no phase for A<A_{h} due to the hardcore restrictions among particles. The hexagonal phase, which exists only at A_{h}, converts immediately to the zigzag 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 zigzag 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 hardcore 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 zigzag 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
A=\sin{2\varphi},\qquad\psi=\cot{\varphi}.  (32) 
For \varphi=\varphi_{h}=\pi/3, one can check that we get the closepacked 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 AA_{h},
\displaystyle\varphi\varphi_{h}  \displaystyle=  \displaystyle AA_{h}+\sqrt{3}(AA_{h})^{2}\frac{2}{3}(AA_{h})^{3}+{\cal O}(% AA_{h})^{4},  
\displaystyle\psi\psi_{h}  \displaystyle=  \displaystyle\frac{4}{3}(AA_{h})+\frac{16}{3\sqrt{3}}(AA_{h})^{2}+\frac{112}% {9}(AA_{h})^{3}+{\cal O}(AA_{h})^{4}. 
These expansions were checked also numerically. The fact that the expansions are analytic in the deviation AA_{h} implies that there is no phase transition in the limit of closepacked hard disks. We conclude that the firstorder transition between the hexagonal and zigzag phases, which exists for any finite n, disappears for n\to\infty.
A similar analysis can be made for 4chain structures with qualitatively similar results.
5 Conclusion
In this paper, we have studied the 2D groundstate Wigner structures of classical particles with standard 126 and n6 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 (1chain) lattice to doubleperiodic zigzag (2chain), 3chain and 4chain structures (Fig. 2).
Section 3 deals with the standard 126 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 secondorder transitions between the rhombic/square and square/rectangular phases are shown to be of meanfield type. As is seen in Fig. 3, the inclusion of the zigzag (2chain), 3chain and 4chain 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 zigzag 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 3chain and 4chain 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 n6 (n being an even integer) LJ potential is the subject of section 4; here, we restrict ourselves to Bravais and zigzag structures. In the case of simple Bravais lattices n ranges from 8 to 36 (Fig. 5), while the firstorder transition between the hexagonal and zigzag phases is studied for n ranging from 8 to 42 (Fig. 6). We studied analytically the interesting limit n\to\infty, where the hexagonal closedpacked lattice of hard disks exists only at A_{h}=\sqrt{3}/2 (there is no phase for A<A_{h} due to the hardcore restrictions among the particles) and the zigzag phase dominates in the whole region A>A_{h}. It was shown that the parameters of the zigzag phase are analytical functions of the deviation AA_{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 firstorder phase transition between the hexagonal and zigzag structures, which exists for any finite value of n, disappears for n\to\infty.
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}{(n1)!}\int_{0}^{\infty}{\rm d}t\ t^{n1}\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]:
z_{\nu}(y)=\int_{0}^{1/\pi}\frac{{\rm d}t}{t^{\nu}}\exp\left(\frac{y}{t}% \right),\qquad y>0.  (A.2) 
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 126 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(246\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^{n1}\sum_{j=1}^% {\infty}\left[z_{2n/2}\left(\frac{4j^{2}}{\delta}\right)+z_{2n/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^{n1}\sum_{j,k=1}^{\infty}\left[1+(1)^{j}(1)^{k}\right]z_% {2n/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{(j1/2)^{2}}{\delta}+% (k1/2)^{2}\delta\right]\Bigg{\}}, 
\displaystyle v_{n}^{\rm rec}(\Delta)  \displaystyle=  \displaystyle\frac{1}{2(\frac{n}{2}1)!}\Bigg{\{}2\pi^{n1}\sum_{j=1}^{\infty}% \left[z_{2n/2}\left(\frac{j^{2}}{\Delta}\right)+z_{2n/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_{2n/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}{nm}\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(N1)!}\Bigg{\{}2\pi^{2N1}\sum_{j=1}^{\infty}\bigg{[}% \cos(2\pi jc_{1})z_{2N}\left(\frac{j^{2}}{\psi}\right)  (A.7)  
\displaystyle+\cos(2\pi jc_{2})z_{2N}\left(j^{2}\psi\right)\bigg{]}+\frac{\pi% ^{N}}{N1}  
\displaystyle+4\pi^{2N1}\sum_{j,k=1}^{\infty}\cos(2\pi jc_{1})\cos(2\pi kc_{2% })z_{2N}\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{14\kappa}}\bigg{)}}{\sqrt{14\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{14\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
References
 [1] LennardJones 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
DelageSantacreu 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ánMoya 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