Table Of Content

Enhancement of Wigner crystallization in quasi low-dimensional solids. G. Rastelli1, P. Quémerais2, and S. Fratini2 1 Istituto Nazionale di Fisica della Materia and Dipartimento di Fisica Università dell’Aquila, via Vetoio, I-67010 Coppito-L’Aquila, Italy and 2 Laboratoire d’Etudes des Propriétés Electroniques des Solides, CNRS BP 166 - 25, Avenue des Martyrs, F-38042 Grenoble Cedex 9, France (Dated: January 25, 2006) 6 0 Thecrystallizationofelectronsinquasilow-dimensionalsolidsisstudiedinamodelwhichretains 0 thefullthree-dimensionalnatureoftheCoulombinteractions. Weshowthatrestrictingtheelectron 2 motion to layers (or chains) gives rise to a rich sequence of structural transitions upon varying the particle density. In addition, the concurrence of low-dimensional electron motion and isotropic n a Coulomb interactions leads to a sizeable stabilization of the Wigner crystal, which could be one of J themechanisms at theorigin of thechargeordered phases frequentlyobserved in such compounds. 5 2 I. INTRODUCTION ferent chemical units are so anisotropic that the carrier ] motion is effectively restricted to two-dimensional (2D) l atomic layers, or one-dimensional (1D) chains. Yet, the e Despite being a well established concept in the - physicsofinteractingelectrons,directevidenceofWigner Coulombforcesretaintheir three-dimensionalcharacter, r t crystallization1hasbeenreportedunambiguouslyinonly being long-ranged and isotropic. In such systems, inter- s . a limited number of systems, namely electrons at the layer(inter-chain)interactionscannotbeneglected,lead- at surface of liquid helium,2,3,4 and in semiconductor het- ing eventually to a full three-dimensionalorderingof the m erostructures of extreme purity.5 In both cases, a two charges.11,15 This suggests why quasi low-dimensional - dimensional electron gas (2DEG) is realized at an inter- solids are a particularly favorable ground for the obser- d face between two media, for which the jellium model of vation of Wigner crystallization: the electron-electron n thehomogeneouselectrongasconstitutesagoodapprox- interactions have the same behavior as in bulk three- o imation. dimensional systems, but the kinetic part is strongly re- c duced by the effective lowering of dimensionality. Re- [ Alternatively,thechargeorderingphenomenaobserved at low temperatures in a number of solids have been ei- minding that a Wigner crystal arises from the competi- 1 tion between potential and kinetic energy, this results in ther interpreted as some form of Wigner crystallization, v a sizeable stabilization of the crystal as compared with or ascribed to the presence of long-ranged Coulomb in- 8 the usual 3D case.42 7 teractions. These include the one-dimensional organic 5 salts TTF-TCNQ,6 (TMTTF)2X,7,8 (DI-DCNQI)2Ag,9 Asimilarconclusionisreachedbyobservingthat,even 1 the ladder cuprate compounds Sr Cu O 10,11 and comparedto purely low-dimensionalsystems such as the 14 24 41 0 chain compounds Na CuO ,12 as well as the lay- 2DEG mentioned above, the Wigner crystal phase could 1+x 2 6 eredsuperconductingcuprates14,15 andpossiblythetwo- be stabilized in quasi low-dimensional solids due to the 0 dimensional BEDT-TTF organic salts.13 For such sys- presence of additional interlayer interactions. This topic / t tems, the jellium model is a priori a rather crude mod- has been analyzed in the literature in the framework of a m elization,andtheconceptofWignercrystallizationmust bilayer quantum wells, i.e. constituted of two coupled begeneralizedtoaccountforothercompetingeffectssuch 2D electron systems, where it has been shown that, de- - d as the periodic potential of the underlying lattice, chem- pending on the strength of the interlayer forces, the or- n icalimpurities, structuraldefects, magnetic interactions, dering pattern can differ from the hexagonal structure o etc. In narrow band solids, for instance, the interplay expected in a single layer.19,20 More importantly, it was c withthehostlatticeofionscanstronglyaffectthecharge found21,22,23 that at interlayer separations comparable : v orderingpatternespeciallyathighlycommensurateband withthemeaninterparticledistance,themeltingdensity i fillings.6,16 Nevertheless, when the radius of localization israisedbyafactorof3withrespecttothepure2Dcase, X oftheparticlesislargerthanthetypicalion-iondistance, which makes a factor as large as 102 when appropriately ar the host lattice can be replaced to a good accuracy by scaled to the 3D situation. an effective continuous medium, restoring de facto the In this work, we model quasi two-dimensional (one- validity of the jellium model.17,18 dimensional) systems as periodic arrays of conducting Setting aside the important problem of the commen- layers(wires)embeddedina three-dimensionalbulk ma- surability with the host lattice, and neglecting disorder terial, where the electrons interact through isotropic and other effects that can certainly play a role in the long-range Coulomb forces. We show that, upon vary- compounds under study, we come to the following ob- ing the particle density or the interlayer (interwire) sep- servation: a common feature shared by the experimen- aration, the Wigner crystal undergoes several structural tal systems listed above is that they are all quasi low- transitions in order to minimize its energy compatibly dimensional solids, i.e. they are bulk three-dimensional with the given geometrical constraints. We then give a (3D)compoundswherethetransferintegralsbetweendif- semi-quantitative estimate of the melting density for the 2 different structures previously identified, based on the powers of 1/r1/2.24,25 The leading term, proportional to s Lindemann criterion, which confirms the stabilization of 1/r ,correspondsto the MadelungenergyE ofEq.(2). s M the crystallized phase expected from general grounds. In free space, it attains its minimum value E = BCC The paper is organized as follows: In Section II, we 0.89593/r (inatomicunits)foraBodyCenteredCubic s introduce a model for the crystallization of electrons in −(BCC) Wigner crystal.26 The second term in the expan- ananisotropicenvironmentandthe method forcalculat- sion, proportional to 1/r3/2, is the zero point energy of s ing the crystal energy in the harmonic approximation, the particle fluctuations in the harmonic approximation, which includes the classical Madelung energy and the which also depends on the selected crystal structure. It zero-pointvibrationalenergyofthecollectiveexcitations. is negligible for r , and remains smaller than the s This is applied to the case of quasi two-dimensionalsys- Madelung term by→typ∞ically an order of magnitude at tems, for which the structural/melting phase diagram is r 100. Nonetheless, it can play an important role in s determined. The validity of the present approximation det∼ermining the relative stability of the different crystal scheme is checked at the end of Section II by analyz- structures,especiallywhenapproachingthemeltingden- ing a system of two coupled layers,for which our results sity. Higher orders in the energy expansion24,25 include compare positively with the numerical results available anharmonic (1/r p with p 2) and exchange terms of s in the literature. An analogous discussion for quasi one- the form e c√rs, which we ≥shall neglect in the following − dimensional systems is reported in Section III, by treat- discussion. ingexplicitelythecasewheretheconductingchainsform Uptoquadraticorderinthe displacements,ourmodel a square array. The main results are summarized in sec- Hamiltonian reads: tion IV. p2 e2 H =NE + i + (~u ~u )Iˆ (~u ~u ) (4) M i j ij i j 2m 4 − − i i,j=i II. WIGNER CRYSTALLIZATION IN LAYERED X X6 SOLIDS whereIˆ isa2 2matrixcharacterizingthedipole-dipole ij × interactions, given by (α,β =x,y): A. Model and approximations 3R~ R~ δ Iˆ = ij,α ij,β αβ (5) Letusconsiderasystemofelectrons(orholes)ofden- ij αβ Rij 5 − Rij 3 (cid:16) (cid:17) | | | | sity n = (4πr3/3) 1 in a strongly anisotropic environment, such thsat t−he particle motion is constrained to withR~ij =R~i R~j. ThemostgeneralelementaryBravais − equallyspacedatomiclayers(atdistanced),butremains lattice compatible with a givenlayeredstructure is iden- isotropic within the layers. To ensure charge neutral- tifiedbyacoupleofbasisvectorsdescribingthe ordering ity, we assume a uniform 3D compensating background within the planes, A~1 = (a1,0,0), A~2 = (a2x,a2y,0), of opposite charge. The hamiltonian for N crystallized and a third vector A~ = (a ,a ,d) which sets the rel- 3 3x 3y particles in a volume V is given by: ative shift (a ,a ) between two equivalent 2D-lattices 3x 3y onneighboringplanes. Otherstructures,withmorethan N p2 one particle per unit cell, are possible in principle, but H =NE + i +V (1) M d will not be considered here. 2m Xi=1 Due to the additional lengthscale d introduced by the layeredconstraint,the crystalenergyisnolongerafunc- The first term tionofr alone. Itsdependenceonthelatticegeometryis s e2 1 d~r bestexpressedbyintroducingadimensionlessparameter EM = n (2) γ,whichmeasuresoftherelativeimportanceofinterlayer 2 " i Ri − ZV r # and intralayer interactions. It is defined as the ratio be- X tween the mean interparticle distance in the planes and is the Madelung energy of the given lattice structure (in the interlayer separation, namely γ = √πr /d. Here the thermodynamic limit, N,V , boundary effects s,2D are negligible and all particles be→com∞e equivalent). The rs,2D defines the 2D density parameter in the individual layers, related to the bulk r by r2 = 4r3/3d. The secondtermis the two-dimensional kinetic energyofthe s s,2D s firsttwotermsofthelow-densityexpansion,correspond- localized particles and the last term accounts for the in- ing respectively to the Madelung energy and the zero- teractionsduetotheplanar displacements~u =(u ,u ) i xi yi point fluctuation energy in the quadratic model (4) can of the electrons around their equilibrium positions: be written in compact form as: e2 1 1 A(γ) B(γ) V = (3) E = + . (6) d 2  R~ +~u R~ ~u − R~ R~  rs rs3/2 i=j i i j j i j X6 − − − (cid:12) (cid:12) (cid:12) (cid:12) It should be noted that an effective mass m∗ =m and Expanding the(cid:12) last term for sm(cid:12)all d(cid:12)isplacem(cid:12)ents re- a dielectric constant κ = 1 can be straightforw6ardly in- (cid:12) (cid:12) (cid:12) (cid:12) 6 sults in a series expansion for the energy Eq. (1) in cluded in the model through a redefinition of the Bohr 3 radius a a = a κ(m/m ), unit energy me4/~2 m e4/κ2B~2→, an∗Bd densBity para∗meter r r (m /m)/→κ. H -0.8952 ∗ s → s ∗ CR Hereafter, energies and lengths will therefore be ex- S -0.8954 pressed in terms of these effective units, characterizing R CR the hostmedium. Amuchmorecomplexsituationarises Rh in systems with a frequency-dependent dielectric screen- -0.8956 ing, leading to the formation of polarons, for which the reader is referred to Refs.15,27. A -0.8958 -0.8960 B. Minimization of the Madelung energy H CR S R CR Rh H CR Rh R -0.8962 Following the hierarchy of the series expansion intro- ducedabove,westartbysearchingforthelayeredconfig- 0 1 2 3 4 5 6 uration which minimizes the electrostatic repulsion be- γ tween the particles, which is appropriate in the limit FIG. 1: Madelung coefficient A in atomic units for differ- of large r . The calculation is performed by standard s ent crystal structures constrained to a layered environment, Ewaldsummationtechniques,whichsplittheslowlycon- as a function of the anisotropy ratio γ. The different curves vergent series in Eq. (2) into two exponentially converg- correspondtodifferentplanarconfigurations: hexagonal(H), ing sums.28 Given the interlayer separation d and the square (S), centered rectangular (CR), rectangular (R) and bulk density n (or, alternatively, given the pair of di- rhombic (Rh). The interlayer orderings in the simplest cases mensionless parameters γ and rs) we are left with 4 free at low γ are sketched below the curves (for the R and CR minimization parameters: 2 for the inplane structure, 2 structures, the stacking varies as indicated by the double ar- for the interlayer ordering. rows). The resulting three-dimensional Wigner crystal re- The result of the minimization for the Madelung co- duces to a perfect BCC at the points marked by filled dots, efficient A in the range 0 < γ < 6 is illustrated in Fig. whose energy is indicated by thehorizontal arrow. 1. Twodistinct regimescanbe identified. Inthe limit of large separations (γ . 1), the coupling between the lay- ersisweak,andtheresultingplanarpatternishexagonal, (S) for γ . 2, as expected for large interlayer separa- with a staggeredinterlayerordering,i.e. the particles on tions, where the relative ordering is fully determined by the neighboring layers falling on top of the centers of the coupling between two adjacent planes, and indeed the triangles. The sharp rise of the Madelung constant coincides with what is found in bilayer systems19 (see in this regime is due to the fact that the compensat- Section IIE below). At larger values of γ, the interac- ing background is distributed homogeneously in three- tions beyond the nearest planes become relevant, which dimensional space, which penalizes strongly anisotropic makesthesimplestaggeredorderingunfavorable. Forin- charge distributions.43 stance, a staggered/non-staggeredtransition takes place Uponreducingthe interlayerseparationsothatγ &1, within the rectangular phase at γ =2.46, corresponding the increasing interlayer interactions make the hexago- to a relative sliding of the planar structures on adjacent nal pattern energetically unfavorable. Above γ = 1.15, planes in the direction of the long bonds (indicated by a more isotropic ordering of the charges is stabilized, the double arrow in Fig. 1). which presents a centered rectangular (CR) structure in Remarkably, each of the phases identified above con- the planes. Further increasing γ leads to a sequence of tainsaspecialpointγ∗ wheretheidealBCCstructure— structureswhoseplanarpatternsarerespectivelysquared which has the lowest possible Madelung energy in three (S, in the interval 1.32 < γ < 2.13), rectangular (R, dimensions— is itself compatible with the layered con- 2.13 < γ < 2.84), centered rectangular (CR, 2.86 < straint. The different planar configurations identified γ < 4.31), a generic rhombic, or oblique phase (Rh, above then correspond to the different ways of cutting 4.31<γ <4.45),then hexagonalagain,and soon. Such a BCC by anarrayofequally spacedlayers. Suchpoints phases are all connected by continuous structural tran- are easily calculated by setting the distance d=2π/K, sitions, with the exception of the hexagonal structure, with K any reciprocal lattice vector, and correspond| t|o which is attained through a discontinuous change of the γ∗ = 21/4, 2, 21/433/4,21/453/4, 2 33/4, etc...Similarly, crystalparameters. Notethatintheverynarrowinterval the higher relative minima visible in Fig. 1 correspond 2.84<γ <2.86,agenericstructurewithrhombicplanar to different orientations of the same three-dimensional symmetry is stabilized, which allows to evolve continu- Face Centered Cubic (FCC) ordering. ously from the rectangular to the centered rectangular Away from such special points, the overall charge dis- patterns (not shown). The sequence ofstructuraltransi- tribution remains very isotropic in all the region γ & 1, tions goes on at larger values of γ. as testified by the extremely small deviations of the TheinterlayerorderingisshownatthebottomofFig.1. Madelung energy from the ideal case, ∆E . 10 4/r . M − s It is staggeredfor the first three patterns (H), (CR) and Such small energy variations, however, refer to the op- 4 timal structures obtained at different values of γ, which a) does not mean that the electrostatic repulsion between H 1.00 CR the carriers is irrelevant in the determination of the 0.98 S charge ordering patterns in real systems: in a given R 0.96 compound, where both the interlayer distance and the Rh. 0.94 density are fixed, one should rather compare the energies of two competing phases at fixed γ. For exam- B 0.92 ple, enforcing a hexagonal symmetry at γ = 2, where 0.90 the optimal structure is squared, would cost an energy 0.88 ∆E 0.015/r 200K at r = 20, which is compa- M s s ∼ ∼ rable with the typical charge ordering energy scales in 0.86 H S R C R H Rh H solids. Yet, since the Madelung energy is determined by 0.84 the interactions with a large number of (distant) neigh- bors, the structures found here are expected to be rel- atively soft against local deformations. The situation is b) different regarding global symmetry changes, as can re- 2.8 sultfromthe inclusionofaperiodic potentialofcompet- ing symmetry,whichcouldstronglymodify the sequence 2.6 and order of the structural transitions, possibly favoring the appearance of alternative phases.29,30 1 M- 2.4 C. Zero point fluctuation energy 2.2 The next term in the series expansion of the ground 2.0 state energy Eq. (6) corresponds to the quantum zero 0 1 2 3 4 5 point fluctuations of the particles around their equilib- γ rium positions, in the harmonic approximation. It is negligibleatlargers (lowdensity), butitbecomesquan- FIG. 2: a) Zero point vibrational term B for the different titatively important at lower r , where it can slightly structures identified in Fig. 1, within their ranges of me- s modifythesequenceofphasesidentifiedinthepreceding chanicalstability. Thesequenceofstructureswiththelowest Section. Upon further reducing r , this term eventually vibrational energy is indicated at the bottom. The arrow s drives the quantum melting of the crystal, that will be indicates the value (2/3)B(3D) = 0.887, where B(3D) is the vibrational energy of a BCC crystal in vacuum; b) Inverse analyzed in the next Section. Thecalculationofthefluctuationtermproceedsasfol- momentM−1 of theDOS,which isproportional tothemean electronic fluctuation u2 Eq.(12) for the same structures. lows. The harmonic model Eq.(4) is diagonalized by in- h i The arrow indicates thevalue for a BCC in free space. troducing the normal modes q s,~k 1 ~ui = √N εˆs,~keı˙~k·R~iqs,~k (7) (D is the number of branches, corresponding to the di- Xs,~k mensionality of the electron motion) as well its dimen- sionless moments: where εˆ are the two-dimensional polarization vectors s,~k n (the electrons oscillate within the planes) and the vec- ω M = dω ρ(ω) (10) tor ~k runs through the Brillouin zone of the three- n ω Z (cid:18) P(cid:19) dimensional reciprocal lattice. This yields two branches s=1,2ofcollectivemodeswitheigenfrequenciesω ,so with ω2 = 3e2/(mr3) the usual 3D plasma frequency. s,~k P s thatthe vibrationalenergyperparticlecanbe expressed With these definitions, the vibrational energy in Eq.(6) as: is seen to be directly proportionalto the first momentof E = 1 ~ωs,~k (8) the DOS, with V N 2 D√3 Xs,~k B(γ)= M1(γ) (11) 2 It is useful to introduce the normalized density of states (DOS)ofthecollectivemodes,thatwewritehereingen- The usual 3D case in vacuum is recovered by restoring eral as: the out-of-plane oscillations in Eq. (4), and by setting D = 3 in Eq. (11). For example, for the BCC structure D ρ(ω)= 1 δ(ω ω ) (9) wefindM(3D) =0.511,whichyieldsthewellknownvalue DN − s,~k 1 Xs=1~kXBZ B(3D) =1.33.25,31 ∈ 5 The analysis of the frequency spectrum shows that each given structure has a limited interval of mechani- rs CR R cal stability: for certain geometries, the dynamical ma- trix acquires negative eigenvalues around some critical 10000 CR (Rh) wavevector k , corresponding to purely imaginary col- c lective frequencies which drive the crystal unstable (this 1000 (Rh) Rh phenomenon also exists in free space, where FCC and the simple cubic structure are known to be intrinsically unstable). For example, in the intervalof γ under study, 100 H S R CR H a structure with hexagonal symmetry is only stable for (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) H γ <1.32, 3.5<γ <4.95 and 5.05<γ <5.8. (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) 10 (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) We have calculated the fluctuation term B(γ) for the (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) Liquid different symmetric structures (H, R, S, CR) identified in the previous section, within their respective intervals 1 of mechanical stability, as well as for the rhombic phase 1 2 3 4 5 at 2.84 < γ < 2.86 and γ > 4.31, which is shown in γ Fig. 2.a. As for the Madelung energy, two essentially different regimes can be identified. For γ . 1, the elec- FIG. 3: Structural phase diagram of the Wigner crystal in a layered environment, based on the total energy (6), as a tron motion is mostly determined by the Coulomb in- function of the anisotropy ratio γ and the bulk density pa- teractions within the layers (interlayer forces are neg- rameter r . The labels are the same as in previous figures. ligible) and the collective modes of the pure 2D case s The solid (dashed) lines are for structural transitions where are recovered. If normalized by an appropriate “two- thecrystalparametersevolvediscontinuously(continuously). dimensional plasma frequency” ω2 = e2/mr3 , the 2D s,2D Theboldlinesindicatemechanicalinstabilities, accompanied first moment in the hexagonal phase tends to the con- by a discontinuity of the crystal energy. The melting line is stantvalueM1,2D = 0.814.32 Goingbacktothe present determined by solving Eq. (13). For the hatched region, see three-dimensionalunits,however,wherethemomentsare text. normalized as in Eq.(10), the fluctuation term diverges at large separations as B(γ) π1/431/2M /2γ1/2. In 1,2D theregimeγ &1,ontheother≃hand,thefluctuationterm to slightly modify the range of mechanical stability of flattensaroundavaluewhichroughlycorrespondsto2/3 each phase. ofthefluctuationinfreespace,indicatedbythe arrowin Another fundamental property of the system, which Fig. 2.a. This followsfromthe factthat onlythe oscilla- gives valuable informations on the collective vibrations tions along 2 of the 3 space directions are allowed,as we of the particles, is the mean electronic fluctuation u2 . can see explicitely from Eq. (11). In the harmonic approximation, this quantity is propor- (cid:10) (cid:11) Thestructuralphasediagramresultingfromtheanaly- tionaltotheinversemomentoftheDOSofthecollective sisofthe totalenergy(6), including thevibrationalterm modes, defined in Eq. (10): (11), and taking into account the ranges of mechanical 1 1 DM sfitrasbtiolibtyseorfvathtieodniffisetrhenatt,pahpaasrets,frisomreptohreteddisianpFpeiga.ra3n.cTehoef u2 = N 2ωk,s = 2√−31rs3/2 (12) k,s the CR phase from certain intervals, which is penalized (cid:10) (cid:11) X by its higher vibrational energy than the H phase, the where again we keep track of the explicit dependence locus of the structural transitions does not change much on the dimensionality D. As can be seen in Fig. 2.b, with rs. The sequence of phases identified in Fig. 1, it increases as each phase approaches the boundaries of based on the analysis of the Madelung energy, is recov- its stability range. This is because the mechanical in- ered at extremely large values of rs. On the other hand, stabilities are approached via a softening of a branch of the vibrationalterm affects the structural transitions al- collective modes, causing an increase of the DOS at low ready at rs . 1000. This is due to the fact that, even frequency and, through Eq. (10), of the inverse moment thoughtheelectrostatictermA/rs isstilllargerthanthe M 1. Alocalincreasealsooccursatthepointswherethe zero-point fluctuation energy B/r3/2, the latter under- sta−ggered interlayer ordering is lost (see e.g. the maxi- s goes much larger relative variations among the different mum at γ =2.46 within the R phase in Fig. 2.b). phases. Belowr 100,thephasediagramisentirelyde- From analogous arguments, it follows from Eqs. (10) s ∼ termined by the minimization of the vibrational energy and(11)thatthe vibrationalenergygenerallyattainsits (see Fig. 2.a). As was stated above,however,the overall minimum value close to mechanical instabilities. Within shapeofthephasediagramdoesnotdependmuchonr , the present approximate framework, this can cause the s the transitions being essentially determined by the pa- total energy to jump discontinuously at the instability rameter γ. Let us also remark that the vibrational term point when the next stable phase is attained, which cor- is much less influenced than the Madelung term by the responds to the bold lines is Fig. 3. For example, the specific interlayer arrangements, whose effect (if any) is hexagonal lattice becomes unstable at γ > 1.33 and, for 6 r .100, the transition to the square phase is accompa- isdirectlyreflectedinFig. 2.binareducedvalueofM s 1 nied by a small jump in energy. Such discontinuities can ascomparedtothecorrespondingvalueinfreespace,an−d in principle be avoided by allowing for Bravais lattices shouldnotbeconfusedwiththetrivialdimensionalfactor with more than one electron per unit cell (the resulting D,thatwastakenoutexplicitelyfromEq. (12). Itisdue internal structure could then be assimilated to some lo- tothefactthat,assoonasthecubicsymmetryislost,the cal tendency to electron pairing33,34). Note also that it restoringforcesinducedbythedipole-dipoleinteractions is precisely close to mechanical instabilities, where u2 Eq. (5) are not equivalent in the three space directions, h i is largest, that the neglected anharmonic corrections are sothatthe electronfluctuationbecomesanisotropic(the expected to be most important. Their consequences on observed shrinking of the planar spread would occur at the structuralphase diagrampresentedhere deserve fur- the expense of increasing the out-of-plane fluctuations, ther theoretical study. which are anyhow suppressed in the model). To give an example, taking an averagevalue M 2.4 and the as- 1 pectratio =√π forthesquareplan−ar≃orderingyieldsa C D. Melting of the crystallized state criticalvaluerc 4.9√d. Comparableresults(within s,2D ≃ few percent) are found for the other structures. In this section, we analyze the melting of the crystal- lizedstatebymakinguseoftheLindemanncriterion,ac- cordingtowhichatransitiontoaliquidphasetakesplace 50 whenthespread u2 attainssomegivenfractionδ ofthe (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)H(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)C(cid:0)(cid:1)(cid:0)(cid:1)R(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) 45 (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) R nearest-neighbor distance a . We take δ = 0.28 from (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) Ref.35, which is a(cid:10)pp(cid:11)ropriaten.fno.r the quantum melting of 40 (cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1) S (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) both 2D and 3D Wigner crystals. Solving the equation 35 (cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1) H hu2i/an.n. =δ in terms of the density parameter rs,2D 30 (cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1) in the planes leads to: D (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) p s,2 25 (cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1) r (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) M (γ) 20 (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) rsc,2D = 2δ2−12(γ)d1/2 (13) 15 (cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1) C 10 (cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1) Liquid where = a /r is an aspect ratio relating the n.n. s,2D 5 C nearest-neighbordistancetothedensityparameterinthe 0 planes, and the implicit condition γ =√πrc /d holds. s,2D 0 5 10 15 20 25 30 35 40 45 50 Note that for structures with rectangular symmetry, the d Lindemanncriterionmustbemodifiedtoaccountforthe existence of two nonequivalent near-neighbor distances. FIG.4: PhasediagramoftheWignercrystalinalayeredenvi- Here we use a simple generalization which consists in ronment,asafunctionoftheinterlayerdistanced(inunitsof ∗ replacingan.n. withthe averageofthe twoshortestnear- the effective Bohr radius aB). The discontinuity close to the neighbor distances, and which reduces to the ordinary three-phasecriticalpoint(S,H,liquid)isduetothedifferent aspectratios ofthetwocompetingstructures. Thehatched criterion for the square and hexagonal structures. A C areaisaregionpossiblycharacterizedbyananisotropicliquid check of the validity of such generalized Lindemann cri- behavior (see text). The shaded area corresponds to γ > 6 terionwillbe giveninSectionIIE,bydirectcomparison and has not been studied. with independent theoretical results on bilayer systems. ThemeltingcurvededucedfromEq. (13)forthediffer- entstructuresconsideredhereisillustratedinFigs.3and In the opposite limit of large separations (γ 1), 4. The most important result is that the crystal melt- where interlayerforces become negligible, we recov≪erthe ing can be pushed to higher densities by reducing the usual critical value rc 40 for the 2D hexagonal interlayer spacing, which can already be inferred by ne- Wigner crystal. Notes,2tDhat≃the actual critical value at glecting the weak γ-dependence of the coefficients and finite γ always lies below this asymptotic estimate, con- C M 1 of Eq. (13) in the region γ & 1. The main reason firmingthattheinclusionofinterlayerinteractionscauses to−this is that for γ &1 the electron spreadis essentially a stabilizationofthe crystalphase comparedto the pure governed by three-dimensional Coulomb interactions, as two-dimensionalcase,aswas arguedin the introduction. we can see from the explicit dependence of Eq. (12) on Afewcommentsonthelimitsofvalidityofthepresent thebulkrs,whiletheelectronmotionistwo-dimensional, model are in order. First, the enhancement of Wigner sothattheappropriatenearest-neighbordistanceforthe crystallizationpredictedbyEq. (13)cannotextendindef- Lindemann ratio is proportional to the planar density initely: themeltinglineshouldeventuallysaturateatlow parameter rs,2D =(2/√3d)rs3/2.36 separations when isotropic electron motion and three- In addition, for each given spacing d, the geometrical dimensional screening are restored by interlayer tunnel- confinementleadstoafurtherstabilizationofthe crystal ing processes.23 On the other hand, as was stated in the through a reduction of the spread u2 itself. This effect introduction, replacing the host lattice of ions by an ef- h i 7 r s,2D 10 15 20 25 30 35 40 45 50 70 H 0.32 60 R (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) S 0.28 (cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1) 50 S (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) H (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) 0.24 (cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)R(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1) 40 (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) D (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) 2 0.20 (cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1) rs, 30 R S Rh H (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) H 0.16 (cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)C(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)R(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1) H 20 S Liquid 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Lindemann 10 γ QMC DFT 0 FIG. 5: Lindemann ratio u2 /an.n. as a function of γ 0 10 20 30 40 50 60 70 80 h i (r , upper abscissa) at a given interlayer spacing d = 20. s,2D p d The continuous line is the average Lindemann ratio, the dashed line represents the Lindemann ratio along the direc- FIG. 6: Phase diagram for the symmetric electron bilayer, tion of the closest near-neighbor (see text). The horizontal in terms of the two dimensional density parameter r , as s,2D line sets the critical value for melting. Upon increasing the afunction of theinterlayer spacing d. Thelengthsare scaled electron density, the transition from the crystal to the liquid to atomic units. H,S,R denote respectively the Hexagonal, couldoccurthroughan intermediateanisotropic liquidphase Square , Rectangular phase. Note that the QMC simulation (hatched region, see also Fig. 4). ofRef.23 wasrestrictedtostudyonlytwophases(H,S),while an additional rhombic phase (Rh) could be stabilized in the DFT calculations of Ref.21. fective jellium is allowed providedthat the spread of the electron wavefunction is larger than the ion-ion distance the analysis of the Madelung energy showedthat several a . Fromeq. (12),thecondition u2 &a givesr &4 0 h i 0 s structural phase transitions occur as the distance d be- (r &6)foratypicalvalueofa =3˚A,assumingκ=1 s,2D p0 tweenthe twoplanes is variedwhile keeping the electron andm =m. Belowthisvalue,thediscretenatureofthe ∗ density fixed. At short distances the planar ordering is host lattice should be included, which can further stabi- rectangular, and collapses to the usual hexagonal phase lize the crystallized state, as pointed out in Refs.17,18,38. in the formal limit d 0. This phase evolves continu- Before concluding this section, let us remark that, for ously into a staggered→square structure, when d is of the anisotropic planar orderingssuch as the rectangular and order of the interparticle spacing, which is clearly remi- thecenteredrectangularstructure,twoindependentLin- niscent of the BCC structure observed in 3D space (cf. demannratioscouldinprinciplebedefined(oneforeach the discussion in Section II B). Upon further increasing nonequivalent near-neighbor direction) rather than the d,thelatticeprogressivelydeformsintoarhombicphase, singleaveragecriterionusedsofar. Itwouldthenappear to attain the hexagonal staggered phase expected in the that the melting along the short bonds is much easier limit of independent layers. Including the zero-point en- than along the long bonds, due to the closer overlap be- ergy of the collective excitations as in Eq. (6) raises the tween the electron wavefunctions. This phenomenon is energy of the rhombic phase, which therefore disappears illustratedinFig. 5,andcouldimplyatendencytowards fromthephasediagramatsufficientlyhighdensity,leav- ananisotropic(or“striped”)liquidphase,whichisgener- ing the other transitions essentially unchanged. allynotruledoutbytheisotropicnatureoftheCoulomb Wehaveanalyzedthequantummeltingofthedifferent repulsion (see also the hatched regions in Figs. 3 and Wigner crystalstructuresrealizedin suchbilayersystem 4).37 The results reported in Fig. 5 also indicate a pos- by making use of the Lindeman criterion discussed in sible reentrant behavior, although no conclusive answer the preceding Section. We see from Fig. 6 that both can be given at this level of approximation. the sequence of phases and the critical melting densities obtained within the present quadratic approximation are in satisfactory agreement with the more sophis- E. Symmetric electron bilayer ticatednumericalresultsofRefs.21,22,23(themeltingdensity is slightly underestimated as compared with QMC, In this section we analyze a system composed of two but quite similar to the DFT result). It is interesting to coupledelectroniclayers,inordertocheckthevalidityof see that the same trends observed in the preceding Sec- ourapproachbydirectcomparisonwithavailableDensity tionforthelayeredsolidsarealreadypresentinthesingle Functional Theory21 and Quantum Monte Carlo based bilayer. In particular, reducing the interlayer separation calculations22,23. Intheearlyworkonclassicalbilayers,19 leads to a sensible stabilization of the crystal compared 8 to the isolated layers. This is clear in Fig. 6, where the the melting line always lies below the critical value 50 rc 40 of a purely 2D Wigner crystal. Note also s,2D ≃ II that, contrary to Ref.22, we find that the enhancement 40 I of Wigner crystallization is slightly more pronounced in aninfinitearrayoflayersthaninasingleelectronbilayer. 30 BCT D 1 s, III. WIGNER CRYSTALLIZATION IN QUASI r 20 ONE-DIMENSIONAL SOLIDS 10 We now extend our analysis to the case of quasi one- dimensional solids, which we model as periodic arrays Liquid of conducting wires. Following the general arguments 0 presented in the previous Section, the enhancement of 0 10 20 30 40 50 Wigner crystallization in this case should be even more d pronounced than in the two-dimensional case, because FIG.7: Phasediagramforathree-dimensionalWignercrystal of the suppression of electronic motion in two transverse embeddedin asquarearray of1-dimensional wires, ofside d. directions rather than one. The effect is even more dra- Lengthsarescaledtoeffectiveatomicunits. Forthedefinition matic if we consider that a quantum crystal with gen- ofthephasesIandII,seetext. Theshadedregioncorresponds uine long-range order cannot be realized in a pure one- to γ >10 (not studied). dimensionalsystem,39 while it is stabilizedif we account forthelong-rangeCoulombinteractionsbetweencarriers on different wires.44 We shall consider here a square array of wires for il- a square array of wires, has the lowest Madelung energy lustrative purposes, although the specific arrangements in the whole range 0 < γ < 2.83, with the two special occurring in real solids (rectangular, rhombic) can be values γ = √2 and γ = 2 corresponding respectively ∗ ∗ treated case by case. Assuming a simple ordering of pe- to a BCC and a FCC. For γ > 2.83 the minimum con- riod a within the wires and an interwire distance d, the figuration becomes less symmetric, with c/a = 1/2 but 6 most general elementary three-dimensional Bravais lat- the ratio b/astill lockedto the value 1/2up to γ =6.99. tice compatible with the given geometrical constraint is This phase is denoted (I) in Fig.7. Beyond γ = 6.99, a described by the following basis vectors: Aˆ = (0,0,a), second structural transition occurs leading to a generic 1 Aˆ = (d,0,b), Aˆ = (0,d,c). The volume of the 3D uni- phase (II) with both b/a = 1/2 and c/a = 1/2. Other 2 3 6 6 tary cell is V = ad2 4πr3/3, the anisotropy ratio is transitionscantakeplaceatlargervaluesofγ,withinthe c ≡ s now defined as γ =a/d and the 1D density parameter is genericphaseII.Thesequenceofphasesdoesnotchange r = a/2. As in the layered case, we take a compen- upon inclusion of the vibrational term. s,1D sating positive charge distributed uniformly in the bulk. ByapplyingtheLindemannruleweobtainaparamet- The analysis presented in the preceding Section can be ric formula for the melting curve analogous to Eq. (13): repeated here following the same steps: i) calculation 2/3 of the structure with the lowest Madelung energy upon 1 M (γ) varying the anisotropy ratio; ii) calculation of the corre- rsc,1D = (128π)1/3 −δ12 d2/3 (14) (cid:20) (cid:21) sponding vibrational energies; iii) determination of the melting curve via the Lindemann criterion. The gener- with the implicit condition γ = rc /2d. The conse- s,1D alization is straightforward,and we only report here the quencesofgeometricalconfinementevidencedinthe lay- main results. ered case are recovered here. The electron spread along Thestucturalphasediagram(Fig. 7)isclearlylessrich the wires is again governed by the three-dimensional thaninthelayeredcase,becauseoncethedensityandthe plasma frequency [cf. Eq. (12)], due to the isotropic interwire distance d are fixed, only the relative ordering nature of the Coulomb interactions, while the nearest- between the electronic crystals on neighboring wires re- neighbor distance here scales with r (2π/3)r3/d2. s,1D ≡ s mainstobe determined,correspondingtothe pairofpa- Further stabilization of the crystallized state is achieved rametersb andc. Inthe limit γ 0,the interwireinter- through a reduction of the electron spread along the → actionsvanishandthelimitofisolatedwiresisrecovered: wires, revealed by an inverse moment M which is 1 the Madelung constant A diverges due to the isotropic typically 50% lower than the value in vacu−um. Its γ- distribution of the jellium, as explained previously (cf. dependence for γ & 1 is quite flat (not shown), except footnote 43). In this limit the interwire ordering is stag- in the vicinity of the transition at γ = 2.85, where it gered, with b/a = c/a = 1/2, corresponding to a body raises due to the mode softening discussed in Section II centered tetragonal (BCT) lattice in three-dimensional C. Replacing the average value M 2 into Eq. (14) 1 space. The BCT structure, everywhere compatible with yieldsrc 1.2d2/3,correspondin−gt≃oanevenstronger s,1D ≃ 9 enhancement of Wigner crystallization than in the lay- isolated units is recovered at large separations (γ 1), ≪ ered case (see Table I). In the opposite anisotropic limit an overall isotropic ordering of the charges is achieved γ 1, M diverges as in the case of an isolated wire for γ & 1, when the interactions between different units 1 (cf≪. footn−ote 43), so that the Wigner crystal is never becomeimportant. Inthiscase,three-dimensionalstruc- stabilized (rc ). tures as close as possible to the ideal case of a BCC s,1D →∞ areformed, leadingto acascadeof structuraltransitions ∗ ∗ which can be tuned by varying the particle density, or γ crystal melting d=8a d=20a B B the distance d itself. In addition to this rich phase di- layers √πr /d rc 4.9 d1/2 γ &1 14 21 s,2D s,2D ≃ agram, the presence of isotropic Coulomb interactions 40 γ 1 in such anisotropic compounds results in a strong sta- ≃ ≪ wires 2rs,1D/d rsc,1D ≃1.2 d2/3 γ &1 5 9 bilization of the charge ordered phases, possibly up to γ 1 densities of practical interest, where the characteristic →∞ ≪ energy scales of the Wigner crystal can become compa- TABLE I: Definition of the anisotropy ratio γ, approximate rable with other relevant scales in the solid. Although it melting lines obtained for quasi two-dimensional and quasi is clearthat the interplaywith severalother factors such one-dimensional systems, and specificvalues obtained at two as the periodic lattice potential,6,16,17,18,30,38 chemical different interlayer(interwire) distances d, expressed in units impurities,40 polarons15,27 or magnetic interactions14,41 ∗ of the effectiveBohr radius a (right columns). B should be considered for an accurate description of real materials, the long-range Coulomb interactions appear in light of the present study as a key ingredient to un- derstand the charge ordering phenomena in quasi low- IV. CONCLUSIONS dimensional systems. We have investigated the Wigner crystallization of ACKNOWLEDGMENTS electrons in quasi low-dimensional compounds, where the carrier motion is effectively low-dimensional, while the Coulomb interactions are assumed long-ranged and isotropic. The system properties are found to depend WethankS.Ciuchiforcriticalandconstructivediscus- cruciallyon the ratio γ ofthe meaninterparticlespacing sions. G.R.thanksthekindhospitalityofCNRS-LEPES withintheconductingunits (layersorchains)tothe sep- Grenoble (France) and finantial support by MIUR-Cofin arationdbetweenunits. Whilethebehaviorexpectedfor 2004/2005matching funds programs. 1 E. P. Wigner, Phys. Rev. 46, 1002 (1934), E. P. Wigner, (2005) Trans. Faraday.Soc. 34, 678 (1938). 17 H. Falakshahi et al., Eur. Phys. J. B 39, 93 (2004) 2 R.S.Crandall, R.Williams, Phys.Lett.A 34404 (1971). 18 B. Valenzuela, S. Fratini and D. Baeriswyl, Phys. Rev. B 3 C.C.Grimes,G.Adams,Phys.Rev.Lett.42,795(1979). 68, 045112 (2003) 4 E. Y. Andrei, Two-Dimensional Electron Systems on He- 19 G. Goldoni, F. M. Peeters, Phys. Rev.B 53, 4591 (1996). lium and other Cryogenic Substrates, (Kluwer Academic 20 I. V. Schweigert, V. A. Schweigert, F. M. Peeters, Phys. Publ. 1997), pags. 245-279 and refs. therein. Rev. Lett 82, 5293 (1999); I. V. Schweigert, V. A. 5 J. Yoon, C. C. Li, D. Shahar, D. C. Tsui, M. Shayegan, Schweigert, F. M. Peeters, Phys.Rev.B 60, 14665 (1999) Phys. Rev.Lett. 82, 1744 (1999). 21 G.Goldoni,F.M.Peeters,Europhys.Lett.37,293(1997). 6 J. Hubbard,Phys. Rev.B 17, 494 (1978). 22 L. S´wierkowski, D. Neilson and J. Szyman´ski, Phys. Rev. 7 F. Nad, P. Monceau, C. Carcel, and J. M. Fabre, Phys. Lett. 67, 240 (1991) Rev. B 62, 1753 (2000); P. Monceau, F. Ya. Nad and S. 23 G.Senatore,F.RapisardaandS.Conti,Int.J.Mod.Phys. Brazovskii, Phys. Rev.Lett. 86, 4080 (2001) B 13, 5-6, 479 (1999) 8 D. S.Chow et al., Phys. Rev.Lett. 85, 1698 (2000) 24 W. J. Carr, Phys. Rev. 122, 1437 (1961). 9 K.HirakiandK.Kanoda,Phys.Rev.Lett.80,4737(1998); 25 W.J.Carr,R.A.Coldwell-Horsfall, andA.E.Fein,Phys. Y. Nakazawa et al., Phys.Rev.Lett. 88, 076402 (2002) Rev.124, 747 (1961). 10 P. Abbamonteet al., Nature 431, 1081 (2004) 26 C. B. Clark, Phys. Rev. 109, 1133 (1958). 11 A. Rusydiet al., cond-mat/0511524 27 S. Fratini, P. Quémerais, Eur. Phys. J. B 14, 99 (2000); 12 P. Horsch,M. Sofin,M. Mayr, and M. Jansen, Phys.Rev. idem,Eur.Phys.J.B29,41(2002);G.Rastelli,S.Ciuchi, Lett. 94, 076403 (2005) Phys. Rev.B 71, 184303 (2005) 13 M. Dressel and N. Drichko,Chem. Rev.104, 5689 (2004) 28 P. P. Ewald, Ann.Phys. 64, 253 (1921). 14 H. C. Fu, J. C. Davis, D.-H. Lee, cond-mat/0403001 (un- 29 E. Ya. Sherman,Phys. Rev.B 52, 1512 (1995) published) 30 E.V.TsiperandA.L.Efros,Phys.Rev.B57,6949(1998) 15 G. Rastelli, S. Fratini, P. Quémerais, Eur. Phys. J. B 42, 31 E.L.Pollock,J.P.Hansen,Phys.Rev.A8,3110(1973);R. 305 (2004) Mochkovitch, J. P. Hansen,Phys. Letts. A 73, 35 (1979). 16 D. Baeriswyl and S. Fratini, J. Phys IV France 131, 247 32 L. Bonsall, A. A. Maradudin, Phys. Rev. B 15, 1959 10 (1977); R. C. Gann, S. Chakravarty G. V. Chester, Phys. high momentum states available (T = 1.11/r2, 0.5/r2 s s,2D Rev.B 20, 326 (1979). and 0.11/r2 in 3D, 2D and 1D respectively) while the s,1D 33 K. Moulopoulos, N. W. Ashcroft, Phys. Rev. B 48, 11646 Madelung energy in the opposite classical limit is less de- (1993) pendent on dimensionality [E = 0.896/r for the body 3D s 34 Yu. G. Pashkevich and A. E. Filippov, Phys. Rev. B 63, centeredcubiclattice,E =1.11/r forthehexagonal 2D s,2D 113106 (2001) lattice]. Incidentally, the crystallization transition in both 35 H.Nagara,Y.NagataandT.Nakamura,Phys.Rev.A36, 3D(r 100)and2D(r 40)takesplaceatthesame s s,2D 1859 (1987). valueof≃theratio E/T 81.≃ 36 G.Rastelli,P.QuémeraisandS.Fratini,J.Phys.IVFrance 43 Theusualvalueforthe≃Madelungenergyofthehexagonal 131, 277 (2005) lattice in two dimensions would be recovered at large in- 37 A.A.Slutskin,V.V.Slavin,andH.A.Kovtun,Phys.Rev. terlayerseparationsbyconsideringalayeredcompensating B 61, 14184 (2000) charge.Thiswouldnotalterthesequenceofcrystalstruc- 38 Y. Noda, M. Imada, Phys. Rev.Lett. 89, 176803 (2002) tures, since it would add a constant term to the energy 39 H. J. Schulz,Phys. Rev.Lett. 71, 1864 (1993) dependingonly on the anisotropy ratio γ. 40 S.T.ChuiandB.Tanatar,Phys.Rev.Lett.74,458(1995) 44 This can be understood by observing that the logarith- 41 B.P.Stojkovic,Z.G.Yu,A.R.Bishop,A.H.CastroNeto micdivergenceofthemeanspread u2 dωω−1 arising h i∼ and N.Grnbech-Jensen,Phys.Rev.Lett.82, 4679 (1999); in 1D is washed out in the present case where the wires V. M. Pereira, J. M. B. Lopes dos Santos, A. H. Castro areembeddedinathree-dimensionalenviroRnmentandmo- Neto, cond-mat/0505741 mentum integrations run overthe3D Brillouin zone. 42 Thekineticenergyevaluatedinthenoninteractinglimitis strongly reduced inlower dimensions, wherethereareless