Free energy analysis for a system of interacting particles arranged in Bravais lattice
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] J u l Free energy analysis for a system of interacting particles arranged in Bravais lattice
B. I. Lev, V. B. Tymchyshyn, and A. G. Zagorodny
Bogolyubov Institute for Theoretical Physics, National Academy of Science, Metrolohichna St. 14-b, Kyiv 03680, Ukraine
We propose a method of free energy calculation for a system of interacting particles arranged ina Bravais lattice. It will be shown how to treat divergences for infinite unbounded systems with“catastrophic” potentials like Coulomb and compare two of them. Besides, we show that currentmethod may be used not only for essentially classical systems, but for some quantum as well.Two systems are considered: grains in dusty plasma and electrons on the liquid-helium surface.For dusty-plasma parameters of grain’s lattice are calculated by numerical solution of obtainedequations. Electrons on the liquid-helium surface are analyzed to get localization distance. Besides,conditions for existence of electron lattice are found.
PACS numbers: 52.27.Lw, 71.10.Pm, 02.30.Mv
INTRODUCTION
There are many soft-matter systems, such as grainsin dusty plasma, colloids in various solvents, surfactantsolutions, etc., that exhibit self-organization and rear-rangement in crystalline structures under certain condi-tions. Often their inter-particle potential is long-rangeand similar to the Coulomb one. Usually, this kind ofpotentials cause divergences (that’s why they are oftencalled “catastrophic”) and complicates consideration [1–4]. On the other hand, mentioned systems are very inter-esting due to their application to the studies of a varietyof peculiar phenomena in different fields of science [5–7].It seems to be the most challenging problem to treatCoulomb-like systems with high concentrations of inter-acting particles [8]. When concentration increases, onecan observe crystallization-like phenomena and transi-tion between different lattice symmetries [9–13].In current contribution we will treat concrete systemsalong with general formalism development. They werechosen as manifesting a lot of interesting effects. For ex-ample, dusty plasma may serve as perfect media for theexperimental investigation of classical fluids and solidsalong with colloids [13–22]. The second system is “elec-trons on the liquid-helium surface” [23–26]. It may un-dergo Wigner crystallization [27–33] and transition be-tween triangular and square lattices [34, 35]. More-over, homogeneous distribution of electron density is notalways stable, and there are critical parameters whenspatial structures, especially periodic deformations andmulti-electron dimples, are formed [36].But our original goal is some general description ofparticles arranged in a Bravais lattice when interparticlepotential may be “catastrophic”. It is rather difficult,since traditional methods of statistical mechanics can-not be applied to systems with Coulomb-like interactionsbecause of free-energy divergences. Moreover, partitionfunction can be exactly evaluated only for few modelsystems of interacting particles in the thermodynamicallimit [37–40]. Thus, basing on [35, 41–43], we will use al-ternative approach that employs Hubbard-Stratonovichrepresentation of the partition function [44]. Besides, we use approach involving probability distribution functionfor the whole lattice as in [1], but construct it from single-particle distribution functions with a special procedure.Namely, we will expand and significantly supplement re-sults obtained in [35] for two-dimensional systems. Newresults for three-dimensional case, that is far more com-plicated than two-dimensional, will be obtained.To be more rigorous: our goal is to find expressionfor free energy of particles arranged in a Bravais latticeeven if inter-particle potential is “catastrophic”. Gen-eral equations will be applied to two concrete systems:three-dimensional — grains in dusty plasma, and two-di-mensional — electrons on the liquid-helium surface. Infirst case we’ll minimize this energy to find out whichlattice should be observed experimentally. And for elec-trons localization distance will be found. Besides, thiswill show that presented method allows easily reduce di-mension number from three to two.Present article is organized as follows. It can be contin-gently divided into two parts: development of techniqueand its application. Sections I and II focus on generalexpression for free energy for any Bravais lattice and anypotential. The former shows that not only classical, butlots of quantum systems may be treated this way. Be-sides, an expression for “entropy part” of free energy iscalculated. And in section II interaction energy for pari-cles arranged in a Bravais lattice is calculated.Both this sections are equipped with examples (I A,II A–II C) that show how to use developed formalism.This examples are rather useful, because they are laterreused when treating physical systems.Concrete applications of developed formalizm startwith section III. Here crystal of dust particles is con-sidered. In III it is shown that hexagonal close packingseems to be the lattice we expect to be seen in exper-iment. This is in agreement with experiments [9] andcomputer simulations [10].System of electrons on the liquid-helium surface is an-alyzed in IV. In IV localization distance is calculated andcompared to classical result [27]. For low temperatures T → I. STATISTICAL DESCRIPTION OFINTERACTING PARTICLES SYSTEM
First of all let’s outline class of systems we will beable to consider. Since everything is clear with classicalsystems, we will give heed to quantum ones.In this section we aim some general consideration of in-homogeneous system of interacting particles [35, 42]. Norsystem, neither interaction potential will be specified tokeep generosity and extendability. As result, we will getfree energy in some general form and ready for “apply-ing” translation symmetry (section II). Of course, alongthis section we will do some assumptions about systemunder consideration. Together they will form restrictionson quantum systems we can treat with proposed method.In current model macroscopic states of the system aredescribed by occupation numbers. Besides, we supposethat it’s Hamiltonian has form: H = X s ε s n s + 12 X s,s ′ V ss ′ n s n s ′ . (1)Here ε s is the additive part of particle energy (usuallyit is kinetic energy, but it can be particle’s energy inexternal field as well.), s indicates particle state, V ss ′ isinteraction energy between particles in states s and s ′ , n s is the occupation number of state s . We neglect quan-tum correlations, so only essentially classical systems areconsidered.Partition function for this kind of system will be: Z = X { n s } exp( − βH ) , where summation is performed over all possible states { n s } of the system.Now we use some properties of Gaussian integrals overauxiliary fields, i.e. Hubbard-Stratonovich transforma-tion [44, 45]:exp ν ϑ X s,s ′ ω ss ′ n s n s ′ == ∞ Z −∞ Dϕ exp ν X s n s ϕ s − ϑ X s,s ′ ω − ss ′ ϕ s ϕ s ′ ! , with Dϕ = Q s dϕ s / p det(2 πβω ss ′ ). Second-order de-pendence on occupation numbers is avoided by carryingthem to introduced field and partition function can bewritten: Z = ∞ Z −∞ Dϕ exp X s ( iϕ s − βε s ) n s − β X s,s ′ (cid:0) V − ss ′ ϕ s ϕ s ′ (cid:1) . For canonical ensemble number of particles can be fixedby using Cauchy equation:12 πi I ξ P s n s − N − dξ = 1 , that leads to N -particle partition function: Z N = 12 πi I dξ Z Dϕ exp " − β X s,s ′ V − ss ′ ϕ s ϕ s ′ −− ( N + 1) ln ξ s X { n s } [ ξ exp( iϕ s − βε s )] n s . Performing summation over occupation numbers: Z N = 12 πi I dξ Z Dϕ exp[ − βF ( ϕ, ξ )] ,βF ( ϕ, ξ ) = 12 X s,s ′ V − ss ′ ϕ s ϕ s ′ ++ δ X s ln (cid:0) − δξe − βε s + iϕ s (cid:1) + ( N + 1) ln ξ. (2)Here type of statistics is incorporated into δ — it equalsto +1 for Bose-Einstein, 0 for Maxwell-Boltzmann and − ξ = exp( βµ ) for chemical activ-ity. Equation (2) contains all information about probablestates of the system and corresponds to the sequence ofequilibrium states with regard to their weights. We areinterested in asymptotic value of partition function, butwant to avoid using perturbation theory. Thus domainis extended to complex plane and saddle-point method isapplied.Dominant contribution is made by states satisfying ex-trema condition: δβFδϕ = δβFδξ = 0 . Variation of (2) yields expression for saddle-point states:1 β X s ′ V − ss ′ ϕ s ′ − iξe − βε s + iϕ s − δξe − βε s + iϕ s = 0; (3a) X s ξe − βε s + iϕ s − δξe − βε s + iϕ s = N + 1 . (3b)For a certain state expression f s = ξe − βε s + iϕ s − δξe − βε s + iϕ s , (4)from (3b) can be treated as an average occupation num-ber. Now saddle-point states (that may be interpreted asthermodynamically stable distributions) can be obtained.Equation (3a) contains inverse matrix, that is quiteinconvenient for us, because its calculation is a chal-lenging mathematical problem even for simple poten-tials [42, 46, 47]. This problem can be avoided if weuse (4) and perform inverse transformation ϕ s = iβ X s ′ V ss ′ f s ′ , β X s,s ′ V − ss ′ ϕ s ϕ s ′ = − β X s,s ′ V ss ′ f s f s ′ . Than free energy can be written βF [ f, ξ ] = − β X s,s ′ V ss ′ f s f s ′ − δ X s ln(1 + δf s )+ ( N + 1) ln ξ ( f ) . (5)In terms of canonical ensemble we can write from (4)ln ξ ( f ) = 1 N X s f s [ β ( ε s + E s ) + ln f s − ln(1 + δf s )] ,E s = X s ′ V ss ′ f s ′ , and substitute last expression into free energy (5) torewrite its expression without chemical potential. F [ f ] = X s f s ε s + 12 X s,s ′ V ss ′ f s f s ′ ++ 1 β X s [ f s ln f s − ( f s + δ ) ln(1 + δf s )] | {z } F ent . (6)First two terms in (6) are kinetic and potential energyrespectively, the last one is a contribution due to entropy F ent (should be equal to zero if T = 0).Now let’s do a quick test of obtained results. If a grandcanonical ensemble with fixed chemical potential is con-sidered, we can get from (4) f s = 1 e β ( ε s − µ s ) − δ , (7)which is generalization of the well-known distributionwith a chemical potential µ s = µ − E s . It is obvious that the “saddle point” and “mean field” ap-proximations are equivalents in this case. If an ideal gas( µ s ≡ µ ) is considered, we can obtain classical statisticaldistributions.Before winding up this section let’s review the results.We aimed to find some general expression for free energyin a mean-field approximation. As the one, equation (6)should be pointed out (and in some sense (7) that can beused to obtain particles distribution in grand canonicalensemble). Besides, we have got some restrictions onthe class of systems we will be able to consider: specialform of Hamiltonian (1) and negligible quantum cross-correlations. A. Example: two-dimensional system of localizedFermi particles — electrons on the liquid-heliumsurface
Before we shift our attention to potential energy ofparticles arranged in a lattice, let’s do some example.With this subsection we aim two things: give some ideaon how to use the above formalism and secondly — cal-culate contribution due to entropy for electrons on theliquid-helium surface (will be used in IV).Electrons on the liquid-helium surface have two degreesof freedom only [27, 28]. Though we set δ = − d = 2 (system is two-dimensional). Thenequation (6) can be written βF [ µ ] = Z d ~p d ~r (2 π ~ ) β~p m e β ( ~p / (2 m ) − µ ( ~r )) + 1 ++ β Z Z d ~p d ~r (2 π ~ ) d ~p ′ d ~r ′ (2 π ~ ) ×× V ( | ~r − ~r ′ | ) (cid:0) e β ( ~p / (2 m ) − µ ( ~r )) + 1 (cid:1)(cid:0) e β ( ~p ′ / (2 m ) − µ ( ~r ′ )) + 1 (cid:1) −− Z d ~p d ~r (2 π ~ ) ln (cid:16) e β ( ~p / (2 m ) − µ ( ~r )) + 1 (cid:17) e β ( ~p / (2 m ) − µ ( ~r )) + 1 ++ ln (cid:16) e − β ( ~p / (2 m ) − µ ( ~r )) + 1 (cid:17) e − β ( ~p / (2 m ) − µ ( ~r )) + 1 . (8)Here integration sign R means integration over the wholephase space.Since last expression is quite cumbersome we shouldsimplify it somehow. As first approximation we suppose ε s = ε ( p ) = p / (2 m ) . Actually, in the presence of an external field, dispersionrelation should be a bit different. But we know, thatsystem of Fermi particles, we will be interested in, ishighly degenerated [27]. So quadratic form for dispersionrelation looks quite reasonable.Using last equation and introducing thermal length λ T = p π ~ β/m, (9)we get from (7) following expression ρ ( ~r ) = Z d ~p (2 π ~ ) e β ( ~p / (2 m ) − µ ( ~r )) + 1= πλ T ln (cid:16) e βµ ( ~r ) (cid:17) , (10)which connects chemical potential µ ( ~r ) and particlesdensity ρ ( ~r ).Now we return to simplification of (8). Performingintegration of (8) over momentum and taking (10) and(9) into account one may get (see appendix VI A) F [ ρ ] = 12 Z Z d ~r d ~r ′ V ( | ~r − ~r ′ | ) ρ ( ~r ) ρ ( ~r ′ )++ πβλ T Z d ~r (cid:20) Li (cid:16) e − π − λ T ρ ( ~r ) (cid:17) − π (cid:21) ++ λ T πβ Z d ~rρ ( ~r ) . (11)One may notice that in case of Bose statistics we willloose the Bose-condensation effects due to this integra-tion. But we aim to use this equations for Fermi particles,so there is nothing to be worried about.We know that electrons on the liquid-helium surfaceexist in forms of fluid or Wigner crystal [29–31]. Weare interested in case when electrons are strongly local-ized, i.e. their coordinates are quite determined. Butsince temperature differs from zero, we expect that par-ticle’s position will fluctuate near its equilibrium loca-tion. Though we assume following distribution function ρ sp for s ingle p article (somewhat analogous to form fac-tor in QFT) ρ sp ( ~r ) = 1(2 πs ) d/ e − r / ( s ) , (12)where s is dispersion, or localization distance by physicalmeaning, and d = 2 because system is two-dimensional .Equation (12) is normal distribution which seems tobe quite reasonable. At least this distribution is validfor ground state of quantum harmonic oscillator which isthe simplest approximation for particle fluctuating nearits equilibrium location.Last assumption we do is that “Gaussian” (12) is very“sharp” and every particle lays outside localization radiusof any other particle. This means that s is much smallerthen other characteristic distances in this system. So freeenergy per one particle can be written: F sp = 12 Z Z d ~r d ~r ′ V ( | ~r − ~r ′ | ) ρ sp ( ~r ) ρ ( ~r ′ )++ λ T π s β − π s βλ T . | {z } F ent (13)Here F ent is contribution to free energy by ent ropy and temperature. First summand in F ent is obtainedfrom direct integration of ρ (see last term in (11)). Toget second summand (penultimate term in (11)) we use“sharpness” of (12) and approximateLi (cid:16) e − π − λ T ρ sp ( ~r ) (cid:17) ≈ ( , r ≤ s ; π / , r > s. (14) We aim to reuse this and some of subsequent equations forsystems with different number of dimensions. In such cases wewill write them as functions of d (dimensionality) even if it isredundant at the point we introduce the equation first time. With equations (13) and (11) we finish our considera-tion. We fulfilled our goal to show an example of usingsection’s I formalism. Moreover, now we have expressionfor entropy part F ent of free energy (13), which we uselater in IV. II. PARTICLES ARRANGED IN A LATTICE
In section I we found some expressions for free energyof inhomogeneous system of interacting particles (6). Inthis section we will be interested in second term of (6),i.e. potential energy. More accurately, we aim to findmethod of potential energy calculation for system of iden-tical (from inter-particle’s potential point of view) parti-cles arranged in a Bravais lattice. During this section wewill show how to use translation symmetry to get in somecases more convenient expressions for potential energy.Last remark before we dive in. We aim to find appro-priate expressions for both two- and three-dimensionalcases. But to avoid reduplication all along the section IIonly three-dimensional case will be treated as far morecomplicated. For two-dimensional case we will only pro-vide results and short description where difference comesfrom.Let V designates one subset of exact cover of R (fig-ure 1) with similar domains containing strictly one par-ticle (relative position of the particle supposed to be thesame in every domain) . Then free energy of the particlelocalized near Bravais lattice site can be written as (weneglect kinetic energy in (6) since particles are localizedand potential energy should be bigger) F int = Z Z Z V Z Z Z R V ( | ~r − ~r ′ | ) ρ ( ~r ) ρ ( ~r ′ ) d ~r ′ d ~r, (15a) F s = Z Z Z V Z Z Z V V ( | ~r − ~r ′ | ) ρ ( ~r ) ρ ( ~r ′ ) d ~r ′ d ~r, (15b) F sp = F int − F s + F ent . (15c)Equation (15) needs some explanations. Interactionenergy is presented here with F int − F s . This “splitting”is just another way of saying we want to integrate over R \ V . If one collects both expressions together and uses RRR R \ V = RRR R − RRR V he will get regular integrationover R \ V for ~r ′ .In other words, since particle is localized near latticesite, we take some “part of space” V around it. Thenby integration we calculate interaction energy betweenparticle in V and all the rest particles in R \ V .One may get a good intuitive idea about integrationdomain R \ V through discrete case. When we calculate For example, Wigner-Seitz cells provide this kind of covering.But they often have too complicated shape and we will use par-allelepipeds instead (Figure 1).
Coulomb energy for i -th particle we write P j = i e i e j /r ij .Integration over R \ V is somewhat analogous to j = i for discrete case.We will further call F int (15a) int eraction energy and F s (15b) s elf-interaction energy . This “splitting” hasone more convenience. Basing on physical backgroundof considered potential we may want or not want to com-pensate self-interaction. Obviously, Coulomb potentialshould have self-interaction compensated — potential en-ergy is zero if there is only one particle in a system(considered later in example II B). On the other hand,self-interaction for effective potentials caused by influencethrough medium should be uncompensated, because evenin absence of other particles medium will influence thisone. Surface distortion is a good example. If particle hasdistorted a surface it has already done some contributionto system’s free energy independently on other particles(considered later in example II C). Of course, it interactswith distortions created by other particles. But at thesame time it interacts with “its own” surface distortionand this may be treated as “self-interaction”.Now lets consider (15) more carefully. All 3D Bra-vais lattices can be treated as constructed from paral-lelepipeds. We will use this property and use paral-lelepiped as integration region V . y xz ~a~b ~cα β c α c Figure 1: One cell from exact cover V . Vectors ~a , ~b and ~c are base vectors for Bravais lattice. Angle between ~a and ~b is supposed to be α . Angle between ~c and XY plane is β c and its projection on XY plane and ~a is α c .One particle is shown in the center.Next thing we want to do is expressing probabilitydistribution function ρ for the whole R through single-particle probability distribution ρ sp . Basing on the in-verse lattice vectors (appendix VI B) we will provide fol-lowing decomposition ρ ( ~r ) = X ~k ∈ Z ρ ~k f ~k ( ~r ) , (16a) ρ ~k = ¯ ρ Z Z Z V f ∗ ~k ( ~r ) ρ sp ( ~r ) d ~r, (16b) f ~k ( ~r ) = e π i ( ~k T ˆ G~r ) , (16c) ˆ G = a − cot( α ) a − cot( β c ) sin( α − α c ) a sin( α )0 csc( α ) b − cot( β c ) sin( α c ) b sin( α )0 0 csc( β c ) c . (16d)Here ∗ designates complex conjugation. It may be no-ticed that all vectors are treated as column vectors. ¯ ρ is mean particle density here. To get better intuitiveunderstanding how to compute (16) one may check ex-ample II A for ρ sp from (12).It may be shown that presented series are actuallyFourier series for ρ sp and ρ is periodic along vectors ~a , ~b and ~c (appendix VI B). So we may treat ρ ( ~r ) as functionon R “composed” from single particle distributions “ar-ranged in a lattice”. This principle is demonstrated byFig. 2. Since ρ has the same symmetry as the lattice andFigure 2: Two-dimensional analogy of transition fromsingle-particle density probability function ρ sp ( ~r ) tomany particle ρ ( ~r ). Function ρ ( ~r ) is defined for all R with two-dimensional analogy of (16).locally is a good approximation for ρ sp it will be treatedas probability distribution function for the whole lattice.Matrix ˆ G comprises all g eometry of the lattice underconsideration. Performing minimization of F sp over en-tries of ˆ G , one may find which lattice can be formed undercertain conditions. We will use later this type of analysisfor dust crystal in section III.If we use equations (16) we can express (15a) in termsof ρ sp (appendix VI D) F int = 1¯ ρ X ~k ∈ Z | ρ ~k | V ~k , (17a) V ~k = ¯ ρ Z Z Z R f ~k ( ~r ′ ) V ( | ~r ′ | ) d ~r ′ , (17b)where ¯ ρ is mean particle density.Before moving further, one more thing should be em-phasized. This equations are useful because they allowus to deal with “catastrophic” potentials like Coulomb.One may notice that for any descending potential only1¯ ρ | ρ | V = 4 π ¯ ρ ∞ Z V ( r ) r dr, can be diverging and only if potential is “catastrophic”.But this term does not contain “lattice geometry”, itdepends on mean particle density only. This means,we can compare two lattices with equal mean particledensity even if inter-particle potential is “catastrophic”(somewhat analogous to renormalization in QFT).With equations (17), (16) and (15) we achieve our goaland finish consideration of interaction energy for parti-cles arranged in a lattice. Following examples show howto use obtained equations for screened Coulomb (exam-ple II B) and capillary interaction (example II C). Theyboth use ρ obtained with ρ sp from (12) (example II A). A. Example: calculating distribution function ρ forGaussian single-particle distribution ρ sp We suppose that all particles are arranged in a lat-tice and their positions are quite determined. But sincetemperature differs from zero we can expect that particleposition fluctuates near its lattice site. Using considera-tions from I A we suppose it to be (12) with d = 3.Introduced constant s means dispersion of presenteddistribution. From physical point of view it can betreated as localization distance. Later, in section IV,we’ll find its value and compare it with localization dis-tance obtained in different works.Using that ρ sp ( r ) is a very “sharp” Gaussian, we ex-pand integration limits to R in (16b) and perform inte-gration (appendix VI C) ρ ~k = ¯ ρe − π s ~k T ˆ G ˆ G T ~k . (18)As mentioned previously, we are interested in two-dimensional case as well. One may check, that thereexists some sort of analogy of (16) for two dimensions[35]. Difference will appeare in changing summations andintegrations P Z → P Z , RRR R → RR R and usage of dif-ferent matrix ˆ G = /a − cot( α ) /a α ) /b . (19)Performing calculations as in appendix VI C with (12)and d = 2, one will get ρ ~k = ¯ ρe − π s ~k T ˆ G ˆ G T ~k , (20)where ¯ ρ is two-, not three-dimensional density as in (18).Equations (18) and (20) are aimed result of this exam-ple. Later they will be used in calculations of potentialenergy for Coulomb (II B) and capillary (II C) interac-tions. B. Example: calculating F int − F s for screenedCoulomb potential In section III grains in dusty plasma will be consideredand in section IV we’ll set our eyes on electrons on the liquid-helium surface. In both cases we have to deal withCoulomb potential — screened for dust particles and reg-ular for electrons. The latter will be treated as screenedCoulomb with screening distance approaching infinity.During this example we aim to find F int − F s for screenedCoulomb potential using methods of section II. Consid-eration of physics of this two systems will be postponeduntil sections III and IV.Screened Coulomb potential can be written as V ( r ) = q e − r/λ D r , (21)where λ D is Debye screening distance. It may be shown(appendix VI E) that expression for V ~k should be V ~k = 4 π ¯ ρq /λ D + 4 π ~k T ˆ G ˆ G T ~k . (22)At this point we already have V ~k (22) and ρ ~k (18). Sowe can write F int from equation (17a) F int = X ~k ∈ Z e − π s ~k T ˆ G ˆ G T ~k π ¯ ρq /λ D + 4 π ~k T ˆ G ˆ G T ~k . (23)Compared to straight-forward method of computing, thisone has better convergence. For big values of λ D we willsee, that first summands in (23) descend as 1 / | ~k | . Butif we write series containing interaction energy betweenthis particle and any other, first summands will descendonly as 1 / | ~k | .In section II we have mentioned that there is one com-pensatory term we need to take into account. This termwill exclude Coulomb self-interaction. If we omit techni-cal details (see appendix VI F) explicit form of F s can beobtained F s = 12 √ πs ∞ Z e − r ′ / (4 s ) V ( r ′ ) r ′ dr ′ , (24)or substituting screened Coulomb potential (21) as V andapproximating for the case when s is much smaller thenscreening distance λ D F s = q √ πs − q λ D . (25)Expressions we’ve got for F int and F s are a bit incon-venient. So we use that s is small compared to distancesin the lattice (“sharp peak” approximation, section II)and perform some approximation.First of all we introduce mean distance calculated fromparticle’s density l = 1 √ ¯ ρ , (26)and split sum (23) into two parts. Further we will sup-pose lattice under consideration as not too much degen-erated. In other words, we will assume that it may be de-scribed as some “deformation” of cubic lattice (appendixVI G) F int − F s = √ πq l X ~k =0 e − ~k T ˆ G − T ˆ G − ~k/l ~k T ˆ G − T ˆ G − ~k/l ++ q πl X ~k =0 e − l ~k T ˆ G ˆ G T ~k l ~k T ˆ G ˆ G T ~k − √ πq l + 4 π ¯ ρq λ D . (27)One may see that last term is equal to energy in caseof evenly distributed charge. Other terms are correctionsdue to charge pointness and comprise all information onlattice geometry.It is very useful to notice, that this is the only sum-mand that approaches infinity if λ D → ∞ . This meansthat (27) can be used for different lattice comparison evenwithout screening. We will only need to measure energystarting from 4 π ¯ ρq λ D as “ground-level”.One more thing to notice. Coulomb is a long-rangeinteraction and localization distance s should not playany significant role in approximation. Here we see thatit is really absent. Rather different may be a picturefor short-range interactions or effective interactions withuncompensated F s (example II C). C. Example: calculating F int for capillaryinteraction Last example we will consider is a rough estimation ofcapillary interaction energy. This results will be usedin section IV. More discussion on physics of the phe-nomenon and references we postpone until IV.It is known, that electrons deform liquid-helium sur-face when pressed against with electric field and thiscauses effective interaction with potential [34] V ( r ) = − q E πσ K ( r/l ) , (28)where qE is force caused by electric field, σ — sur-face tension of liquid helium, r — distance betweenparticles, K — modified Bessel function, and l = p σ/gρ He [34] — capillary length that depends on thefluid properties only.Since we know that capillary interaction is essen-tially two-dimensional (particles should be localized onsome surface, no three-dimensional analogue) we willuse (17) with (19) and (20). Performing calculations (ap-pendix VI H) we get V ~k = − q E σ ¯ ρ /l + 4 π ~k T ˆ G ˆ G T ~k . (29) So we can immediately write F int = − X ~k ∈ Z ¯ ρq E e − π s ~k T ˆ G ˆ G T ~k σ (cid:16) /l + 4 π ~k T ˆ G ˆ G T ~k (cid:17) . (30)For our purposes (see section IV) a very rough approx-imation would be enough. So we use the same trick as inappendix VI G, but perform integration instead of sum-mation and dealing with theta functions (see appendixVI I) F int = q E πσ ln ( s/l ) . (31)One should notice, that s ≪ l , so logarithm value isless then zero and thus whole expression (31) is less thenzero. It coincides with the fact that capillary interactionmakes particles attract each other [32–34]. III. GRAINS IN DUSTY PLASMA, DUSTCRYSTAL AND ITS LATTICE
It is known that grains in dusty plasma interact andeven show self-organization in form of melting and crys-tallization of dust crystalls [9–13].We aim to apply methods developed in section IIto this system. As a first approximation we supposethat grains interact only as charged particles (screenedCoulomb potential). Obtained equation (27) (exampleII B) allows to calculate energy of any lattice just insert-ing correct matrix ˆ G and performing summation. Butmuch more interesting is finding the lattice with minimalenergy when particles density ¯ ρ is fixed. Obtained resultcan be verified experimentally.Both sums in (27) are highly convergent and this isvery convenient for numeric calculations. Thus we willlimit summation with | ~k | ≤
4. This limit was foundexperimentally: taking more summands does not changecalculated minimal lattice and free energy. On the otherhand, taking less results in wrong lattice.And one more thing to mension. Searching for latticewith minimal energy does not need calculation of fullexpression for F int − F s . We need to minimize only partswith ˆ G . Moreover, they both contain q/l , which means,that results of minimization do not depend of particlecharge.Minimization was performed with standard functionFindMinimum of Wolfram Mathematica v.9. We use (27)as function to be minimized. Particles density is fixed ina following way. We consider1 l = ¯ ρ = 1 abc sin( α ) sin( β c ) , (32)suppose a = l (1+ δa ), b = l (1+ δb ) and express c throughother parameters to keep constant charge density. Per-forming this one may notice that l cancel out in expo-nents and denominators (27). ~b ~b ′ ~a~b ′′ ~a Figure 3: Same lattice in different representations. Ithas translation symmetry with respect to two vectors ~a and ~b . Upper figure shows different choice for secondtranslation vector ~b and ~b ′ . Using reflections one mayshow that lattice on bottom figure is the same as on theupper one. So ~b ′′ is possible choice as well. Whenperforming numerical calculations additionalrestrictions are used to overcome this uncertainty.Before considering results, one more remark should bedone. Mapping between set of all possible parameter val-ues (translation vectors) and all possible Bravais latticesis not a bijection. It is rather a surjection — one latticemay be described with different sets of parameters (dif-ferent translation vectors). For minimization algorithmsto be not confused, we provide some restrictions.Figure 3 shows uncertainty when choosing parametersfor lattice description for 2D case. One may see that itcan be eliminated if we claim that projection of ~b on ~a isless then a half of ~a . In terms of parameters this may beexpressed as 0 ≤ (1 + δb ) cos( α ) ≤ δa . (33)Restrictions for 3D case are more sophisticated. Butthey may be achieved similarly as (33). Actually (33)should be kept for 3D as is. We only need to add analo-gous considerations for ~c .First of all we introduce projection of vector ~c on XY plane (see figure 1) and take particles density into ac-count (substituting c from (32)) c x = c cos( β c ) cos( α c ) = l cot( β c ) cos( α c )(1 + δa )(1 + δb ) sin( α ) ,c y = c cos( β c ) sin( α c ) = l cot( β c ) sin( α c )(1 + δa )(1 + δb ) sin( α ) . Obviously this projection should fall into parallelogramconstructed on ~a and ~b . This condition is an obvious ex-tension of two-dimensional case to third dimension. But ~b ~c ′ − ~c ′ ~a ~b ~c − ~c~a Figure 4: Projection of the same lattice in differentrepresentations on XY plane. Particle from layer above XY plane is shown with circle and below the plane withfilled circle. Particles in XY plane are not explicitlyshown. If particle from above the plane is projected intothe upper half of ab parallelogram a reflection withrespect to XY plane can be used to get representationwhere this particle is in lower half of parallelogram.we can strenghten this condition by using reflection withrespect to XY plane (figure 4). We can always assumethat ~c is projected into lower part of parallelogram con-structed with ~a and ~b .As result we get following conditions0 < c y l < (1 + δb ) sin( α )2 c y cot( α ) l < c x l < c y cot( α ) l + (1 + δa )Now we can present results of numerical minimizationin form of figure 5 and following table.Figure 5: Part of lattice obtained from minimization.As is seen, 3-D lattice consists of interleaving planeswith 2-D hexagonal (triangular) lattices in each. Theonly difference between planes is shift: particles fromone plane project onto “empty spaces” in another.Translation vectors are shown as well.Looking at figure 5, one my see, it is very similar tohexagonal close packing (HCP) lattice. HCP consists oftwo interleaving planes with triangular 2-D lattice in eachas well. So we may check, if this is really HCP lattice bycomparing their parameters.Parameter Calculated Hexagonal closevalue packing (HCP) α ≈ . π/ ≈ . α c ≈ . π/ ≈ . β c ≈ . (cid:16)p / (cid:17) ≈ . a ≈ . √ ≈ . b ≈ . √ ≈ . IV. ELECTRONS ON THE LIQUID-HELIUMSURFACE AND THEIR LOCALIZATIONDISTANCE
In this section we consider electrons on the liquid he-lium surface, calculate localization distance and compareit to obtained in other works.In the presence of an external field electrons can bepressed against the helium surface with significant force.But they cannot go through this surface, because quan-tum effects push them out [27, 32, 33]. Thus they deformsurface and this deformation, in its turn, changes interac-tion potential between electrons. The problem of findingexplicit expression of this potential is solved by adding toCoulomb interaction effective capillary interaction. Thelateral for two electrons on the helium surface was calcu-lated in [34].We aim to find localization distance for electrons on theliquid helium surface. To achieve this goal we will mini-mize single-particle’s free energy F sp over s (localizationdistance, see example I A). So first of all we construct F sp and find ∂F sp /∂s .In section I we developed a method for treating somequantum systems. Example I A shows how to use it fortwo-dimensional system of Fermi-particles. So we con-sider F ent from (13) as free energy of one particle causedwith entropy in the system.Potential energy is more complicated. In section IIwe developed methods for its calculation. There aretwo parts of potential energy for electrons on the liq-uid helium surface [27, 32–34]: Coulomb interaction andcapillary interaction. For Coulomb interaction we canfind a shortcut — it was shown in [35] that this part ap-proximately does not depend on localization distance s .Besides, one may think on approximate expression forCoulomb interaction from example II B (27). It does notcontain s , and the same behavior we may expect fromCoulomb interaction for 2D lattice. Since we are inter-ested in derivative ∂F sp /∂s , consideration of this partof potential energy may be omitted. More interesting iscapillary interaction presented in example II C, approx-imation (31). It contains explicit dependence on s andthus should be taken into account. Bringing it all together we get ∂F sp ∂s = − λ T π βs − π s βλ T + q E σπs , independently on lattice we are considering.Substituting electron charge q → ¯ e and solving equa-tion ∂F/∂s = 0 we get s = 3 ~ β ¯ e E π σm − s − σ π e E β ! (34)If temperature tends to be small T → β → ∞ ) s = r πσ ~ m ¯ e E . (35)This result does not differ from the classical one pre-sented in [27].If we consider (34) more carefully, some interesting re-sult may be obtained. Let’s consider low densities so thatany value of s may be treated as much lower than inter-electron distance l . Then we may notice that solution for s exists if E T ≥ σπ k √ e ≈ dyncm · K . (36)Last equation means that there is no lattice if tem-perature T is too high, or if field E is too weak. Thisstatement is in agreement with our physical intuition. V. CONCLUSIONS
In conclusion we want to make some overview of resultswe have obtained.In section I we started with a statistical descriptionof inhomogeneous system of interacting particles. It wasshown that variety of quantum systems may be describedwith no much difference from classical ones. We foundan expression for free energy (6). The main feature ofthis equation appeares in expression for “entropy part”of free energy. Obtained result is applicable to a widevariety of systems with different number of dimensions.To make a step aside from pure theory we consideredan example I A of two-dimensional system of Fermi par-ticles. This result was later reused in section IV.In section II we stated that particles are arranged ina lattice. It seems, that adding translation symmetry toconsideration we can get some simplification for poten-tial energy expressions (16) and (17). One may noticethat for any “catastrophic” potential there is only onediverging term. But this term does not contain informa-tion on “lattice geometry”, it depends on mean particledensity only. This means, we can compare two latticeswith equal mean particle density even if inter-particlepotential is “catastrophic”.0In this approach all potentials are classified into twogroups: ones with self-interaction (often effective interac-tions through medium) and without. Examples II B andII C demonstrate how to use proposed technique in bothcases. Later they are reused for treating grains in dustyplasma (section III) and electrons on the liquid-heliumsurface (section IV).In section III we considered grains in dusty plasma.We started with potential energy obtained in examle II Band minimized it over all possible lattices. Since thereare no good approximations for theta-function it was toocomplicated to perform this minimization analytically.Thus we used numerical computations. This compu-tations have some features presented in a real-life caseIII. Results obtained in III show that we should expecthexagonal close packing (HCP) lattice for grains in dustyplasma. This coincides with numerical simulations in [10]and experiments [9].Last section IV shows how to use developed methodsfor electrons on the liquid-helium surface. It is shownthat in some cases much more simple (but more roughas well) compared to section III approximation may beobtained. This approximation is used to calculate lo-calization distance and then compared to classical result[27]. When temperature T → VI. APPENDICES
During this article we tried to provide rigorous expo-sition of presented ideas, but without physical sole beinglost in lots of equations. Thus, calculations that heavilyrely on mathematical transfomations only, were moved toappendices. Reader, interested only in “physical part”,may omit them.
A. Integration over momentum in βF [ µ ] . We start with (8). First two summands can be easilyintegrated if we mention1 e x + 1 = e − x e − x , but last summand should be integrated in terms of specialfunctions, namely dilogarithm [48]Li ( z ) = − z Z ln(1 − z ) z dz. If we mention that 11 + e x = 1 −
11 + e − x . we would get βF [ µ ] = m π ~ β Z Z d ~r d ~r ′ V ( | ~r − ~r ′ | ) ×× ln (cid:16) e βµ ( ~r ) (cid:17) ln (cid:16) e βµ ( ~r ′ ) (cid:17) ++ m π ~ β Z d r h βµ ( ~r ) ln (cid:16) e βµ ( ~r ) (cid:17) −− β µ ( ~r )2 − Li (cid:16) − e − βµ ( ~r ) (cid:17) − π (cid:21) , (37)which can be simplified further.First of all we will use Landen’s identity [48]Li (1 − z ) + Li (cid:18) − z (cid:19) = −
12 ln ( z ) ,z ∈ C \ ] −∞ ; 0] (38)with z = 1 + e − βµ ( ~r ) and expand ln (cid:0) e − βµ ( ~r ) (cid:1) . Thisleads to βF [ µ ] = m π ~ β Z Z d ~r d ~r ′ V ( | ~r − ~r ′ | ) ×× ln (cid:16) e βµ ( ~r ) (cid:17) ln (cid:16) e βµ ( ~r ′ ) (cid:17) ++ m π ~ β Z d r (cid:20) Li (cid:18)
11 + e βµ ( ~r ) (cid:19) − π (cid:16) e βµ ( ~r ) (cid:17)(cid:21) . Now we substitute (10) and (9) into the last expressionand get (11).
B. On the properties of f ~k ( ~r ) . In (16) we present how to calculate probability distri-bution function for the whole system basing on single-particle probability distribution function. It is based onthe decomposition of ρ sp using f ~k ( ~r ) (16c) as basis func-tions in inverse lattice space. Here we add some con-sideration about connection between Fourier series andpresented one.Suppose we have a function g ( ~r ) defined on the unitcube U ≡ [0; 1] × [0; 1] × [0; 1]. This function can beexpressed in terms of Fourier series g ( ~r ) = X ~k ∈ Z g ~k e π i ( ~k T ~r ) ,g ~k = Z Z Z U g ( ~r ) e − π i ( ~k T ~r ) d ~r. (39)Let’s consider a bijection ˆ G (16d) from V to U and itsinverse ˆ G − (see Fig. 1 for geometrical reasoning)ˆ G − = a b cos( α ) c cos( α c ) cos( β c )0 b sin( α ) c sin( α c ) cos( β c )0 0 c sin( β c ) (40)1It may be seen that rows of ˆ G are components of in-verse lattice basis vectors. We are interested in case g ( ~r ) = ρ sp (cid:16) ˆ G − ~r (cid:17) . Since we know that ∀ ~r ∈ V : ρ sp (cid:16) ˆ G − ˆ G~r (cid:17) = ρ sp ( ~r ) one can immediately write ρ sp ( ~r ) = X ~k ∈ Z g ~k e π i ( ~k T ˆ G~r ) . (41)And same way we consider second equation from (39) g ~k = Z Z Z V ρ sp (cid:16) ˆ G − ~r (cid:17) e − π i ~k T ˆ G ˆ G − ~r d (cid:16) ˆ G − ~r (cid:17) J h ˆ G − ~r i , (42)where J is Jacobian. Here ˆ G − ~r is treated as new vari-able with domain V .Jacobian J h ˆ G − ~r i = abc sin( α ) sin( β c ) = 1 / ¯ ρ is ac-tually the volume of V or inverse mean particle density.Designating g ~k as ρ ~k in (41) and (42), changing vari-able in (42) ˆ G − ~r → ~r and designating exponent in (41)as f ~k (it can be seen that exponent in (42) is f ∗ ~k , where ∗ designates complex conjugation) we immediately getequations (16). Moreover, obtained result means that allproperties of Fourier series can be applied to decomposi-tion (16).Last thing to mention are periodical properties of pre-sented series. We are going to expand domain to R sothat ρ ( ~r ) will be defined everywhere in the space. Sonow we need to explore behavior of this function. From(16d) one may see, if l ∈ Z , m ∈ Z and n ∈ Z :ˆ G (cid:16) ~r + l~a + m~b + n~c (cid:17) = ˆ G~r + l~e x + m~e y + n~e z , (43)where ~e x , ~e y and ~e z are unit vectors along coordinateaxis. With regard to (16c) f ~k (cid:16) ~r + l~a + m~b + n~c (cid:17) = f ~k ( ~r ) (44)and from (16a) ρ (cid:16) ~r + l~a + m~b + n~c (cid:17) = ρ ( ~r ) . (45)Last equation justifies our view of connection between ρ and ρ sp as it is presented on Fig. 2. C. On the expression of ρ ~k for Gaussian ρ sp We suppose that ρ sp ( ~r ) is either 0 everywhere in R \ V ,or at least negligibly small. If so, we can change integra-tion limits in (16b) to infinite.At this point we are interested in specific form of ρ sp (12), thus it is explicitly substituted into (16b). Be-sides we designate ~ K = 2 π~k T ˆ G and rewrite expression in Cartesian coordinates changing multiple integral to theproduct of integrals ρ ~k = ¯ ρ (2 πs ) / Y j =1 ∞ Z −∞ e − r j / (2 s ) − i K j r j dr j . Last expression can be integrated if we use followingrelation [49] + ∞ Z −∞ e − p x ± qx dx = √ πp e q / (2 p ) , ℜ (cid:0) p (cid:1) > . As result we will get ρ ~k = ¯ ρ Y j =1 e −K j s / dr j . Changing product of exponents to exponent of sum andmentioning that P j K j = 4 π ~k T ˆ G ˆ G T ~k we immediatelyget (18). D. On the expression of F int . Lets start with (15a). Inner integral has infinite bor-ders so we may rewrite this expression as follows F int = Z Z Z V Z Z Z R V ( | ~r ′ | ) ρ ( ~r ) ρ ( ~r + ~r ′ ) d ~r ′ d ~r. (46)Since ρ ( ~r ) is real it can be replaced with complex conju-gate ρ ∗ ( ~r ) without changing F int . Now we can substitute ρ from (16a) and mention that f ~k ( ~r + ~r ′ ) = f ~k ( ~r ) f ~k ( ~r ′ )(16c). So if we designate V ~k as in (17b) it can be written F int = 1¯ ρ X ~k ∈ Z ρ ~k V ~k X ~k ′ ∈ Z ρ ∗ ~k ′ Z Z Z V f ~k ( ~r ) f ∗ ~k ′ ( ~r ) d ~r. Last equation can be simplified if we perform integra-tion. From appendix VI B we expect orthogonality of f ~k functions. Let’s mention from (16c) f ~k ( ~r ) f ∗ ~k ′ ( ~r ) = f ~k − ~k ′ ( ~r ) and we will get Z Z Z V f ~k ( ~r ) f ∗ ~k ′ ( ~r ) d ~r = abc sin( α ) sin( β c ) δ ~k,~k ′ . (47)Here δ ~k,~k ′ = δ k x ,k ′ x δ k y ,k ′ y δ k z ,k ′ z — product of three Kro-necker’s delta.Mentioning that abc sin( α ) sin( β c ) = 1 / ¯ ρ , where ¯ ρ ismean particle density (one particle per V , see fig. 1 forgeometrical reasoning) we get (17a).2 E. V ~k for screened Coulomb potential Let’s consider equations (16c) and (16d) in sphericalcoordinates f ~k ( ~r ) = e π i r (cid:0) g ~k cos( θ )+ g ′ ~k cos( ϕ − δϕ ~k ) sin( θ ) (cid:1) ,g ~k = k x sin( α c − α ) a tan( β c ) sin( α ) − k y sin( α c ) b tan( β c ) sin( α ) + k z c sin( β c ) ,g ′ ~k = s k x a + (cid:18) k y b sin( α ) − k x cot( α ) a (cid:19) , and rewrite V ~k from (17b) V ~k = ¯ ρ ∞ Z drV ( r ) r π Z dθ sin( θ ) ×× π Z dϕe π i r (cid:0) g ~k cos( θ )+ g ′ ~k cos( ϕ − δϕ ~k ) sin( θ ) (cid:1)| {z } I ( θ ) . Integrating over ϕ we get I ( θ ) = 2 πe π i g ~k r cos( θ ) J (cid:16) πg ′ ~k r sin( θ ) (cid:17) . Then we substitute last expression into equation for V ~k V ~k = 2 π ¯ ρ ∞ Z V ( r ) r π Z e π i g ~k r cos( θ ) ×× J (cid:16) πg ′ ~k r sin( θ ) (cid:17) sin( θ ) dθdr. We will consider screened Coulomb potential (21),which means integration over r can be performed. Onemay use following relation [49] ∞ Z e − αx J ν ( βx ) x ν +1 dr = 2 α (2 β ) ν Γ (cid:0) ν + (cid:1) √ π ( α + β ) ν +3 / , ℜ ( α ) > |ℑ ( β ) | , ℜ ( ν ) > − V ~k to get V ~k = π Z π ¯ ρq (cid:0) /λ D − π i g ~k cos( θ ) (cid:1) sin( θ ) dθ (cid:18)(cid:0) /λ D − π i g ~k cos( θ ) (cid:1) + (cid:16) πg ′ ~k sin( θ ) (cid:17) (cid:19) / . Changing variable t = cos( θ ) and performing integrationover t we will get V ~k = 4 π ¯ ρq /λ D + 4 π g ~k + 4 π g ′ ~k . Last expression is equivalent to (22). One may check thisby expanding g ~k and g ′ ~k . F. Calculation of F s for screened Coulombpotential. We will start with equation (15b) and assume that ρ sp is a very “sharp” function (section I A). With this as-sumption we can change integration limits to infinite andthis, in turn, will allow us to perform variable exchangeas we did in (46) F s = Z Z Z R Z Z Z R V ( | ~r ′ | ) ρ sp ( ~r ) ρ sp ( ~r + ~r ′ ) d ~r ′ d ~r. (48)First of all we perform some mathematical transfor-mations with (12) (in Cartesian coordinate system) andwrite the following ρ sp ( ~r ) ρ sp ( ~r + ~r ′ ) = 18 π s e − [ ~r + ~r ′ / ] /s − ~r ′ / ( s ) . Now we can substitute this expression into (48) and per-form integration over ~r . Writing result in spherical co-ordinate system and integrating over angles we immedi-ately get (24).Since we know explicit expression for V (21) we cansubstitute it into (24) and perform integration using def-inition of the complementary error function [50]erfc( x ) = 2 √ π ∞ Z x e − t dt. (49)As result we will get F s = q √ πs (cid:18) − s √ πλ D e s /λ D erfc( s/λ D ) (cid:19) . (50)Since we may obviously suppose that s is very small com-pared to screening distance, then s ≪ λ D . So we canapproximate last equation and as result we will get (25). G. Approximation of F int − F s for screenedCoulomb potential We start with expression (23) for F int . But one maysee that it contains expressions s ˆ G ˆ G T . Components ofˆ G ˆ G T matrix in (23) are proportional to different productsof inverse distances in a lattice, e.g. 1 /a , 1 / ( ab ), and soon (see (16d)). Since s ≪ a , s ≪ b and s ≪ c (assump-tion about “sharp Gaussian”) we expect (cid:12)(cid:12)(cid:12) s ˆ G ˆ G T (cid:12)(cid:12)(cid:12) ≪ | ~k | . In this section we will try to rewriteequations to get even better convergence then we have.First of all we rewrite expression for F int as follows F int = 2 ¯ ρq π e s /λ D ∞ Z πs s e − s / (2 πλ D ) X ~k ∈ Z e − s ~k T ˆ G ˆ G T ~k d s . s takes values much smallerthen mean distance (26) as well as much bigger. Thuswe split this integral into two. First integration we willprovide from 2 πs to l and second from l to infinity.We assume that lattice is not too degenerated (sec-tion II B), which means that angles α and β c should notsignificantly differ from π/
2. So we may expect elementsof s ˆ G ˆ G T to be less than 1 for first integral and greaterthan 1 for second. If so, then second integration may beperformed and we end up with the sum that convergesmuch better than (23). F int = 2 ¯ ρq π e s /λ D l Z πs s e − s / (2 πλ D ) X ~k ∈ Z e − s ~k T ˆ G ˆ G T ~k | {z } Θ ( s ˆ G ˆ G T ) d s ++ 4 π ¯ ρq X ~k ∈ Z e s /λ D − l ( / [2 πλ D ] + ~k T ˆ G ˆ G T ~k )1 /λ D + 4 π ~k T ˆ G ˆ G T ~k . From (16d) one may check with Sylvester’s crite-rion [51] that ˆ G (16d) is a positive-definite matrix. Ob-viously, we expect ˆ G T and their product ˆ G ˆ G T to bepositive-definite matrices as well. Since s takes only pos-itive values s ˆ G ˆ G T is positive as well. This means high-lighted sum in the last equation is a Riemann theta func-tion [52] at point z = 0, so we appropriately designate itwith Θ (cid:16) s ˆ G ˆ G T (cid:17) .Using modular transformation [52]Θ( z ; ˆ A ) = π d/ p det ˆ A Θ( ˆ A − z ; ˆ A − ) , (51)where d is number of dimensions (3 in this case) we getΘ (cid:16) s ˆ G ˆ G T (cid:17) = π / s det ˆ G Θ (cid:16)
0; ˆ G − T ˆ G − / s (cid:17) . We used the fact that det ˆ G = det ˆ G T and that mul-tiplication by s is equal to multiplication by diagonalmatrix which has all elements equal to s . Explicit ex-pression for ˆ G − may be used from (40). Besides, wemay notice thatdet ˆ G = 1 abc sin( α ) sin( β c ) = ¯ ρ, (52)and rewrite expression for F int as follows F int = 2 √ πq e s /λ D X ~k ∈ Z l Z πs e − s / (2 πλ D ) − | ˆ G − ~k | / s d ss | {z } I ( ~k ) ++ 4 π ¯ ρq X ~k ∈ Z e s /λ D − l (cid:16) / [2 πλ D ] + | ˆ G T ~k | (cid:17) /λ D + 4 π ~k T ˆ G ˆ G T ~k . Now we can perform integration and find I (cid:16) ~k (cid:17) . Sinceresult is very complicated we first introduce two helperfunctions φ ( L ) = e − L / (2 πλ D ) L − √ πλ D erfc (cid:18) L πλ D (cid:19) ,f α ( L ) = e α | ˆ G − ~k | / ( πλ D ) erfc (cid:12)(cid:12)(cid:12) ˆ G − ~k (cid:12)(cid:12)(cid:12) L + αL πλ D , and present result in terms of this functions I (cid:16) ~k = 0 (cid:17) = φ (2 πs ) − φ ( l ) ,I (cid:16) ~k = 0 (cid:17) = √ π ( f − ( l ) − f − (2 πs )) / (cid:12)(cid:12)(cid:12) G − ~k (cid:12)(cid:12)(cid:12) + √ π ( f +1 ( l ) − f +1 (2 πs )) / (cid:12)(cid:12)(cid:12) G − ~k (cid:12)(cid:12)(cid:12) . Using precise expression for F s (50) one may see F s = 2 √ πq e s /λ D φ (2 πs ) . So subtracting F s from F int simply means neglection ofthis term.Our previous calculations do not rely on any approxi-mation, but expression for I (cid:16) ~k = 0 (cid:17) is very complicated,so we need to perform one. Before we do, one shouldshow it is acceptable to approximate this sum. Subse-quent calculations are divided into two parts: uniformconvergence proof and approximation itself. Proof of uniform convergence . We will use Cauchycriterion for sequence of functions f i ( x ) in domain E toachieve this goal: ∀ ε > ∃ N ∀ m ≥ n > N : ∀ x ∈ E : | P mi = n f i ( x ) | < ε . Further P ~k ∈ Z I (cid:16) ~k (cid:17) acts as asequence for convergence checking and entries of ˆ G − aswell as s , l and λ D are treated as variables from domain ofadmissible parameters. We have already examined prop-erties of ˆ G − and will not get into details again.First of all we use s < l < λ D . This relation holdsfor physical system under consideration and we will useit in all subsequent calculations. This means f α ( l ) >f α (2 πs ) and thus all terms in sum will be positive andwe can avoid using absolute value. The second thing thatappears is following bounding0 < X ~k ∈ Z \ ~ I (cid:16) ~k (cid:17) < X ~k ∈ Z \ ~ f − ( l ) + f +1 ( l ) (cid:12)(cid:12)(cid:12) G − ~k (cid:12)(cid:12)(cid:12) √ π. Then we use relation erfc( x ) ≤ e − x for x > f ± ( l ). Expanding squares in exponentsand performing simplification we get X ~k ∈ Z \ ~ I (cid:16) ~k (cid:17) < X ~k ∈ Z \ ~ exp (cid:18) − | ˆ G − ~k | l − l π λ D (cid:19)(cid:12)(cid:12)(cid:12) ˆ G − ~k (cid:12)(cid:12)(cid:12) . √ π < k ~x k =1 k ˆ A~x k = √ λ min , where λ min is minimal eigenvalue of ˆ A ∗ ˆ A and ∗ designates Hermitian adjoint [54]. We know thatdet ˆ G − = 1 / ¯ ρ > G − ∗ ˆ G − has nonzero eigen-values. Taking the minimal one λ min and designating g = √ λ min we claim (cid:12)(cid:12)(cid:12) ˆ G − ~k (cid:12)(cid:12)(cid:12) ≥ g (cid:12)(cid:12)(cid:12) ~k (cid:12)(cid:12)(cid:12) .Now we should change multidimensional summationto one-dimensional. Since all summands are positive wecan rearrange them in any convenient order. Thus we willperform summation over “cube surfaces”. First of all let’sdesignate K n = { ~k ∈ Z | | k x | ≤ n ∧ | k y | ≤ n ∧ | k z | ≤ n } .Then total sum can be expressed as X ~k ∈ Z \ ~ I (cid:16) ~k (cid:17) = + ∞ X n =1 X ~k ∈ K n \ K n − I (cid:16) ~k (cid:17) , where number of summands in K n \ K n − can be easilycalculated as (2 n + 1) − (2 n − ≡ n + 2. Minimallength for index vector in this set of indices is n .Collecting together previous calculations we get X ~k ∈ Z \ K m I (cid:16) ~k (cid:17) < + ∞ X n = m +1 (24 n + 2) exp (cid:16) − g n l − l π λ D (cid:17) gn . Further proof is supposed to be obvious and thus we avoidit. Proof of uniform convergence for first summand in F int can be performed much more easily by presentedscheme and thus is avoided as well. Instead one usefulremark shoud be done. Series converges better when g/l is bigger. Taking physical meaning of this variables intoaccount we get following statement: if we project vectors ~a , ~b and ~c on axis X , Y and Z respectively (see figure 1)and divide it by mean distance in lattice √ ¯ ρ , we willget measure of how good current series converges, or byphysical meaning — how close is current lattice to cubic. Approximation.
We take into account s ≪ l ≪ λ D ,neglect all terms containing l/λ D and s/λ D in erfc, andapproximate erfc( x → ∞ ) ∼ e − x / ( x √ π ) [55]. Obviously f α ( l ) ≫ f α (2 πs ) and thus all f with 2 πs arguments areneglected. Besides, we neglect all small terms in expo-nent and denominators. To make this procedure moreclear we state (cid:12)(cid:12)(cid:12) ˆ G − ~k (cid:12)(cid:12)(cid:12) /l > l/λ D . It means screeningdistance is big enough not to feel lattice deviation fromcubic. This leads us to equations φ ( l ) ≈ l ,I (cid:16) ~k = 0 (cid:17) ≈ l (cid:12)(cid:12)(cid:12) ˆ G − ~k (cid:12)(cid:12)(cid:12) e − | ˆ G − ~k | /l . To get the first one, identity erfc( x ) = 1 − erf( x ) andapproximation erf( x → ∼ x/ √ π [55] were used.Now approximate equation can be written (27). H. V ~k for capillary interaction Let’s consider two-dimensional version of (16c) and(16d) in polar coordinates f ~k ( ~r ) = e π i rg ~k cos( ϕ − δϕ ~k ) ,g ~k = s k x a + (cid:18) k y b sin( α ) − k x cot( α ) a (cid:19) , and rewrite V ~k from (17b) V ~k = ¯ ρ ∞ Z drV ( r ) r π Z dϕe π i rg ~k cos( ϕ − δϕ ~k ) | {z } I . Integrating over ϕ we get I = 2 π J (cid:0) πg ~k r (cid:1) . Then we substitute last equality and expression for V ( r )from (28) into equation for V ~k . It is known [49], if a > b > ∞ Z x J ( ax )K ( bx ) dx = 1 a + b . (53)With this equation we easily obtain V n,m for capillaryinteraction part V ~k = − q E σ ¯ ρ π g ~k + 1 /l . (54)Comparison with (29) shows that they are equal. I. Approximation of F int for capillary interaction We start with expression (30) for F int . First of all werewrite it as follows F int = − ¯ ρq E σ e s /l ∞ Z s s e − s /l X ~k ∈ Z e − π s ~k T ˆ G ˆ G T ~k d s . Then we change summation to integration P ~k ∈ Z → RR R . There is a method of integration for this kind offunctions [51]. As result we get F int = − ¯ ρq E σ e s /l ∞ Z s s e − s /l d s s det ˆ G . For two-dimension density it is true det ˆ G = ¯ ρ . Besides,we perform integration over s and get F int = − q E πσ e s /l Γ (cid:0) s /l (cid:1) , where Γ is gamma-function. Since w know that s ≪ l last equation can be approximated. We use knownapproximation [55] Γ ( x → ≈ − γ − ln( x ), where γ is Euler–Mascheroni constant. We neglect this constantand exponent. The rest can be written as (31).5 [1] Y. Bilotsky, Advances in Materials Science and Applica-tions, vol. 2, iss. 4, pp. 127–137, Dec. (2013).[2] L. N. Kantorovich and I. I. Tupitsyn, J. Phys.: Condens.Matter, “Coulomb potential inside a large finite crystal,”vol. 11, pp. 6159–6168.[3] L. P. Buhler, R. E. Crandal, J. Phys. A: Math. Gen.,“On the convergence problem for lattice sums,” vol. 23,pp. 2523–2528, (1990).[4] D. Borwein, J. M. Borwein, R. Shail and I. J. Zucker, “En-ergy of static electron lattices,” J. Phys. A: Math. Gen.,vol. 20, pp. 1519–1531, (1988).[5] V. E. Fortov et al., Phys. Rep. , 1 (2005).[6] H. Lowen, Phys. Rep. , 249 (1994).[7] G. E. Morfill, H. M. Thomas, U. Konopka, and M. Zuzic,Phys. Plasmas , 1769 (1999).[8] E. Helpand and F. H. Stilinger, J. Chem. Phys. , 1232(1968).[9] B. Klumov, G. Joyce, C. R¨ath, P. Huber et all, Struc-tural properties of 3D complex plasmas under micrograv-ity conditions , EPL, Vol.92, No.1, p.15003, (2010).[10] B. A. Klumov, G. E. Morfill,
Structural Properties ofComplex (Dusty) Plasma upon Crystallization and Melt-ing , JETP Letters, Vol.90, No.6, pp. 444–448, (2009).[11] B. I. Lev and A. G. Zagorodny, Phys. Lett. A , 158(2009).[12] B. I. Lev, V. B. Tymchyshyn, and A. G. Zagorodny, Con-dens. Phys. , 593 (2009).[13] H. Thomas, G. E. Morfill, V. Demmel, J. Goree, B. Feu-erbacher, and D. Mohlmann, Phys. Rev. Lett. , 652(1994).[14] J. H. Chu and Lin I, Phys. Rev. Lett. , 4009 (1994).[15] A. Melzer, T. Trottenberg, and A. Piel, Phys. Lett. A ,301 (1994).[16] S. V. Vladimirov, S. A. Khrapak, M. Chaudhuri, andG. E. Morfill, Phys. Rev. Lett. , 055002 (2008).[17] H. Ikezi, Phys. Fluids. , 1764 (1986).[18] A. Melzer, A. Homann, and A. Piel, Phys. Rev. E , 2757(1996).[19] A. G. Sitenko, A. G. Zagorodny, and V. N. Tsytovich, AIPConf. Proc. , 311 (1995).[20] S. A. Brazovsky, Sov. Phys. JETP , 715 (1975).[21] B. I. Lev and H. Yokoyama, IJMP B (International Jour-nal of Modern Physics) , 4913 (2003).[22] H. Totsuji, T. Kishimoto, and C. Totsuji, Phys. Rev. Lett. , 3113 (1997).[23] P. Leiderer Two-Dimensional Electron Systems by Ed.E.Y. Andrei, (Springer Netherlands 1997).[24] V. B. Shikin and P. Leiderer, Sov. Phys. JETP Lett., ,92, (1981).[25] D. C. Tsui, H. L. Stormer, and A. C. Gossard, Phys. Rev.Lett. , 1559 (1982).[26] R. B. Laughlin, Phys. Rev. Lett. , 1395 (1983).[27] V. S. Edelman, Uspekhi Fizicheskikh Nauk , 227(1980).[28] T. Ando, A. Fowler, and F. Stern, Rev. Mod. Phys. ,437 (1982).[29] W. Wigner, Phys. Rev. , 1002 (1934).[30] C. C. Grimes, and G. Adams, Phys. Rev. Lett. , 795 (1979).[31] P. M. Platzman and H. Fukuyama, Phys. Rev. B , 3150(1974).[32] I. Skachko, Phase diagram of a 2-dimensional electronsystem on the surface of liquid helium , Ph.D. thesis, Rut-gers, The State University of New Jersey (2006).[33] D. K. Lambert,
Electrons on the surface of liquid helium ,Ph.D. thesis, Lawrence Berkeley Laboratory (1979).[34] M. Haque, I. Paul, and S. Pankov, Phys. Rev. B ,045427 (2003).[35] B. I. Lev, V. P. Ostroukh, V. B. Tymchyshyn, andA. G. Zagorodny, Statistical description of the systemelectrons on the liquid helium surface , Eur. Phys. J. B,87:253 (2014).[36] L. P. Gor‘kov and D. M. Chernikova, JETP Lett. , 119(1973).[37] D. Ruelle, Statistical Mechanics: Rigorous Results ,W.A.Benjamin, New York, (1969).[38] A. Isihara,
Statistical Mechanics , Academic Press, NewYork-London, (1971).[39] K. Huang,
Statistical Mechanics , J. Wiley and Sons, NewYork, (1963).[40] R. Baxter,
Exactly Solved Models in Statistical Mechan-ics , Academic Press, New York, (1982).[41] Y. D. Bilotsky and B. I. Lev, Teor. Math. Fiz. , 120(1984).[42] B. I. Lev and A. Y. Zhugaevych, Phys. Rev. E, , 6460(1998).[43] B. I. Lev and A. G. Zagorodny, Statistical description ofCoulomb-like systems, Phys. Rev. E , 061115 (2011).[44] R. L. Stratonovich, Soviet Physics Doklady, Vol. 2,(1958).[45] J. Hubbard, Phys. Rev. Lett. , 77 (1959).[46] S. F. Edwards and A. Lennard, J. Math. Phys. , 778(1962).[47] S. Samuel, Phys. Rev. D , 1916 (1978).[48] A. Erd´elyi, W. Magnus, F. Oberhettinger, and F. G. Tri-comi, Higher Transcendental Functions , Vol. 1., McGraw-Hill, New York, (1953).[49] I. S. Gradshteyn and I. M. Ryzhik,
Table of integrals, se-ries, and products , Amsterdam: Elsevier / AcademicPress, (2007).[50] A. Erd´elyi, W. Magnus, F. Oberhettinger, and F. G. Tri-comi,
Higher Transcendental Functions , Vol. 2. Malabar,Krieger, (1981).[51] R. Bellman,
Introduction to Matrix Analysis , Society forIndustrial and Applied Mathematics, (1970).[52] R. Bellman,
A Brief Introduction to Theta Functions ,Holt, Rinehart and Winston, (1961).[53] M. Chiani, D. Dardari, M. K. Simon, New ExponentialBounds and Approximations for the Computation of Er-ror Probability in Fading Channels, IEEE Transactionson Wireless Communications, 4(2), pp 840–845, (2003).[54] C. Meyer,
Matrix Analysis and Applied Linear Algebra ,SIAM, (2000).[55] M. Abramowitz, I. Stegun,