Wigner crystals for a planar, equimolar binary mixture of classical, charged particles

We have investigated the ground state configurations of an equimolar, binary mixture of classical charged particles (with nominal charges $Q_1$ and $Q_2$) that condensate on a neutralizing plane. Using efficient Ewald summation techniques for the calculation of the ground state energies, we have identified the energetically most favourable ordered particle arrangements with the help of a highly reliable optimization tool based on ideas of evolutionary algorithms. Over a large range of charge ratios, $q = Q_2 / Q_1$, we identify six non-trivial ground states, some of which show a remarkable and unexpected structural complexity. For $0.59 \lesssim q<1$ the system undergoes a phase separation where the two charge species populate in a hexagonal arrangement spatially separated areas.


I. INTRODUCTION
The identification of the ordered ground state configurations of classical charged particles is known in literature as the Wigner problem 1 . In two dimensions and at vanishing temperature these charges form a hexagonal lattice [2][3][4] . In rather recent investigations, this problem has been extended: one example is the bilayer problem where the charges are confined between two planes, which are separated by a finite distance 5,6 . Due to the availability of closed, analytic expressions for the potential energy of this particular system, the complete set of its ground state configurations could be identified.
In the present contribution we return to the single layer problem and consider an equimolar, binary mixture of charged particles (with nominal values Q 1 and Q 2 ), that condensate at vanishing temperature on a neutralizing plane. For given values Q 1 and Q 2 the ordered equilibrium configurations of the particles at vanishing temperature are imposed by the requirement that the potential energy is minimized. With the help of efficient and highly accurate Ewald summation techniques 7 , the lattice sum of this system can be evaluated for any particle arrangement on an arbitrary, two-dimensional lattice. Employing suitable optimization techniques, the parameters of these lattices are then optimized in such a way as to minimize the lattice sum. In this contribution we have used an optimization tool that is based on ideas of evolutionary algorithms (EAs) 8,9 . Within this concept, any possible two-dimensional lattice is considered as an individual, to which a fitness value is assigned.
These individuals are then exposed on the computer to an artificial evolution: via creation and mutation operations a large number of individuals is produced; in the former procedure a pair of new individuals is created from a pair of parent individuals which are selected according to their fitness values. Along this evolution only the best, i.e., the fittest, individuals are expected to survive and are thus retained. Bearing in mind that we are looking for the individual (= structure) with the lowest lattice sum, we assign a high fitness value to an energetically favourable ordered structure. EA-based optimization algorithms have turned out to be highly efficient and reliable tools for identifying ordered equilibrium structures in a broad variety of condensed matter systems, in general 9-12 , and for quite a few two-dimensional systems, in particular [13][14][15][16][17] .
In total we have identified six non-trivial ground states. They are characterized by a broad variety of structural complexity, which is the result of an intricate competition between the interactions of the two charges. Introducing the charge ratio q as the only relevant parameter of the system (defined as q = Q 2 /Q 1 with 0 ≤ q ≤ 1 due to simple consideration and due to symmetry arguments) we identified two structures that show a remarkable stability over relatively large q-ranges: (i) one of them can be described via two intertwining, commensurate square sublattices, one of them populated by charges Q 1 , the other by charges Q 2 ; this structure shows, in addition, among all ground states the highest energy gain with respect to a suitably defined reference state; (ii) in the other ordered equilibrium configuration strongly distorted, but symmetric hexagonal tiles cover the entire two-dimensional space hosting in their interiors pairs of the weaker charges. For 0.59 q < 1 the system undergoes a phase separation where the two spatially separated phases are represented by hexagonal lattices populated by either species of charges. The remaining four non-trivial ground states are dominated by distorted, asymmetric hexagonal arrangements of charges Q 1 , hosting in their interior pairs and triplets of charges Q 2 . For the limiting values, i.e., q = 0 and q = 1, we obtain the expected hexagonal particle arrangements.
The paper is organized as follows. In the subsequent section we briefly introduce our model system. Section 3 is dedicated to the methods we have used: both the Ewald summation technique as well as our optimization tool (based on ideas of evolutionary algorithms) are briefly summarized; further we introduce a suitable state of reference for our energetic considerations. In section 4 we thoroughly discuss the results. The contribution is closed with concluding remarks.

