The rich polymorphic behavior of Wigner bilayers

Self-assembly into target structures is an efficient material design strategy. Combining analytical calculations and computational techniques of evolutionary and Monte Carlo types, we report about a remarkable structural variability of Wigner bilayer ground states, when charges are confined between parallel charged plates. Changing the inter-layer separation, or the plate charge asymmetry, a cascade of ordered patterns emerges. At variance with the symmetric case phenomenology, the competition between commensurability features and charge neutralization leads to long range attraction, appearance of macroscopic charges, exotic phases, and non conventional phase transitions with distinct critical indices, offering the possibility of a subtle, but precise and convenient control over patterns.

The self-assembly of colloidal systems opens the way to the synthesis of materials that considerably widen the class of known natural crystals, among which opals or butterfly wings. From an academic perspective, these new, complex structures allow for detailed and original studies of fundamental processes like nucleation, glass transition or low dimensional statistical physics [1,2]. Skillfully combined with progress in particle synthesis, self-assembly has led to a wealth of applications such as patterned magnetic systems or bandgap materials used in displays, optical devices, photochemistry and biological sensors [3][4][5][6]. Taking advantage of targeted selforganization requires a fine control of interactions between the entities under study. This tailoring is achieved in most cases either by a) introducing some patchiness on the colloids, b) increasing the complexity of the problem, by considering e.g. mixtures instead of pure systems, c) affecting the solvent through various additives (polymeric, electrolytic, etc.), d) introducing an external field, be it electric, magnetic, laser-optical or stemming from the interactions with a patterned substrate [7][8][9][10][11][12]. As fruitful as they have turned out to be, these strategies in general do not allow for convenient in situ changes of the obtained ordered structures, so that it is challenging to probe and tune their variety in a simple fashion, by controlling an external parameter. Relinquishing the four routes above, we consider here a pure classical system of charges and show that the simplest form of external control -confinement in a slab -induces an unexpected structural variability, which in turn opens the way for a precise structural control. We shall focus on energy minimizing configurations, relevant when the kinetic energy is small compared to the Coulombic potential energy, and where the charges are forced into a bilayer configuration, thereby creating a particular realization of a so-called Wigner crystal.
Wigner crystals were first predicted by the eponymous physicist in the 1930s for electrons in a metal [13], where they have actually never been observed. Instead, their occurrence has been reported in the 1970s for electrons at helium interfaces [14]. Found in neutron stars and in the interior of white dwarfs, they have subsequently been evidenced in semi-conductors [15][16][17][18], graphene [19], quantum dots, trapped ionic plasmas or other dusty plasmas [20], and in the colloidal realm [1]. While the symmetric setup is now completely understood [21][22][23][24][25], very little is known for the asymmetric one [26]. How do charges organize into the bilayer structures that spontaneously form in our problem? Upon answering this question, we will treat the general asymmetric situation (no mirror symmetry between the two layers), which turns out to be considerably more complex than the symmetric case.
The question addressed is the following. Consider an ensemble of mobile point charges ("particles") interacting via a 1/r potential, confined between two parallel plates bearing uniform charge densities σ 1 e and σ 2 e, with −e the (elementary) charge of the mobile particles; the system as a whole is electroneutral. What is the energyminimizing arrangement of particles for a fixed plate-toplate distance d? The Earnshaw theorem [27] provides a first clue: given that a classical system of point charges under the action of direct (i.e., not image) electrostatic forces alone cannot be in an equilibrium configuration, the particles are expelled from the slab interior, and have to stick to the confining plates. Numerical and analytical work have furthermore shown that when σ 1 = σ 2 , staggered configurations arise on each plate, which -depending on d -can be rectangular/square, rhombic or hexagonal [21][22][23][24][25], see also the line A = σ 2 /σ 1 = 1 in Fig. 1.
Our interest focuses on the asymmetric case (σ 1 = σ 2 ), first on the cornucopia of ordered structures that appear as energy-minimizing, and second, on the distinct properties that characterize these new phases. Indeed, macroscopic charges do emerge, which result in a long range attraction between the plates. In addition, different universality classes are probed by changing solely the interplate separation, and overcharging is reported in some pocket of the phase diagram. By a combination of complementary analytical and computational techniques of evolutionary and Monte Carlo type [40], our goal is to unravel these properties, while charting out the phase diagram.
Without loss of generality, we assume σ 1 > 0. We introduce the asymmetry parameter A = σ 2 /σ 1 , and consider A ∈ [0, 1] [28]. We then define the dimensionless distance η = d (σ 1 + σ 2 )/2. Our system is thus entirely specified by η and A. We further introduce the surface particle densities n 1 and n 2 and the order parameter x = n 2 /(n 1 +n 2 ). Electroneutrality imposes σ 1 +σ 2 = n 1 +n 2 . In general n i = σ i (i = 1, 2) and thus each of the plates as a whole (i.e. particles plus surface charge density) is charged. φ(z) = −2πe(σ 1 − σ 2 )z, 0 < z < d. If local neutrality holds for both plates, then n i = σ i (i = 1, 2) and we find x = x neutr ≡ A/(1 + A). This should be the case when d → ∞, since violating local neutrality would result in a macroscopic electric field at large d, with divergent energy.
Upon changing η in the symmetric case (A = 1), it is known that a sequence of five phases (denoted I to V) emerges, consisting of two equivalent, "ideal" (i.e. undistorted) structures on plates 1 and 2, shifted with respect to one another. For A = 1, each of the plates is locally neutralized, and x = x neutr = 1/2. Increasing η, the hexagonal Wigner monolayer (phase I) is found at η = 0, then a bilayer with rectangular arrangements on each plate (structure II), which transforms into a square lattice (structure III). A staggered rhombic arrangement (phase IV) and a staggered hexagonal lattice (structure V) are subsequently observed, see Fig. 1 [25]. All transitions are continuous, except IV → V, of first order [25].
The analytical work proceeds with the derivation of new series representations for the Coulombic energies of undistorted structures [40]. This yields the exact Coulombic energy of the structures considered. On the other hand, the numerical work is two-pronged: a first technique, inspired from evolutionary algorithms (EA), identifies the optimal periodic structures among all those that have less than 40 particles per unit cell [40]; the second line of attack consists in extensive Monte Carlo simulations [29,40] on much lager systems (∼ 4000 particles per unit cell). To quantify order and identify the complex patterns formed, it is indispensable to introduce suitable probes: besides the population index x, we have used the bond orientational order parameters of symmetry n = 4, 5, 6, 7, 8, 10, 12, 18 and 24, as defined from a Voronoi construction [40]. These parameters Ψ (L) n take unit value under perfect n-fold ordering, and have been computed in four different variants: restricting to particles in layer 1 (a choice referred to with index L = 1), in layer 2 (L = 2), by projecting layer 1 perpendicularly onto layer 2 (L = 3), or finally by studying the neighbors in layer 1 of a given particle in layer 2 (index L = 4). While all four variants lead to compatible results, it appears that the latter choice, Ψ (4) n , with n = 4, 5 and 6 is particularly conducive to investigating the phase behavior. These parameters have been used to construct the phase diagram in Fig. 1, in conjunction with a precise RGB scheme [32]. Fig. 1 gathers the results for about 35000 state points obtained with EA computations. The computational cost for MC simulations is about 200 times higher than for EA, due to the complex treatment of long ranged interactions in quasi-2D systems (see Ref. [30] and references therein). Consequently, a smaller number of state points can be explored with MC simulation and a careful selection of those has to be operated to optimize resources [40]. MC simulations show that the structures obtained following the EA route are stable; more generally a complete agreement EA-MC is reported [31].
The first noticeable feature revealed by Fig. 1 is that at small distance η, it is always favorable for the charges to stick to the plate of highest surface charge (i.e. plate 1). Thus, in the black region of Fig. 1, the classic hexagonal Wigner monolayer is realized, with x = 0 (structure I). For large asymmetry (smaller A), the monolayer stability is, expectedly, augmented. Regions where the system either remains in phase I or partly populates layer 2, are separated by a curve in the (η, A)-plane, denoted by η c (A) [or conversely A c (η)] shown in Figure 1. For A 0.4085, this curve separates phase I from a family of phases that are denoted by I x : the latter originates from the monolayer (i.e., a hexagonal lattice α on plate 1 with spacing a) by picking a fraction x of particles in a hexagonal arrangement to relocate them on plate 2 (with thus spacing b > a). An illustration is provided in Fig.1. For these particular structures, a complete analytical analysis can be achieved and simple geometric considerations imply that only a discrete set of x-values is compatible with this geometric constraint: x = 1/(j 2 + jk + k 2 ), with non-negative integers j and k such that j + k > 1: x ∈ {1/3, 1/4, 1/7, 1/9, · · · }. As x → 0, these values become essentially dense, so that we can consider x in this regime as a quasi-continuous variable. A sufficient condition for instability of phase I is that it becomes energetically favorable to extract one particle from the monolayer, keeping all others in position. This leads to an upper bound for η c (A), shown by the thick curve in Fig. 1, quite close to the boundary obtained by the EA algorithm. On this curve, the increase in the potential energy of a tagged particle shifted from plate 1 to plate 2 is balanced by a decrease in the particle's interaction energy (correlation term). For A 0.4085, phases competing with phase I originate from a different mechanism. This family of phases, denoted by V x , is made up of two triangular (hexagonal) structures on the two plates (lattices α and β) with some shift, see Fig.1 for the special cases x = 1/2 (leading to structure V), and also x = 1/4. Note that when the rescaled distance η → ∞, one expects structure V x with x = x neutr = A/ (1 + A). Considering x as a continuous variable, one can calculate analytically the location of the transition line η c (A), depicted in Figure 1, along with its EA counterpart. We present the essence of the calculation, which sheds light on the critical behavior. For a given A-value, the energy difference between structures I and I x can be written in a small-x expansion as where f (η) is a closed expression in η, which also depends on A, and λ 1 [40]. The order parameter x vanishes at the transition, i.e., at a point where η = η c , fixed by the condition f (η c ) = 0. This relation yields the continuous curve η c (A) in Fig. 1 and can be viewed as a locus of critical points. Besides, expanding f for η > η c up to linear order in η gives access to the critical index. Together with the extremum condition of (1) with respect to x, this leads to the prediction x ∝ (η − η c ) β with β = 2/3 [40]. This exponent differs from its Ginzburg-Landau theory counterparts, based on an energy expansion that is analytic in the order parameter. Here, the long-range nature of Coulomb interaction breaks analyticity. Numerical results are fully compatible with β = 2/3, not only for A 1 where it is admissible to neglect lattice distortions (see Fig. 2), but for all A values, along the full curve η c (A). Transitions I→I x and I→V x therefore share the same non-standard exponent β = 2/3. On the other hand, the other transitions such as II → III and III → IV are standard: there, the analytical treatment is rigorous, and leads to mean-field continuous transitions, with β = 1/2 [40]. Thus, fixing A, it is remarkable that a sequence of transitions with distinct critical indices takes place when increasing η. Fig. 1 provides the stability domain of structures I x . Moving away from the η c (A)-curve by increasing η, a cascade of I x -phases emerges, each of them corresponding to an x-value specified above and associated to a plateau on the left hand side of Fig. 2. A snapshot of structure I 1/4 is shown in figure 1. Noteworthy is the honeycomb lattice (H phase) on plate 1, structure I 1/3 ; the corresponding plateau in Fig. 2 for 0.07 < η < 0.17 is of significant extension. Comparing the boundaries of the H-phase, evaluated via the numerical and the analytical methods shows an excellent agreement for 0 < η 0.30, confirming thereby the absence of distortions of the optimal structure within that region. For η > 0.45, the analytical approach still establishes H as the most stable phase (see Fig. 1), while EA and Monte Carlo identify novel intricate snub-square or pentagonal structures, see below.
Phases V x extend over a significant area in the (A, η)plane of Figure 1. The agreement numerical/analytical is fair, with discrepancies arising from the emergence of complex, distorted structures. At large distances though, where structures α and β are undistorted, some exact statements can be put forward [40]: (i) for A < 1 the plates (plus ions in contact) remain charged at any finite distance; only as η → ∞ and/or A → 1, does x approach the electro-neutral value x neutr ; (ii) while the always attractive inter-plate pressure is short-range (exponential) for the symmetric case, it becomes long-range whenever A = 1, decaying like 1/η 2 [33]. A detailed analysis of structures V x reveals that at intermediate distances, they can accommodate significant distortions, leading to new structures coined DV x in Fig. 1 [34].
We now focus on the vicinity of the symmetric line A = 1, where the structures that prevail for A = 1 do exist in some parameter range (A larger than 0.9). These are the (x = 1/2)-phases II, III and IV. As the symmetric structures II to IV are undistorted, both the analytical and the numerical approaches predict regions of stability that are in perfect agreement (see the vertical lines in the upper part of Fig. 1). While for phases III and IV, no generalizations to x-values different from 1/2 have been identified, phases II x emerge in a small pocket of the (A, η)-plane. The numerical EA approach indicates a continuous II x → II transition. The fact that several undistorted II x structures can be classified via well-defined sequences of alternating rows of particles in the two layers (see Fig. 1 where the II x=2/5 arrangement is depicted), opens the possibility of an analytical, exact energy calculation [40]. The numerical and analytical routes provide fully consistent results for the stability of these phases.
The upper part of the phase diagram is the locus of a rather unexpected phenomenon of charge reversal. While for A = 1 each plate plus ions in contact is electro-neutral at all distances (x = 1/2), the majority of identified states is characterized by undercharging: plate 2 (with density 0 < σ 2 < σ 1 ) plus the (negative) ions in contact carry a net positive charge. Thus plate 2 attracts less mobile charges than required for neutrality, so that x < x neutr . This is somewhat expected but is no longer the case for A 1, where overcharging takes place: the most weakly charged plate attracts more ions than necessary for neutrality so that x > x neutr , see Fig. 2 [35].
Finally, we report more exotic phases, starting with the snub type. The regular snub square structure S 1 shown in Fig. 1  square lattice [36], while particles in layer 2 arrange in a square lattice (with slight deformations as η grows). Since particles in layer 1 have five nearest neighbors, the S 1 phase can be quantified via the five-fold order parameter Ψ (1) 5 together with Ψ (2) 4 [40]. Interestingly, the undistorted snub square lattice is an Archimedean tiling [36,37]. Such a geometry is amenable to an analytical treatment [40]. Surprisingly, another snub square structure, denoted S 2 , can be identified with x = 1/3 as well. In contrast to the S 1 phase, it shows stronger deformations, which decrease Ψ (2) 5 . Both structures occupy relatively small regions in Fig. 1, where pentagonal (P) phases are also reported [38].
We have considered a charged bilayer system governed by two parameters only: the charge asymmetry A between the parallel plates, and the slab width η, more prone to experimental tuning. The competition of slit confinement with Coulomb interactions leads to a plethora of ordered bilayers with phase transitions pertaining to different universality classes. In light of the simplicity of the model, the complexity and variability of emerging phases is striking. Patterns emerge as a tradeoff between the commensurability of structures, and incomplete charge neutralization (the latter effect being quantified by x − x neutr ). Fig. 1 summarizes our main findings. Besides possible experimental confirmation in quantum wells [17], semiconductor bilayers [16,18], bilayer graphene [19], or ionic plasmas [39], other relevant perspectives deal with the inclusion of an ionic hard core, which would frustrate several of the arrangements put forward, the study of disordered or patterned substrates, as well as the analysis of dynamical processes and elementary excitations.
with L = 1 and L = 2 would be more relevant here, but would delineate the same region of existence. It should be kept in mind however that structures Vx are expected at large distance η, and are not of particular interest as such.
[33] Likewise, the approach to the large distance energy turns from exponential (when A = 1) to algebraic (when A = 1).
[34] Such structures can be defined by their symmetry group [40], by their large 6-fold bond orientational order parameters for both layers, their significant 4 or 8-fold order for layer 1, and 5 or 10-fold order for layer 2 due to the distortions of the triangular lattices.
[35] The curve that separates undercharged (x < xneutr) from overcharged (x > xneutr) states coincides with the lower boundary of structures II, III and IV in Fig. 1 Fig. 1 does not clearly distinguish between S2 and P patterns. Making use of parameters Ψ (2) n clarifies the boundaries [32].