II. MODEL
We consider an equimolar mixture of classical charges with the particles being confined to a planar (i.e., two-dimensional ) geometry. The point charges (with nominal values Q 1 and Q 2 ) are located at positions r i and r j and interact via an unscreened Coulomb interaction with r ij = |r i − r j |.
Since the total number density (i.e., number of particles per unit area), ρ, can be scaled out via the distances, its actual quantity is irrelevant for further considerations. In the equimolar case we obtain for the partial number densities ρ 1 = ρ 2 = ρ/2.
For convenience we introduce the parameter q = Q 2 /Q 1 , i.e., the ratio between the two types of charges. Since negative values of q lead to a divergent potential energy and taking into account the symmetry q ↔ 1/q we can restrict ourselves to the range 0 ≤ q ≤ 1; thus we assume that charge Q 1 is stronger than charge Q 2 . Note that we recover the classical Wigner problem for q = 1.
To compensate for the charges, we introduce a uniform, neutralizing background on the plane, specified by a charge density σ, which is given by

III. METHOD
The present contribution is dedicated to a complete identification of the ground state configurations of an equimolar mixture of point charges, i.e., the ordered equilibrium structures at vanishing temperature. Following the basic laws of thermodynamics, the particles will arrange under these conditions in an effort to minimize the corresponding thermodynamic potential. For our system (i.e., fixed particle number N and density ρ), we have to minimize the potential energy, which reduces at vanishing temperature to the lattice sum of the ordered particle configuration.
Among the numerous optimization schemes available in literature, we have opted for an optimization algorithm that is based on ideas of EAs 8,9 . Our choice is motivated by the fact that this strategy has turned out to be highly successful in related problems for a wide variety of soft matter systems, including in particular problems in two-dimensional geometries [13][14][15][16][17] . For a comprehensive presentation of this optimization algorithm and of the related computational and numerical details, we refer the reader to 9,18 .
The quantity that has to be minimized is the lattice sum of an ordered particle configuration. Taking into account the long-range character of the interactions (1), this quantity can most conveniently be calculated via Ewald sums 7 . For the separation of r-and k-space contributions, we have used the cutoff values r c = 15/ √ ρ and k c = 10 √ ρ, respectively; for the Ewald summation parameter we use α = 0.3. This set of numerical parameters guarantees a relative accuracy of 10 −5 for the evaluation of the internal energy.
In an effort to specify the ordered structures, we introduce for convenience a q-dependent reference energy. For the one component system (i.e., q = 1) the ground state energy (per particle) of point charges, arranged on a neutralizing plate in a hexagonal lattice, is given C M = 1.960515789 being the Madelung constant of this particular particle arrangement.
We extend this expression continuously to q ≤ 1 along the following lines: we imagine the system to be split up into two infinitely large regions, labeled γ = 1 and γ = 2, each of them hosting exclusively the respective charges, Q 1 and Q 2 , and each of them being locally charge neutral. The two regions share a common border. Introducing local number densities, ρ i , for species i in region γ (i = 1, 2 and γ = 1, 2), we arrive at the following relations: Together with (2), we obtain With these values for the local number densities and assuming that the charges will form hexagonal lattices in the respective regions, we obtain -with the help of equation (3) -the total energy per particle for this system In an effort to characterize the emerging ground state configurations, we have evaluated the two-dimensional orientational bond order parameters, Ψ 4 and Ψ 6 19 , defined via N i is the number of nearest neighbours of a tagged particle with index i and Θ j is the angle of the vector connecting this particle with particle j with respect to an arbitrary, but fixed orientation.

IV. RESULTS
In an effort to identify the complete set of ordered ground state configurations of our system specified in Section II, we have performed extensive EA-runs, taking into account up to 20 particles per species and per unit cell. For a given state point, up to 5,000 individuals were created. Calculations have been performed on a discrete q-grid with a spacing of ∆q = 0.01; thus in this contribution ∆q defines the accuracy in the location of the boundaries between ground states. The respective minimum energy configurations were retained as the ground state particle arrangements.
In Figure 1 we display the energy (per particle), E(q), of these ground state configurations and the energy difference, ∆E(q) = [E(q) − E 0 (q)], with respect to the reference energy, E 0 (q), as defined in equation (5). Note that over the entire q-range ∆E(q) is very small, i.e., less than 5 × 10 −3 . The fact that the differences between competing structures are so small is a fingerprint of the long-range nature of the Coulomb interaction. This function is now non-monotonous and shows kinks for particular q-values which can be associated with the structural changes. For convenience, the vertical broken lines in Figure   1 indicate the limits of stability of the six non-trivial identified ground states. The fact that these kinks are sometimes more or less pronounced is related to three issues: (i) the In the adjacent q-range, 0.05 q 0.09, charges Q 1 maintain their hexagonal arrangements (Structure 2, displayed in the right panel of Figure 2). As compared to Structure 1, In the relatively wide range of 0.10 q 0.25 the two species of charges arrange in two intertwining, commensurate square lattices (Structure 3, cf. Figure 3). It has to be emphasized that both sublattices remain perfect over the entire q-range of stability, i.e., Ψ 4 = 1 for 0.10 q 0.25. The structural stability of this particular ground state is also reflected by the fact that Structure 3 is characterized by the highest energy gain compared to the energy value of the reference structure, E 0 (q) (see Figure 1).
Around the value q 0.26, charges Q 1 form hexagonal structural units which are in their shape reminiscent of gems or diamonds. These six-particle rings form in a headto-tail arrangement parallel lanes: adjacent six-particle rings of neighbouring lanes share vertices, while the remaining edges form equilateral triangles. This ground state is referred to as Structure 4; it is depicted in the left panel of Figure 4 and the six-particle rings are highlighted. Each of the six-particle units hosts in its center an essentially equilateral triangle of charges Q 2 . The positions of the six surrounding charges Q 1 are imposed by the condition that the smallest distance of these charges from any of the three inner charges (Q 2 ) has the same value; this requirement induces the particular shape of the six-particle rings. The triangular arrangements that fill up the interstitial space are not populated by The ground state identified for 0.27 q 0.28 (denoted as Structure 5 and depicted in the right panel of Figure 4) differs only in one feature from Structure 4: the lanes, formed by the head-to-tail arrangements of the six-particle rings (highlighted in the corresponding panel) are now antiparallel. In this configuration the neighbouring rings of adjacent lanes share edges, which leads now to the formation of rhombic four-particle arrangements which are again void of Q 2 charges. For 0.59 q < 1.00, we find that E 0 (q) < E(q) (see Figure 1), indicating that the phase separated particle configuration is energetically more stable than any other ordered structure identified by our optimization tool. Due to the limitation in the number of particles per unit cell, the EA-based search for ground state configurations proposes -depending on the number of particles per cell -configurations with increasing complexity, all of them being characterized by an energy value E(q) that is larger than the corresponding value E 0 (q).
Thus the demixed state, formed by two separate hexagonal, ordered regions and each of them being populated by one species of charge, is the ground state in this q-range.
Finally, for q = 1, we recover the one-component hexagonal monolayer (not displayed).

V. CONCLUSIONS
In this contribution we have investigated the ground state configurations of an equimolar, binary mixture of classical charged particles (with nominal charges Q 1 and Q 2 ), that self-assemble on a neutralizing plane. Our investigations are based on reliable Ewald summation techniques which allow an efficient evaluation of the ground state energies (= lattice sums). With the help of reliable optimization tools, which are based on ideas of evolu-tionary algorithms, we are able to identify the ordered ground states of the system: by searching essentially among all possible two-dimensional lattices, this algorithm identifies for a given charge ratio q = Q 2 /Q 1 the energetically most favourable particle arrangement.
Apart from the expected, trivial hexagonal lattices for q = 0 and q = 1, we could identify for 0 < q 0.59 in total six ground state configurations.
Quite unexpectedly, these particle arrangements show a remarkable structural complexity which is the result of the energetic competition between the charge-charge interactions.
Throughout, a pronounced impact of the weaker charges, Q 2 , on the sublattices formed by the stronger charges, Q 1 , could be observed: this holds even if the corresponding q-values are rather small (i.e., q 0.05). Except for a purely square particle arrangement (which is stable over a remarkably large q-range and which shows the highest energy gain among all ground states with respect to a suitably defined reference state -see also the discussion about energies below), the ground states can be described on the basis of asymmetric, sometimes strongly distorted six-particle arrangements formed by charges Q 1 , which host in their interior simple two-or three-particle configurations of charges Q 2 . A deeper insight into the mechanisms that govern the formation of ground states is gained by introducing a phase separated reference state, where the two species of charges populate in hexagonal arrangements spatially separated areas. Comparing the energies of our ground state configurations, E(q), with the energy of this reference state, E 0 (q), we find that these two functions differ by less than 5 × 10 −3 , a fact that represents a characteristic fingerprint of the long-range Coulomb interactions. An analysis of this energy difference as a function of q reveals that transitions from one ground state to an adjacent one become visible as (more or less pronounced) kinks in this function. Based on these considerations we could show that for 0.59 q < 1 the energetically most favourable particle arrangement is the (ideal) phase separated state.

ACKNOWLEDGMENTS
This work was financially supported by the Austrian Research Fund (FWF) under Proj.